Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2016
Lecture topic: Predictor Corrector (multi-step methods)

Clear[tp]
sol=Integrate[(t-tp[0])(t-tp[1])(t-tp[2]),{t,tp[1],tp[2]}]
Do[sol=sol/.tp[i]->h+tp[i-1],{i,2,1,-1}]
Simplify[sol]

Deriving Some Formulas

Adams-Bashforth Two Step Explicit Method

Here, the equation is derived using the Lagrange formalism.

Clear[y,f,approx,t,tp]
Lag[n_,k_,t_]:=((∏_(i=0)^(k-1))((t-tp[i])/(tp[k]-tp[i])) )( (∏_(i=k+1)^n)((t-tp[i])/(tp[k]-tp[i])))

Here is the formula that we use to approximate the y((Subscript[t, 2])):

approx[n_]:=y[tp[1]]+(∑_(k=0)^n)f[tp[k],y[tp[k]]]*(∫_(tp[1])^(tp[2]))Lag[n,k,t]\[DifferentialD]t

approx[1]

approx[1]=approx[1]/.tp[2]->h+tp[1]/.tp[1]->h+tp[0]

Simplify[approx[1]]

This is the Adams-Bashforth Two Step Explicit Method (Equation 6.72).

Adams-Bashforth Explicit Methods

Here is a little module that returns the Adams Bashforth explicit m-step methods, and the local truncation error:

ABExplicitmethod[n_]:=Module[{approx,locerror,m=n-1},

approx=y[tp[n-1]]+(∑_(k=0)^m)f[tp[k],y[tp[k]]]*(∫_(tp[n-1])^(tp[n]))((∏_(i=0)^(k-1))((t-tp[i])/(tp[k]-tp[i])) )( (∏_(i=k+1)^m)((t-tp[i])/(tp[k]-tp[i])))\[DifferentialD]t;

Do[approx=approx/.tp[i]->h+tp[i-1]
,{i,n+1,1,-1}];
 
locerror=((Derivative[n+1][y][μ])/((n)!))(∫_(tp[n-1])^(tp[n]))((∏_(i=0)^m)(t-tp[i]))\[DifferentialD]t;

Do[locerror=locerror/.tp[i]->h+tp[i-1]
,{i,n+1,1,-1}];

Print["Local Truncation Error: ",Simplify[locerror] ];
Return[Simplify[approx]];

]

Here is Equation 6.72, Adams Bashforth Two Step Explicit to compute (Subscript[y, 2])=y((Subscript[t, 0])+2h)

ABExplicitmethod[2]

Here is it rewritten (makes it a bit easier to compare to the text):

(Subscript[y, 2])=(Subscript[y, 1])+(h/2) ( 3 (Subscript[f, 1])-(Subscript[f, 0]))

Here is Equation 6.87, Adams Bashforth Three Step Explicit:

ABExplicitmethod[3]

Rewritten:

(Subscript[y, 3])=(Subscript[y, 2])+(h/12) (5 (Subscript[f, 0])-16  (Subscript[f, 1])+23  (Subscript[f, 2]))

Here is Equation 6.88, Adams Bashforth Four Step Explicit:

ABExplicitmethod[4]

Rewritten:

(Subscript[y, 4])=(Subscript[y, 3])+(h/24) (-9  (Subscript[f, 0])+37  (Subscript[f, 1])-59  (Subscript[f, 2])+55  (Subscript[f, 3]))

Here is Adams Bashforth Five Step Explicit:

ABExplicitmethod[5]

Rewritten:

(Subscript[y, 5])=(Subscript[y, 4])+(h/720) (251  (Subscript[f, 0])-1274  (Subscript[f, 1])+2616  (Subscript[f, 2])-2774  (Subscript[f, 3])+1901  (Subscript[f, 4]))

Adams-Moulton Implicit Methods

Here is a little module that returns the Adams Moulton implicit m-step methods, and the local truncation error. The change from explicit methods is in red.

AMImplicitmethod[n_]:=Module[{approx,locerror,m=n},

approx=y[tp[n-1]]+(∑_(k=0)^m)f[tp[k],y[tp[k]]]*(∫_(tp[n-1])^(tp[n]))((∏_(i=0)^(k-1))((t-tp[i])/(tp[k]-tp[i])) )( (∏_(i=k+1)^m)((t-tp[i])/(tp[k]-tp[i])))\[DifferentialD]t;

Do[approx=approx/.tp[i]->h+tp[i-1]
,{i,n+1,1,-1}];
 
locerror=((Derivative[n+1][y][μ])/((n)!))(∫_(tp[n-1])^(tp[n]))((∏_(i=0)^m)(t-tp[i]))\[DifferentialD]t;

Do[locerror=locerror/.tp[i]->h+tp[i-1]
,{i,n+1,1,-1}];

Print["Local Truncation Error: ",Simplify[locerror] ];
Return[Simplify[approx]];

]

Here is Equation 6.90, Adams Moulton Two Step Implicit:

AMImplicitmethod[2]//Simplify

Rewritten

(Subscript[y, 2])=(Subscript[y, 1])+(h/12) (-(Subscript[f, 0])+8(Subscript[f, 1])+5(Subscript[f, 2]))

Here is Equation 6.95, Adams Moulton Three Step Implicit:

AMImplicitmethod[3]

Rewritten

(Subscript[y, 3])=(Subscript[y, 2])+(h/24) ((Subscript[f, 0])-5 (Subscript[f, 1])+19 (Subscript[f, 2])+9 (Subscript[f, 3]))

Here is Equation 6.96, Adams Moulton Four Step Implicit:

AMImplicitmethod[4]

Rewritten

(Subscript[y, 4])=(Subscript[y, 3])+(h/720) (-19 (Subscript[f, 0])+106 (Subscript[f, 1])-264 (Subscript[f, 2])+646 (Subscript[f, 3])+251 (Subscript[f, 4]))

Here is Adams Moulton Five Step Implciit:

AMImplicitmethod[5]

Rewritten

(Subscript[y, 5])=(Subscript[y, 4])+(h/1440)(27 (Subscript[f, 0])-173 (Subscript[f, 1])+482 (Subscript[f, 2])-798 (Subscript[f, 3])+1427 (Subscript[f, 4])+475 (Subscript[f, 5]))

Adams Bashforth Moulton Four Step Predictor Corrector

The initial value problem we wish to solve is

y'(t)=f(t,y)=y-(t^3)cos(6t)sin(5t), 0≤t≤10, y(0)=0.5

Solution

First, let's solve this using Mathematica:

Clear[f,t,y]
f[t_,y_]=-y-t^3 Cos[6t]Sin[5t];
DSolve[{y'[t]==f[t,y[t]],y[0]==1/2},y[t],t];
Simplify[%]

Clear[sol]
sol[t_]=(1/55383364)((E^(-t)) (27687722-13845841 (E^t) t (3-3 t+(t^2)) Cos[t]+11 (E^t) (360-10797 t-11163 (t^2)+226981 (t^3)) Cos[11 t]+41537523 (E^t) Sin[t]-41537523 (E^t) t Sin[t]+13845841 (E^t) (t^3) Sin[t]+10437 (E^t) Sin[11 t]+33123 (E^t) t Sin[11 t]-669780 (E^t) (t^2) Sin[11 t]-226981 (E^t) (t^3) Sin[11 t]));
plotsol=Plot[sol[t],{t,0,10},PlotStyle->Black]

Runge Kutta Method (Order Four)

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

Clear[f,t,w,h]
f[t_,y_] = -y-t^3 Cos[6t]Sin[5t];

h = 0.1;
num=100;

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[Take[data,{1,num,10}],TableHeadings->{None,{"(Subscript[t, i])","(Subscript[y, i])"}}]

rk4plot = ListPlot[data,Frame->True,FrameLabel->{"t","error"},PlotStyle->{Blue,PointSize[0.02]}];
Show[plotsol,rk4plot,ImageSize->{500,Automatic},BaseStyle->20]

error = Table[{data[[i,1]],Abs[data[[i,2]]-sol[data[[i,1]]]]},{i,1,num}];
rk4errorplot = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{Blue,PointSize[0.02]}];
Show[rk4errorplot,ImageSize->{500,Automatic},BaseStyle->20]

Adams Bashforth Four Step Explicit Method

Here is an implementation of the Adams Bashforth four step explicit method. Notice we need a Runge Kutta method to get the ball rolling (find the first few values). The code in orange is the initialization by using the one step Runge Kutta method, the code in  blue is the Adams Bashforth explicit method.

Clear[f,data,t,w,h]
f[t_,y_] = -y-t^3 Cos[6t]Sin[5t];
h = 0.1;
num=100;

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

t0=data[[1,1]];
t1=data[[2,1]];
t2=data[[3,1]];
t3=data[[4,1]];
w0=data[[1,2]];
w1=data[[2,2]];
w2=data[[3,2]];
w3=data[[4,2]];

Do[
t=t+h;
w=w3+h(55f[t3,w3]-59f[t2,w2]+37f[t1,w1]-9f[t0,w0])/24.0;
t0=t1;
t1=t2;
t2=t3;
t3=t;
w0=w1;
w1=w2;
w2=w3;
w3=w;
data=Append[data,{t,w}],
{i,4,num}]


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

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

error = Table[{data[[i,1]],Abs[data[[i,2]]-sol[data[[i,1]]]]},{i,1,num}];

AB4errorplot = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{Red,PointSize[0.02]}];
Show[AB4errorplot,ImageSize->{500,Automatic},BaseStyle->20]

Adams Predictor Corrector Method

Here is an implementation of the Adams fourth order predictor corrector method. Notice we need a Runge Kutta method to get the ball rolling (find the first few values). The code in orange is the initialization by using the one step Runge Kutta method, the code in  blue is the predictor step using the four step Adams Bashforth explicit method, and the code in red is the corrector step using the three step Adams Moulton implicit method. The only difference from the explicit method is the addition of the red code as a corrector.

Clear[f,t,w,h,data]
f[t_,y_] = -y-t^3 Cos[6t]Sin[5t];

h = 0.1;
num=100;

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

t0=data[[1,1]];
t1=data[[2,1]];
t2=data[[3,1]];
t3=data[[4,1]];
w0=data[[1,2]];
w1=data[[2,2]];
w2=data[[3,2]];
w3=data[[4,2]];

Do[
t=t+h;
w=w3+h(55f[t3,w3]-59f[t2,w2]+37f[t1,w1]-9f[t0,w0])/24.0;
w=w3+h(9f[t,w]+19f[t3,w3]-5f[t2,w2]+f[t1,w1])/24.0;
t0=t1;
t1=t2;
t2=t3;
t3=t;
w0=w1;
w1=w2;
w2=w3;
w3=w;
data=Append[data,{t,w}],
{i,4,num}]

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

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

error = Table[{data[[i,1]],Abs[data[[i,2]]-sol[data[[i,1]]]]},{i,1,num}];

AB4PCerrorplot = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{RGBColor[1,0,1],PointSize[0.02]}];
Show[AB4PCerrorplot,ImageSize->{500,Automatic},BaseStyle->20]

Comparison

Show[rk4errorplot,AB4errorplot,AB4PCerrorplot,PlotRange->{Automatic,{0,12.0}},ImageSize->{500,Automatic},BaseStyle->20]

Notice that the predictor-corrector method is a definite improvement over the Adams Bashforth method, but the Runge Kutta is still better in this example.

Adams Bashforth Moulton Five Step Predictor Corrector

The initial value problem we wish to solve is

y'(t)=f(t,y)=y-(t^3)cos(6t)sin(5t), 0≤t≤10, y(0)=0.5

Adams Predictor Corrector Method

Here is an implementation of the Adams fifth order predictor corrector method. The code in orange is the initialization by using the one step Runge Kutta method, the code in  blue is the predictor step using the five step Adams Bashforth explicit method, and the code in red is the corrector step using the four step Adams Moulton implicit method iterated twice. Iterating twice gives a slight improvement, but realize that this is refining towards the AB five step explicit approximation, not the solution itself. It's better in practice to just decrease the step-size over doing multiple corrections.

Clear[f,t,w,h,data]
f[t_,y_] = -y-t^3 Cos[6t]Sin[5t];

h = 0.1;
num=100;

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

t0=data[[1,1]];
t1=data[[2,1]];
t2=data[[3,1]];
t3=data[[4,1]];
t4=data[[5,1]];
w0=data[[1,2]];
w1=data[[2,2]];
w2=data[[3,2]];
w3=data[[4,2]];
w4=data[[5,2]];

Do[
t=t+h;
w=w4+(1/720) (251 h f[t0,w0]-1274 h f[t1,w1]+2616 h f[t2,w2]-2774 h f[t3,w3]+1901 h f[t4,w4]);

w=w4+(1/720) (-19 h f[t1,w1]+106 h f[t2,w2]-264 h f[t3,w3]+646 h f[t4,w4]+251 h f[t,w]);
w=w4+(1/720) (-19 h f[t1,w1]+106 h f[t2,w2]-264 h f[t3,w3]+646 h f[t4,w4]+251 h f[t,w]);

t0=t1;
t1=t2;
t2=t3;
t3=t4;
t4=t;
w0=w1;
w1=w2;
w2=w3;
w3=w4;
w4=w;
data=Append[data,{t,w}],
{i,5,num}]

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

AB5PCplot = ListPlot[data,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Darker[Green],PointSize[0.02]}];
Show[AB5PCplot,plotsol,ImageSize->{500,Automatic},BaseStyle->20]

error = Table[{data[[i,1]],Abs[data[[i,2]]-sol[data[[i,1]]]]},{i,1,num}];

AB5PCerrorplot = ListPlot[error,Frame->True,FrameLabel->{"t","error"},PlotStyle->{Darker[Green],PointSize[0.02]}];
Show[AB5PCerrorplot,ImageSize->{500,Automatic},BaseStyle->20]

Comparison

Show[rk4errorplot,AB4PCerrorplot,AB5PCerrorplot,PlotRange->{Automatic,{0,2.0}},ImageSize->{500,Automatic},BaseStyle->20]

There is an improvement in using fifth order predictor-corrector, but Runge Kutta still kicks ass.