RK-4 Examples

Andy Long
Mat360
Spring, 2017

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

rk4[f_,init_,tfinal_,h_]:=
Module[
{n=Ceiling[tfinal/h],t=0.,w=init,g1,g2,g3,g4},
{w} ~Join~ 
Table[g1=f[w];
g2=f[w+0.5 h g1];
g3=f[w+0.5 h g2];
g4=f[w+h g3];
(* The g? are derivatives, and there are four of them. We can fit a cubic to those, and then get a fourth derivative. The we use the same step-size control that we used on the Taylor method. How does that sound? 
*)
w=w+h(g1+2. g2+2. g3+g4)/6.
,
{n}
]
]

Rabbits and foxes:

Clear[r,f,g,solution]
rprime[r_,f_]:=2. r-0.01 r f
fprime[r_,f_]:=-f+0.01 r f
g[{r_,f_}]:={rprime[r,f],fprime[r,f]}
solution=rk4[g,{300,150},20.,0.04];
ListPlot[solution[[Range[Length[solution]],{1,2}]],AxesLabel->{"time","rabbits"},Joined->True]
(*part b*)
ListPlot[solution[[Range[Length[solution]],{1,2}]],PlotLabel->"trajectory",AxesLabel->{"rabbits","foxes"}]

The Lorenz System:

f[{x_,y_,z_}]:={-3.(x-y),-x z+26.5 x-y,x y-z}
rk4[f,{0.,1.,0.},20.,0.0004];
Show[Graphics3D[{Line[%]}],Axes->Automatic]

The exponential Problem:

f[x_]:=x
n=10;
h=1/n;
exp=rk4[f,1.,1.,h];
p1=ListPlot[Table[{i*h,Exp[i*h]-exp[[i+1]]},{i,n}]]
p1=ListPlot[Table[{i*h,exp[[i+1]]},{i,n}]];
p2=Plot[Exp[x],{x,0,1}];
Show[p1,p2]
error=exp[[n+1]]-Exp[1]

RK-4 with Step-Size Control -- an attempt!

dividedDifference[list_]:=
Module[
{n=Length[list]},
If[n==1,
b=list[[1]][[2]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
b
]

RK4WithStepSizeControl[f_,vi_,xf_]:=
Module[
{xs,vs,halfh,ltol=0.0001,maxsteps=200,k=0,h,xnew},
h=xf/maxsteps*ltol;
xs={x=0.};
vs={v=vi};
While[x<xf&&k<maxsteps, (* Don't go further than xf, and we've not taken more than the allowed number of step: *)
halfh=h/2;
v1=f[x,v];
v2=f[x+halfh,v+halfh v1];
v3=f[x+halfh,v+halfh v2];
v4=f[x+h,v+h v3];
(* The v? are derivatives, and there are four of them. We can fit a cubic to those, and then get a fourth derivative. Then we use the same adaptive step-size control that we used on the Taylor method. How does that sound? 
*)
apderiv=6*dividedDifference[{{0,v1},{h/3,v2},{2h/3,v3},{h,v4}}];
h=(4! ltol/(Abs[apderiv]))^0.25;
(* Print[h," ",apderiv]; *)
(* check that we don't step past the final value of x: *)
xnew=Min[x+h,xf];
(* adjust h if it happens to be too big, stepping beyond the final x value *)
h=xnew-x;
halfh=h/2;
v1=f[x,v];
v2=f[x+halfh,v+halfh v1];
v3=f[x+halfh,v+halfh v2];
v4=f[x+h,v+h v3];
v=v+h(v1+2.v2+2. v3+v4)/6.;
AppendTo[xs,x];
AppendTo[vs,v];
x=xnew;
k=k+1];
{xs,vs}]

Some physical parameters: radius of the Earth r, and the acceleration due to gravity g

(* here we're starting with the escape velocity *)
r=6378.388; (* in km *)
g=0.00980621;  (* in km/s^2 *)
vi=Sqrt[2. g r] (* in km/s *)
xf=10^6;
f[x_,v_]:=-g r^2/(v (r+x)^2)
f[0,vi]

t4=RK4WithStepSizeControl[f,vi,xf] ;
(* This equation can be solved by separation of variables: *)
vTheoretical[x_]:=Sqrt[2g*r^2/(r+x) +vi^2-2g*r]
p0=Plot[vTheoretical[x],{x,0,xf},PlotRange->{{0,xf},{0,vi}}];
p1=ListPlot[Transpose[t4],PlotMarkers->{Automatic,Medium}]
Show[p0,p1]
(*The foregoing computation is a striking illustration of the power of automatic stepsize control.*)

https://math.okstate.edu/people/yqwang/teaching/.../rungekutta.pdf suggests we should see

function  v  =  f(t,y)
v  =  y-tˆ2+1;
Run this program, we see that starting with h = 0.2, the program outputs
Step  0:  t  =  0.0000,  w  =    0.500000000000000
Step  1:  t  =  0.2000,  w  =    0.829299076923077
Step  2:  t  =  0.4353,  w  =    1.287432405787216
Step  3:  t  =  0.6766,  w  =    1.827289794651997
Step  4:  t  =  0.9264,  w  =    2.448301479233138
Step  5:  t  =  1.1902,  w  =    3.153049280338359
Step  6:  t  =  1.4806,  w  =    3.955581050460808
Step  7:  t  =  1.8537,  w  =    4.952039512278185
Step  8:  t  =  2.0000,  w  =    5.305486816572746

RKF45WithStepSizeControl[f_,vi_,xf_,step_,epsilon_]:=
Module[
{x,v,xs,vs,vone, vstar, halfh,fourth,delta,R,k1,k2,k3,k4,k5,k6,maxsteps=200,k=0,h=step,xnew},
xs={x=0.};
vs={v=vi};
Print["Step ",k,": t=",x,", w=", v];
While[x<xf&&k<maxsteps, (* Don't go further than xf, and we've not taken more than the allowed number of step: *)
xnew=Min[x+h,xf];
h=xnew-x;
halfh=h/2;
fourth=h/4;
k1=h*f[x,v];
k2=h*f[x+fourth,v+k1/4];
k3=h*f[x+3h/8,v+3/32* (k1+3* k2)];
k4=h*f[x+12h/13,v+(1932*k1-7200*k2+7296*k3)/2197];
k5=h*f[x+h,v+(439*k1/216-8*k2+3680*k3/513-845*k4/4104)];
k6=h*f[x+halfh,v+(-8*k1/27+2*k2+-3544*k3/2565+1859*k4/4104-11*k5/40)];
vone=v+25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
vstar=v+16*k1/135+6656*k3/12825+28561*k4/56430-9k5/50+2*k6/55;
R=Abs[vstar-vone]/h;
delta=0.84(epsilon/R)^0.25;
h=delta*h;
If[R≤epsilon,
v=vone;
AppendTo[xs,x];
AppendTo[vs,v];
x=xnew;
Print["Step ",k,": t=",x,", w=", v];
k=k+1,
Print["Rejected h: ",h/delta];
]
];
{xs,vs}]

f[t_,y_]:=y-t^2+1
RKF45WithStepSizeControl[f,0.5,2,0.2,.00001]

r=6378.388; (* in km *)
g=0.00980621;  (* in km/s^2 *)
vi=Sqrt[2. g r] (* in km/s *)
f[x_,v_]:=-g r^2/(v (r+x)^2)
xf=10^6;
t4=RKF45WithStepSizeControl[f,vi,xf,100,10^-8];
(* This equation can be solved by separation of variables: *)
vTheoretical[x_]:=Sqrt[2g*r^2/(r+x) +vi^2-2g*r]
p0=Plot[vTheoretical[x],{x,0,xf},PlotRange->{{0,xf},{0,vi}}];
p1=ListPlot[Transpose[t4],PlotMarkers->{Automatic,Medium}]
Show[p0,p1]
(*The foregoing computation is a striking illustration of the power of automatic stepsize control.*)