Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Oct 2016
Lecture topic: Runge Kutta

Taylor Polynomial of Two Variables

Here is the general Taylor polynomial in two variables:

Clear[f,t0,y0,TaylorSeries];
TaylorSeries[n_,t0_,y0_,t_,y_]:=(∑_(k=0)^n)(1/(k!))(∑_(j=0)^k)Binomial[k,j](((t-t0))^(k-j))(((y-y0))^j)Derivative[k-j,j][f][t0,y0]

Here are the first four expansions:

TaylorSeries[0,t0,y0,t,y]
TaylorSeries[1,t0,y0,t,y]
TaylorSeries[2,t0,y0,t,y]
TaylorSeries[3,t0,y0,t,y]

If we define a function, we can then use this on that function.

f[t_,y_]=-Exp[-((((t-2))^2)/4)-((((y-3))^2)/4)]Cos[2t+y-7];

TaylorSeries[2,2,3,t,y]

Here is a plot of the function:

plotfcn=Plot3D[f[t,y],{t,0,5},{y,1,5},PlotPoints->40]

Here is a plot of the approximation:

plotseries=Plot3D[TaylorSeries[2,2,3,t,y],{t,0,5},{y,1,5},PlotPoints->40,Mesh->False,ColorFunction->"BlueGreenYellow",PlotStyle->Opacity[0.6]]

Show[plotfcn,plotseries,PlotRange->{Automatic,Automatic,{-1,1}},AxesLabel->{"x","y","z"}]

We see that the approximation is close to the surface in the region near (x,y)=(2,3), which is where we based our expansion around. Higher order terms will produce better approximations.

plotseries=Plot3D[TaylorSeries[6,2,3,t,y],{t,0,5},{y,1,5},PlotPoints->40,Mesh->False,ColorFunction->"BlueGreenYellow",PlotStyle->Opacity[0.6]];
Show[plotfcn,plotseries,PlotRange->{Automatic,Automatic,{-1,1}},AxesLabel->{"x","y","z"}]

Deriving Runge Kutta Methods of Order 2: MidPoint, Trapezoidal, Heun’s

Here is Taylor's method of order 2:

Clear[f,h,a1,a2,k,α,β]
TaylorMethod[t0_,y0_,n_]:=Module[{t,y,data,TDerf},

(Subscript[TDerf, 0])[t_,y_]=f[t,y];
t[0]=t0;
y[0]=y0;
Do[
(Subscript[TDerf, k])[t_,y_]=D[(Subscript[TDerf, k-1])[t,y],{t,1}]+D[(Subscript[TDerf, k-1])[t,y],{y,1}]*y'[t]/.y'[t]->f[t,y]
,{k,1,n-1}];

t[i_]:=t[i]=t[i-1]+h;
y[i_]:=y[i]=y[i-1]+(∑_(k=1)^n)((h^k)/(k!))*(Subscript[TDerf, k-1])[t[i-1],y[i-1]];
Return[{t[1],y[1]}];
]

rel1=TaylorMethod[t0,y0,2][[2]]

We want to get rid of the derivatives of f in this expression. This time, let's assume the following, where a1, a2, α and β are to be determined:

k[1]=h*f[t0,y0]
k[2]=h*f[t0+h*α[2],y0+β[2]*k[1]]

Here is the approximation to the new y coordinate:

rel2=y0+a1*k[1]+a2*k[2]

We need to expand this in powers of h:

rel2=Series[rel2,{h,0,2}]

Now, we compare coefficients of powers of h between re1 and rel2:

eq1=Coefficient[rel1,h f[t0,y0] ]==Coefficient[rel2,h  f[t0,y0]]
eq2=Coefficient[rel1,(h^2)Derivative[1,0][f][t0,y0]]==Coefficient[rel2,(h^2)Derivative[1,0][f][t0,y0]]
eq3=Coefficient[rel1,(h^2)f[t0,y0] Derivative[0,1][f][t0,y0]]==Coefficient[rel2,(h^2)f[t0,y0] Derivative[0,1][f][t0,y0]]

Now, solve the system of equations (it is underdetermined, so we can solve for everything in terms of a2):

vals=Solve[{eq1,eq2,eq3},{a1,α[2],β[2]}]

And finally, the Runga Kutta Order 2 result is obtained:

k[1]/.Flatten[vals]
k[2]/.Flatten[vals]
y0+a1*k[1]+a2*k[2]/.Flatten[vals]

Doing it this way, we have actually recovered a bunch of different Runga Kutta Order 2 results, which we get by changing the value of a2.

If a2=1, we recover the Midpoint Method:

k[1]/.Flatten[vals]/.a2->1
k[2]/.Flatten[vals]/.a2->1
y0+a1*k[1]+a2*k[2]/.Flatten[vals]/.a2->1

If a2=1/2, we recover the Explicit Trapezoidal Method:

k[1]/.Flatten[vals]/.a2->1/2
k[2]/.Flatten[vals]/.a2->1/2
y0+a1*k[1]+a2*k[2]/.Flatten[vals]/.a2->1/2

If a2=3/4, we get Heun’s Method:

k[1]/.Flatten[vals]/.a2->3/4
k[2]/.Flatten[vals]/.a2->3/4
y0+a1*k[1]+a2*k[2]/.Flatten[vals]/.a2->3/4

Runge Kutta

Here is an example of using Taylor's method to solve the initial value problem

y'(t)=f(t,y)=y-(t^2)+1, 0≤t≤4, y(0)=0.5

Solution

First, let's solve this using Mathematica:

Clear[f,t,y]
f[t_,y_]=y-t^2+1
DSolve[{y'[t]==f[t,y[t]],y[0]==1/2},y[t],t]
Simplify[%]

sol[t_]=-((E^t)/2)+(((1+t))^2)
plotsol=Plot[sol[t],{t,0,4},BaseStyle->Large,AxesLabel->{"t","y(t)"}]

Runge Kutta Method (Order Four)

Now, let's solve this using Runge Kutta Order Four.

f[t_,y_] = y-(t^2)+1;

h = 0.2;
num=20;

t=0.0;
w=0.5;

data={{t,w}};
Do[
k1=h*f[t,w];
k2=h*f[t+h/2,w+k1/2];
k3=h*f[t+h/2,w+k2/2];
k4=h*f[t+h,w+k3];
w=w+(k1+2k2+2k3+k4)/6.0;
t=t+h;
data=Append[data,{t,w}],
{i,1,num}]

TableForm[data,TableHeadings->{None,{"(Subscript[t, i])","(Subscript[y, i])"}}]

rk4plot = ListPlot[data,Frame->True,FrameLabel->{"t","y"},PlotStyle->{RGBColor[0,0,1],PointSize[0.02]},BaseStyle->20];
Show[rk4plot,plotsol]

error = Table[{data[[i,1]],Abs[data[[i,2]]-sol[data[[i,1]]]]},{i,1,num}];
TableForm[error,TableHeadings->{None,{"(Subscript[t, i])","|(Subscript[y, i])-y((Subscript[t, i]))|"}}]

rk4errorplot = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{RGBColor[0,0,1],PointSize[0.02]},BaseStyle->20]

Taylor's Method n=4

Clear[f,h]
TaylorMethod[t0_,y0_,n_]:=Module[{t,y,data,Derries},
(Subscript[Derries, 0])[t_,y_]=f[t,y];
t[0]=t0;
y[0]=y0;
Do[{(Subscript[Derries, k])[t_,y_]=D[(Subscript[Derries, k-1])[t,y[t]],t]/.y'[t]->f[t,y]},{k,1,n-1}];

t[i_]:=t[i]=t[i-1]+h;
y[i_]:=y[i]=y[i-1]+(∑_(k=1)^n)((h^k)/(k!))*(Subscript[Derries, k-1])[t[i-1],y[i-1]];
Return[{t[1],y[1]}];
]

Clear[t,y,f,h,data]

f[t_,y_] = y-(t^2)+1;
t[0]=0.0;
y[0]=0.5;

h = 0.2;
num=20;
n=1;
k=4;

data={{t[0],y[0]}};
Do[
t[i]=TaylorMethod[t[i-1],y[i-1],k][[1]];
y[i]=TaylorMethod[t[i-1],y[i-1],k][[2]];
data=Append[data,{t[i],y[i]}];
,{i,1,num}]

TableForm[data,TableHeadings->{None,{"(Subscript[t, i])","(Subscript[y, i])"}}]

Here are some plots:

Taylorplot4 = ListPlot[data,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Red,PointSize[0.02]},BaseStyle->20];
Show[Taylorplot4,plotsol]

error = Table[{t[i],Abs[y[i]-sol[t[i]]]},{i,0,num}];
TableForm[error,TableHeadings->{None,{"(Subscript[t, i])","|(Subscript[y, i])-y((Subscript[t, i]))|"}}]

errorplot4 = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{Red,PointSize[0.02]},BaseStyle->20]

Comparison

Show[rk4errorplot,errorplot4]