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,n]
n=1000
b=1.0
a=0.0
h=(b-a)/n*1.0
eta=0
y[i_]:=1-Exp[-i*h]

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

w[0]=eta
i=1
h=(b-a)/n*1.0;
While[i<=n,w[i]=w[i-1]+h(1-w[i-1]);i++]

Here is the Taylor Order-2 approximation:

ww[0]=eta
i=1
While[i<=n,ww[i]=ww[i-1]+h(1-ww[i-1]+(h/2) (-1+ww[i-1]) );i++]

Here is the Taylor Order-3 approximation:

www[0]=eta
i=1
While[i<=n,www[i]=www[i-1]+h(1-www[i-1]+(h/2) (-1+www[i-1]+(h/3) (1-www[i-1])));i++]

Here we plot the approximations, against the true solution;

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:

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]
eta=y[0]
b=1.0
a=0.0
n=1000
h=(b-a)/n*1.0;

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

w[0]=eta
i=1
While[i<=n,w[i]=w[i-1](1+h);i++]

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

ww[0]=eta
i=1
While[i<=n,ww[i]=ww[i-1](1+h(1+(h/2) ));i++]

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

www[0]=eta
i=1
While[i<=n,www[i]=www[i-1](1+h(1+(h/2)(1+(h/3) )));i++]

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

wwww[0]=eta
i=1
While[i<=n,wwww[i]=wwww[i-1](1+h(1+(h/2)(1+(h/3)(1+(h/4) ) )));i++]

Here we plot the approximations, against the true solution;

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:

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]

Our e^x Example: small step-size on the interval [0,1], with error bound.

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]
eta=y[0]
b=1.0
a=0.0
n=10000
h=(b-a)/n*1.0;
bound=h/2*Exp[1](Exp[1]-1)

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

w[0]=eta
i=1
While[i<=n,w[i]=w[i-1](1+h);i++]

Here we plot the approximations, against the true solution;

ListPlot[{Table[y[i],{i,0,n}],Table[w[i],{i,0,n}]},Joined->True]

And now we plot the errors:

ListPlot[{Table[y[i],{i,0,n}]-Table[w[i],{i,0,n}],Table[bound,{i,0,n}]},Joined->True]