Section 6.3: Foundations for RK-4

Andy Long, Spring, 2024

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

The example is from our text, Tea Time Numerical Analysis

This version of RK-4 with fixed step-size is for general non-autonomous f(t,y): that means you can incorporate time into your solutions.

We’ll create the routines, and use the “standard example” to test the various routines (comparing to tables in the text):

Clear[x,y,z,t,f]
f[t_,y_]:= -y/t+t^2
yExact[t_]:=t^3/4+16/t
y0=20.0;
t0=4;
tf=2;
h=-0.25;
pExact=Plot[yExact[t],{t,t0,tf}];

Euler-ode

rkFoundationsEuler[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];
t=t+h; (* update t for the next step *)
w=w+h g1,
{n}
]
]

solnEuler=rkFoundationsEuler[f,y0,t0,tf,h]
pEuler=ListPlot[Table[{t0+(n-1)*h,solnEuler[[n]]},{n,1,Length[solnEuler]}],PlotStyle->Blue];
Show[pEuler,pExact,Axes->Automatic]

Trap-ode

rkFoundationsTrap[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+h,w+h g1];
t=t+h; (* update t for the next step *)
w=w+halfh (g1 + g2),
{n}
]
]

solnTrap=rkFoundationsTrap[f,y0,t0,tf,h]
pTrap=ListPlot[Table[{t0+(n-1)*h,solnTrap[[n]]},{n,1,Length[solnTrap]}],PlotStyle->Red];
Show[pTrap,pExact,Axes->Automatic]

Mid-ode

rkFoundationsMid[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];
t=t+h; (* update t for the next step *)
w=w+h g2,
{n}
]
]

solnMid=rkFoundationsMid[f,y0,t0,tf,h]
pMid=ListPlot[Table[{t0+(n-1)*h,solnMid[[n]]},{n,1,Length[solnMid]}],PlotStyle->Black];
Show[pMid,pExact,Axes->Automatic]

Simp-ode

rkFoundationsSimp[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+h,w+h g2];
t=t+h; (* update t for the next step *)
w=w+h (g1 + 4g2+g3)/6,
{n}
]
]

solnSimp=rkFoundationsSimp[f,y0,t0,tf,h]
pSimp=ListPlot[Table[{t0+(n-1)*h,solnSimp[[n]]},{n,1,Length[solnSimp]}],PlotStyle->Black];
Show[pSimp,pExact,Axes->Automatic]

solnExact=Table[yExact[t0+k h],{k,0,Length[solnEuler]-1}]
TableForm[TableForm[Transpose[{solnEuler,solnTrap,solnMid,solnSimp,solnExact}],TableHeadings->{None,{"Euler-ode","Trap-ode","Mid-ode","Simp-ode","Exact"}}]]

Open and ClOpen methods

rkFoundationsOpen[f_,init_,tinit_,tfinal_,h_]:=
Module[
{n=Ceiling[(tfinal-tinit)/h],halfh=h/2.0,thirdh=h/3.0,twothirdh=2*h/3.0,t=tinit,w=init,wlist=init,g1,g2,g3,g4},
{wlist} ~Join~
Table[
g1=f[t,w];
g2=f[t+thirdh,w+thirdh g1];
g3=f[t+twothirdh,w+twothirdh g2];
t=t+h; (* update t for the next step *)
w=w+halfh (g2+g3),
{n}
]
]

solnOpen=rkFoundationsOpen[f,y0,t0,tf,h]
pOpen=ListPlot[Table[{t0+(n-1)*h,solnOpen[[n]]},{n,1,Length[solnOpen]}],PlotStyle->Black];
Show[pOpen,pExact,Axes->Automatic]

rkFoundationsClosed[f_,init_,tinit_,tfinal_,h_]:=
Module[
{n=Ceiling[(tfinal-tinit)/h],fourthh=h/4.0,thirdh=h/3.0,twothirdh=2*h/3.0,t=tinit,w=init,wlist=init,g1,g2,g3,g4},
{wlist} ~Join~
Table[
g1=f[t,w];
g2=f[t+thirdh,w+thirdh g1];
g3=f[t+twothirdh,w+twothirdh g2];
t=t+h; (* update t for the next step *)
w=w+fourthh (g1+3g3),
{n}
]
]

solnClosed=rkFoundationsClosed[f,y0,t0,tf,h]
pClosed=ListPlot[Table[{t0+(n-1)*h,solnClosed[[n]]},{n,1,Length[solnClosed]}],PlotStyle->Black];
Show[pClosed,pExact,Axes->Automatic]

RK4-ode

Let’s add in Runge-Kutta-4: it’s our target, and it certainly looks like these other methods!

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}
]
]

solnRK4=rk4[f,y0,t0,tf,h]
pRK4=ListPlot[Table[{t0+(n-1)*h,solnRK4[[n]]},{n,1,Length[solnRK4]}],PlotStyle->Black];
Show[pRK4,pExact,Axes->Automatic]

TableForm[TableForm[Transpose[{solnSimp,solnOpen,solnClosed,solnRK4,solnExact}],TableHeadings->{None,{"Simp-ode","Open-ode","ClOpen-ode","RK4-ode","Exact"}}]]

It’s easy to operate on systems:

Here’s the pendulum problem:

Clear[x,y,z,f,c,m,g,l]
(* We're thinking of f as a vector-valued function, obtained by reducing the second order equation for the pendulum to a first order system. The arguments a time t and a location (u,v): *)   
f[t_,{u_,v_}]:={
-c/m u-g/l Sin[v],
u
}
g=9.81; (* m/s2 *)
c=1.1;
l=.31; (* m *)
m=70; (* kg *)
theta0=Pi/3;
thetaprime0=0;

(* We just give an initial value that is also a list, of the proper size: *)
rk4[f,{thetaprime0,theta0},0,10,0.0025];
Show[Graphics[{Line[%]}],Axes->Automatic]

Here is the lovely 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
}
rk4[f,{0.,1.,0.},0,50.,0.0004];
Show[Graphics3D[{Line[%]}],Axes->Automatic]