The Runge-Kutta Methods For ODE

This notebook gives some examples for the Runge-Kutta method for solving ordinary differential equations. We start with Euler's method.

The RK2 method

The idea of RK2 is that we want to choose constants α, β, w1, and w2 so that a formula of the form 

	K1 = h f(t,x)
	K2 = h f(t + αh, x + β K1)
	x(t+h) = x(t) + w1 K1 + w2 K2
	
has the least error possible (in a Taylor series sense). By using the Taylor expansion of f, we computed in class that the system

	α = 1, β = 1, w1= w2 = 1/2 
	
will work. We go ahead and code this method

RK2[f_,x0_,t0_,t1_,n_] := Module[{α,β, w1,w2,h},

h = (t1 - t0)/(n-1);
α = 1;
β = 1;
w1 = 1/2;
w2 = 1/2;

RK2Step[{t_,x_,h_}] := Module[{K1,K2}, 
K1 = h*f[t,x];
K2 = h*f[t+α h, x + β K1];
{t+h,x + w1 K1 + w2 K2, h}
];

  NestList[RK2Step,{ t0,x0,h},n-1][[All,1;;2]]

];

Now we try applying this to the ODE:

	x'(t) = x(t),  t0 = 0, x(0) = 1. 
	
We already know the solution, but we try it with Mathematica just to see how it works:

ExactSoln = (x /. DSolve[{x'[t] == x[t],x[0] == 1},x,t])[[1]]

We now need to use our RK2 code. For this equation, we write

f[t_,x_] := x;
RK2[f,1.0,0.0,1.0,10]

This is easier to understand in a plot:

Manipulate[Show[Plot[ExactSoln[t],{t,0,1}],ListLinePlot[RK2[f,1.0,0.0,1.0,n],PlotMarkers->Automatic],

AxesLabel->{Style["t",Larger,FontFamily->"Helvetica",Bold],Style["x",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Exact and Runge-Kutta 2 Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")],

{{n,10},2,100,1}]

We can go ahead and look at just the result of the calculation as well:

RK2Val[f_,x0_,t0_,t1_,n_] := Module[{α,β, w1,w2,h},

h = (t1 - t0)/(n-1);
α = 1;
β = 1;
w1 = 1/2;
w2 = 1/2;

RK2Step[{t_,x_,h_}] := Module[{K1,K2}, 
K1 = h*f[t,x];
K2 = h*f[t+α h, x + β K1];
{t+h,x + w1 K1 + w2 K2, h}
];

  Nest[RK2Step,{ t0,x0,h},n-1][[2]]

];

RK2Val[f,1.0,0.0,1.0,100]

ExactSoln[1.0]

As you can see, even for this comparatively trivial case, there aren't a lot of correct digits in the solution here. We go ahead and make the plot of correct digits just to see how fast this method is converging to a correct solution.

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];
ListLinePlot[Table[{n,Digits[RK2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct Digits in RK2 Solution of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

We now compare this to the old Taylor2 method just to see that we're really getting the same kind of performance, even though we never explicitly differentiated f.

Taylor2Val[f_,x0_,t0_,t1_,n_] := Module[{Taylor2Step,h,ft,fx,t,x,u,v},

ft[t_,x_]:= D[f[u,v],u] /. {u->t,v->x};
fx[t_,x_]:= D[f[u,v],v] /. {u->t,v->x};

h = (t1 - t0)/(n-1);
Taylor2Step[{t_,x_,h_}] := {t+h,x +h*f[t,x] + (h^2/2)*(ft[t,x] + fx[t,x]*f[t,x]),h};
Nest[Taylor2Step,{ t0,x0,h},n-1][[2]]

];

ListLinePlot[Table[{n,Digits[Taylor2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct digits in Taylor2 solution of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

Since these graphs look almost the same (indeed, when you plot them on top of one another, you can't see a difference), we plot the difference between the number of correct digits in the two methods to see how we do.

ListLinePlot[Table[{n,Digits[RK2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]] - Digits[Taylor2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Difference between # of correct digits in RK2 and Taylor2 solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

The Runge-Kutta 4th order method.

Here we are interested in taking steps in the form:

	x(t+h) = x(t) + w1 K1 + w2 K2 + w3 K3 + w4 K4
	
where

	K1 = h f(t,x)
	K2 = h f(t + α h, x + β K1)
	K3 = h f(t + γ h, x + δ K2) 
	K4 = h f(t + ϵ h, x + ϕ K3)
	
It turns out that 

	α = β = γ = δ = 1/2, ϵ = ϕ = 1, w1 = w4 = 1/6, w2 = w3 = 1/3
	
agrees with our Taylor expansion of x(t+h) through the 4th order. We code this method:

RK4[f_,x0_,t0_,t1_,n_] := Module[{α,β, γ, δ, ϵ, ϕ, w1,w2,w3, w4, h},

h = (t1 - t0)/(n-1);

{α,β,γ, δ, ϵ, ϕ} = {1/2,1/2,1/2,1/2,1,1};
{w1,w2,w3,w4} = {1/6,1/3,1/3,1/6};

RK4Step[{t_,x_,h_}] := Module[{K1,K2, K3, K4}, 

K1 = h*f[t,x];
K2 = h*f[t+α h, x + β K1];
K3 = h*f[t+γ h, x + δ K2];
K4 = h*f[t+ϵ h, x + ϕ K3];

{t+h,x + w1 K1 + w2 K2 + w3 K3 + w4 K4, h}

];

  NestList[RK4Step,{ t0,x0,h},n-1][[All,1;;2]]

];

RK4[f,1.0,0.0,1.0,10]

As before, we'll show this in a plot, this time versus the RK2 solution

Manipulate[Show[Plot[ExactSoln[t],{t,0,1}],ListLinePlot[RK4[f,1.0,0.0,1.0,n],PlotMarkers->Automatic,PlotStyle->{Red}],
ListLinePlot[RK2[f,1.0,0.0,1.0,n],PlotMarkers->Automatic],
AxesLabel->{Style["t",Larger,FontFamily->"Helvetica",Bold],Style["x",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Exact, RK2, and RK4 (red) Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")],
{{n,10},2,100,1}]

We can also check the values for the solver:

RK4Val[f_,x0_,t0_,t1_,n_] := Module[{α,β, γ, δ, ϵ, ϕ, w1,w2,w3, w4, h},

h = (t1 - t0)/(n-1);

{α,β,γ, δ, ϵ, ϕ} = {1/2,1/2,1/2,1/2,1,1};
{w1,w2,w3,w4} = {1/6,1/3,1/3,1/6};

RK4Step[{t_,x_,h_}] := Module[{K1,K2, K3, K4}, 

K1 = h*f[t,x];
K2 = h*f[t+α h, x + β K1];
K3 = h*f[t+γ h, x + δ K2];
K4 = h*f[t+ϵ h, x + ϕ K3];

{t+h,x + w1 K1 + w2 K2 + w3 K3 + w4 K4, h}

];

  Nest[RK4Step,{ t0,x0,h},n-1][[2]]

];

RK4Val[f,1.0,0.0,1.0,10] // InputForm

ExactSoln[1.0] // InputForm

We'll now make a plot of the correct digits in the RK2 amd RK4 methods:

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];
Show[ListLinePlot[{Table[{n,Digits[RK4Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],
Table[{n,Digits[RK2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}]},
Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct Digits in RK2 and RK4 Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]]