The Runge -Kutta Method

As we' ve noted before,  
      Euler' s method is  a fairly intuitive but notterribly efficient means for obtaining numerical approximations 
        forthe solution of an initial value problem
                        
          y' (x)  =  f (x, y),    y ((Subscript[x, 0]))  =  (Subscript[y, 0])
\(The\ Runge - 
        Kutta\ method\ is\ a\ much\ more\ efficient\ 
          \((and\ much\ less\ 
intuitive)\)\ numerical\ method . \ \ The\ 
          following\ routine\ generates\ an\)approximate solution for the above initial value problem using
        Runge - 
        Kutta .    The syntax is the same as it was for Euler' s 
          methodin the demonstration from Section 8.1 .

\(RKMethod[f_, {x0_, y0_}, h_, n_]\  := 
	
\ \ \ \ \ 
      Module[
		\ \ \ \ \ \ \ {RKStep}, 
		\ \ \ \ \ \ \ \ 
        RKStep[{x_, y_}]\  := \ 
          Module[{k1, k2, k3, k4}, 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ k1\  = \ f[x, y]; 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            k2\  = \ f[x + h/2, \ y + \ k1*h/2]; 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            k3\  = \ f[x + h/2, y + \ k2*h/2]; 
										\ \ 
            k4\ \  = \ f[x + h, y + k3*h]; 

            		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            N[{x + h, y + h*\((k1\  + \ 2*k2\  + \ 2*k3\  + \ k4)\)/6}]

            			\ \ \ \ \ ]; 
		
		\ \ \ \ \ \ \ \ 
        Interpolation[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          NestList[RKStep, {x0, y0}, n], 
			\ \ \ \ \ \ \ \ \ \ \ \ 
          InterpolationOrder\  -> \ 1
			\ \ \ \ \ \ \ ]
	\ \ \ \ \ ]
        \)

As before,  
        we' ll try out our method on the initial value problem                  y' (x)  =  y (x), 
                y (0)  =  1and compare the results with the exact solution y 
              (x)  =  (e^x) .

f[x_, y_]  =  y

With a ridiculously large step size h=1, we obtain the following 
approximation of e.

approx  =  RKMethod[f, {0, 1}, 1, 1]

approx[1]

That's two digits of accuracy already!  If we cut the step size by a
factor of ten, we obtain

approx  =  RKMethod[f, {0, 1},  .1, 10]

approx[1]

\(TextForm
    \`which\ is\ accurate\ to\ \((at\ least)\)\ 6\ \(places!\)\ \ To\ obtain\ 
      higher\ accuracy, 

    we' d\ have\ to\ also\ specify\ more\ accuracy\ for\ the\ "N"\ command\ 
      used\)\(TextForm
    \`in\ defining\ RKMethod . \ \ \ It\ turns\ out\ that\ local\ truncation\ 
      error\ is\)\(TextForm
    \`proportional\ to\ \(h\^5\) for\ this\ method\ and\ global\ truncation\ 
      error\ is\)proportional to (h^4) .    That is,  
    for every factor of ten that we shrinkh,  we gain a factor of 10, 
    000  (i . e .  4 more decimal places)in accuracy .    Not bad .

In closing, let's graphically examine the difference 
between the exact and approximate solutions using h=.25 
and x ranging from 0 to 4.

approx  =  RKMethod[f, {0, 1},  .25, 16]

Plot[{approx[x],  Exp[x]}, {x, 0, 4}]

...virtually indistinguishable.

Plot[approx[x] - Exp[x], {x, 0, 4}]

Why do you suppose the difference oscillates so much?