RK-4 General Routine, without step-size control

Andy Long
MAT360
Spring, 2015

First of all, notice how beautifully simple this is. We’re not using any step-size control, and that might make it better -- but this is a pretty good routine for something so simple!

This version of RK-4 with fixed step-size is for general non-autonomous f(t,y):

rk4[f_,init_,tinit_,tfinal_,h_]:=
Module[
{n=Ceiling[(tfinal-tinit)/h],halfh=0.5h,t=tinit,w=init,wlist=init,g1,g2,g3,g4},
{wlist} ~Join~
Table[
g1=f[t,w];
g2=f[t+halfh,w+halfh g1];
g3=f[t+halfh,w+halfh g2];
g4=f[t+h,w+h g3];
t=t+h; (* update t for the next step *)
w=w+h(g1+2. g2+2. g3+g4)/6.
,
{n}
]
]
Clear[x,y,z,t,f]
f[t_,y_]:= t*t+y*y
rk4[f,1,2,2.2,0.1]

Once again, it’s easy to operate on systems. Here again is the Lorenz System, exhibiting chaos:

Clear[x,y,z,f]
f[t_,{x_,y_,z_}]:={
-3.(x-y),
-x z+26.5 x-y,
x y-z
}
(* We just give an initial value that is also a list, of the proper size: *)
rk4[f,{0.,1.,0.},0,20.,0.0004];
Show[Graphics3D[{Line[%]}],Axes->Automatic]