Vector RK Methods

The sweet thing about Mathematica is that we don't even have to really change the code...

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

];

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

];

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:

RK4Step[{f_,{α_,β_,γ_, δ_, ϵ_, ϕ_},{w1_,w2_,w3_,w4_},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];

{f,{α,β,γ, δ, ϵ, ϕ},{w1,w2,w3,w4},t+h,x + w1 K1 + w2 K2 + w3 K3 + w4 K4, h}

];


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


  NestList[RK4Step,{f,{α,β,γ, δ, ϵ, ϕ} ,{w1,w2,w3,w4}, t0,x0,h},n-1][[All,4;;5]]

];

We will now try the system from class notes:

x'(t) = y(t)
y'(t) = 1 - x(t) + 4/x(t)^3

F[t_,x_]  := {x[[2]],1 - x[[1]] + 4/(x[[1]]^3)};

RK4[F,{1.0,1.0},0.0,1.0,10]

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

Manipulate[Show[ListLinePlot[RK4[F,{1.0,1.0},0.0,1.0,n][[All,2]],PlotMarkers->Automatic,PlotStyle->{Red}],
ListLinePlot[RK2[F,{1.0,1.0},0.0,1.0,n][[All,2]],PlotMarkers->Automatic],
AxesLabel->{Style["t",Larger,FontFamily->"Helvetica",Bold],Style["x",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["RK2 and RK4 (red) Solutions of System",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]]

];