In 'Numerical methods for engineers' by Chapra and Canale there was an example in order to prove that Runge-Kutta-Fehlberg method is not always accurate, where given the equation $$\frac{dy}{dx}+0.6y=10e^{-(x-2)^2/[2(0.075)^2]}$$ with initial condition $y(0)=0.5$ the general solution is said to be $$y=0.5e^{-0.6x}$$

That way it is true (as easily seen by the graphs) that the approximate and the true solution are not in close agreement, but why is the error function term ignored by the general solution? If we take in consideration this term then the two solutions are not so different.


Look at a graph of the right side. Or consider that its exponent is $$ -\frac{(x-2)^2}{2\left(\frac3{40}\right)^2}=-\frac{1600(x-2)^2}{18} $$ and $\ln(10^{-17})=-39.14..$ which level the expression will reach for $|x-2|>0.66..$.

Thus from a numerical perspective, the right side does not exist before $x=1.3$ and after $x=2.7$, and on the segment $[0,1.3]$ the solution to the now homogeneous equation is indeed $0.5·e^{-0.6x}$.

Using classical RK4 estimating the error by repeating steps with half the step size (or more precisely, for every two steps one step for the same time interval using the double step size) I get the following picture for a relative error threshold of $10^{-4},10^{-5},10^{-6}$.

RK4 approximate and more exact solutions

with a computation log for the $10^{-4}$ integration of

new step size    h=  0.06666667  at t=  0.00000000       with fac= 10.0000
internal step to t=  0.06666667, y[0]=  0.48039472,      abs. error =  1.99859e-14      rel. error =  2.99788e-13       curve length =     0.066667 
new step size    h=  0.66666667  at t=  0.06666667       with fac= 10.0000
internal step to t=  0.73333333, y[0]=  0.32202024,      abs. error =  1.91098e-09      rel. error =  2.60588e-09       curve length =     0.733333 
internal step to t=  1.40000000, y[0]=  0.21585798,      abs. error =  1.21963e-05      rel. error =  8.71164e-06       curve length =     1.400000 
new step size    h=  0.22417033  at t=  1.40000000       with fac=  0.3363
internal step to t=  1.62417033, y[0]=  0.18869247,      abs. error =   0.00176841      rel. error =   0.00108881       curve length =     1.624170 
new step size    h=  1.15353919  at t=  1.62417033       with fac=  5.1458
internal step to t=  2.77770952, y[0]=  0.09701207,      abs. error =    0.0038085      rel. error =   0.00137109       curve length =     2.777710 
internal step to t=  3.93124872, y[0]=  0.04856108,      abs. error =   0.00825655      rel. error =   0.00210023       curve length =     3.931249 
internal step to t=  5.08478791, y[0]=  0.02430809,      abs. error =    0.0178082      rel. error =   0.00350225       curve length =     5.084788 
internal step to t=  6.23832710, y[0]=  0.01216784,      abs. error =    0.0383644      rel. error =   0.00614979       curve length =     6.238327 
internal step to t=  7.39186630, y[0]=  0.00609082,      abs. error =    0.0826262      rel. error =     0.011178       curve length =     7.391866 
internal step to t=  8.54540549, y[0]=  0.00304887,      abs. error =     0.177942      rel. error =    0.0208231       curve length =     8.545405 
internal step to t=  9.69894468, y[0]=  0.00152616,      abs. error =     0.383207      rel. error =    0.0395102       curve length =     9.698945 
new step size    h=  3.92357535  at t=  9.69894468       with fac=  3.4013
internal step to t= 13.62252003, y[0]=  0.00016007,      abs. error =      1.88675      rel. error =     0.138502       curve length =    13.622520 

where one sees that at $x=1.62417033$, directly before the jump, the step size is increased to $dx=1.15353919$, thus neatly stepping over the bump even in the half step evaluations.

But this is more artefact of the implementation and a rather low error threshold. Already tightening the relative error margin to $10^{-6}$ corrects the solution to be as good as the best numerical solution.

  • $\begingroup$ I forgot to mention that the segment reffered in the example is [0,4] . So the right side that exists after 1.3 and before 2.7 as you said, is totally ignored $\endgroup$ – asdf Dec 5 '16 at 13:00
  • $\begingroup$ One can of course not ignore that substantial bump. What are the absolute and/or relative error margins used in the step size computations? I could imagine that RKF in the segment $[0,1]$ computes an optimal step size $dx$ that is larger than $12$ and thus completely misses the bump also in the internal steps that start at $dx/4$. $\endgroup$ – LutzL Dec 5 '16 at 13:04
  • $\begingroup$ The stepsize is said to be 0.1 but I dont think that ios a problem.As stated in the book : the general solution is $y=0.5e^{-0.6x}$ which is a smooth curve that gradually approaches zero as x increases. In contrast, the particular solution undergoes an abrupt transition in the vicinity of x = 2 due to the nature of the forcing function $\endgroup$ – asdf Dec 5 '16 at 13:22
  • $\begingroup$ oh now I see.. I misunderstood, I thought that the 'particular solution' reffered to the numerical one we get by applying the method. But I guess this reffers to the general solution, right? $\endgroup$ – asdf Dec 5 '16 at 13:24
  • $\begingroup$ A particular solution in the context of linear ODE is any solution of the inhomogeneous equation. The full solution is then the general homogeneous solution plus the fixed particular solution. $\endgroup$ – LutzL Dec 5 '16 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.