Example of Taylor Methods of Higher Order

Andy Long
MAT360
Spring, 2017

An Example:

Here is the true solution to the differential equation 
y’=1-y
with initial condition y(0)=0:

Clear[h,y,w,ww,www]
y[i_]:=1-Exp[-i*h]
Table[Series[y[i],{h,0,5}],{i,0,10,2}]//TableForm

Here is the Euler approximate solution: h+(1-h)w[i-1]

w[i_]:=If[i>0,(*then*)
w[i-1]+h(1-w[i-1]),(*else*) 
0]
Table[Expand[w[i]],{i,0,10,2}]//TableForm

Here is the Taylor Order-2 approximation: h(1-h/2)+(1-h+h^2/2)ww[i-1]

Clear[h]
Series[h(1-h/2)+(1-h+h^2/2)a,{h,0,3}]
ww[i_]:=If[i>0,a=ww[i-1]+h(1-ww[i-1]+(h/2) (-1+ww[i-1]) ),0]
Table[Expand[ww[i]],{i,0,10,2}]//TableForm

Here is the Taylor Order-3 approximation: h(1-h/2+h^2/6)+(1-h+h^2/2+h^3/6)www[i-1]

Clear[a,h]
www[i_]:=If[i>0,a=www[i-1];a+h(1-a+h((1/2) (-1+a) +(h/6) (1+a))) ,0]
Table[Expand[www[i]],{i,0,10,2}]//TableForm

Here we plot the approximations, against the true solution;

n=10;
h=1/n;
ListPlot[{Table[y[i],{i,0,n}],Table[w[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}],Table[ww[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}],Table[www[i],{i,0,n}]},Joined->True]

And now we plot the errors:

n=10;
h=1/n;
ListPlot[{Table[y[i],{i,0,n}]-Table[w[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}]-Table[ww[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}]-Table[www[i],{i,0,n}]},Joined->True]

Our e^x Example: we want to reproduce Exp[x] exactly on the interval [0,1].

Here is the true solution to the differential equation 
y’=y
with initial condition y(0)=1:

Clear[h,y,w,ww,www]
y[i_]:=Exp[i*h]
init=y[0]
Table[Series[y[i],{h,0,5}],{i,0,10,2}]//TableForm

Here is the Euler approximate solution: h+(1-h)w[i-1]

w[i_]:=If[i>0,(*then*)
w[i-1](1+h),(*else*) 
init]
Table[Expand[w[i]],{i,0,10,2}]//TableForm

Here is the Taylor Order-2 approximation:h(1-h/2)+(1-h+h^2/2)ww[i-1]

Clear[h]
Series[h(1-h/2)+(1-h+h^2/2)a,{h,0,3}]
ww[i_]:=If[i>0,a=ww[i-1](1+h(1+(h/2) ) ),init]
Table[Expand[ww[i]],{i,0,10,2}]//TableForm

Here is the Taylor Order-3 approximation: h(1-h/2+h^2/6)+(1-h+h^2/2+h^3/6)www[i-1]

Clear[a,h]
www[i_]:=If[i>0,a=www[i-1](1+h(1+h((1/2) (1 +(h/3) )))) ,init]
Table[Expand[www[i]],{i,0,10,2}]//TableForm

Here is the Taylor Order-4 approximation: (1+ h(1+h/2(1+h/3(1+h/4)) wwww[i-1]

Clear[a,h]
wwww[i_]:=If[i>0,a=wwww[i-1](1+h(1+(h/2)(1 +(h/3)(1+h/4) ))),init]
Table[Expand[wwww[i]],{i,0,10,2}]//TableForm

Here we plot the approximations, against the true solution;

n=10;
h=1/n;
ListPlot[{Table[y[i],{i,0,n}],Table[w[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}],Table[ww[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}],Table[www[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}],Table[wwww[i],{i,0,n}]},Joined->True]

And now we plot the errors:

n=10;
h=1/n;
ListPlot[{Table[y[i],{i,0,n}]-Table[w[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}]-Table[ww[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}]-Table[www[i],{i,0,n}]},Joined->True]
ListPlot[{Table[y[i],{i,0,n}]-Table[wwww[i],{i,0,n}]},Joined->True]