Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Oct 2016
Lecture topic: Initial Value Problems

Taylor's Method

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

Mathematica 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},PlotStyle->Red,BaseStyle->20]

Derivation of Taylor's Method

Our approximation will be a set of points, ((Subscript[t, i]),(Subscript[y, i])), i=0,1,2,3,...,N which is initialized by ((Subscript[t, 0]),(Subscript[y, 0]))=(a,c).

Expand y((Subscript[t, i+1])) as a Taylor series about (Subscript[t, i]):

y((Subscript[t, i+1]))~(∑_(k=0)^n)(((y^((k)))((Subscript[t, i])))/(k!))((((Subscript[t, i+1])-(Subscript[t, i])))^k)
y((Subscript[t, i+1]))~(∑_(k=0)^n)(((y^((k)))((Subscript[t, i])))/(k!))(h^k), where h = (Subscript[t, i+1])-(Subscript[t, i])

Now, from the differential equation, we have:

(y^((1)))(t)=f(t,y(t))
 (y^((2)))(t)=(f^((1)))(t,y(t))
 (y^((3)))(t)=(f^((2)))(t,y(t))
⋮
 (y^((k)))(t)=(f^((k-1)))(t,y(t))

 (y^((k)))((Subscript[t, i]))=(f^((k-1)))((Subscript[t, i]),y((Subscript[t, i]))) (derivative (f^((k-1))) is with respect to t)

Using this result, we have an expression for the solution to the initial value problem in terms of the recursion relations:

h  (user defined)
(Subscript[t, 0])=a
(Subscript[y, 0])=c

(Subscript[t, i+1])=(Subscript[t, i])+h
(Subscript[y, i+1])=(Subscript[y, i])+(∑_(k=1)^n)(((f^((k-1)))((Subscript[t, i]),y((Subscript[t, i]))))/(k!))(h^k)

This is the Taylor Method of Order n. Choosing different values of n gives you different numerical approximation schemes. For example, n=1 is Euler's method which you probably saw in Calculus:

h  (user defined)
(Subscript[t, 0])=a
(Subscript[y, 0])=c

(Subscript[t, i+1])=(Subscript[t, i])+h
(Subscript[y, i+1])=(Subscript[y, i])+h f((Subscript[t, i]),y((Subscript[t, i])))

For n=2 you get:

(Subscript[y, i+1])=(Subscript[y, i])+h f((Subscript[t, i]),y((Subscript[t, i])))+((h^2)/2)(d/(d t))f((Subscript[t, i]),y((Subscript[t, i])))
=(Subscript[y, i])+h f((Subscript[t, i]),y((Subscript[t, i])))+((h^2)/2)((∂/(∂ t))f((Subscript[t, i]),y((Subscript[t, i])))+((∂/(∂ y))f((Subscript[t, i]),y((Subscript[t, i]))))(dy/dt))
=(Subscript[y, i])+h f((Subscript[t, i]),y((Subscript[t, i])))+((h^2)/2)((∂/(∂ t))f((Subscript[t, i]),y((Subscript[t, i])))+((∂/(∂ y))f((Subscript[t, i]),y((Subscript[t, i]))))f((Subscript[t, i]),y((Subscript[t, i]))))

As you take larger n, you get more accurate results for a given step size h, but this must be balanced with the increased computational time and round-off error. Typically, n=4 is as high as anyone goes (but you could go higher!).

Symbolic Taylor's Method of Order n=k

Here is the symbolic way of generating Taylor's Method of order n: The only tricky bit is getting the derivatives of f calculated correctly.

Clear[f,h]
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]}];
]

TaylorMethod[t0,y0,1]
TaylorMethod[t0,y0,2]
TaylorMethod[t0,y0,3]
TaylorMethod[t0,y0,4]

Here is how we can pull out the (Subscript[t, 1]) and (Subscript[y, 1]) parts:

k=1;
TaylorMethod[t0,y0,k][[1]]
TaylorMethod[t0,y0,k][[2]]

Taylor's Method n=1 (Euler's Method)

Now, let's solve this using Taylor's Method.

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;

data={{t[0],y[0]}};
Do[
t[i]=TaylorMethod[t[i-1],y[i-1],n][[1]];
y[i]=TaylorMethod[t[i-1],y[i-1],n][[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:

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

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]))|"}}]

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

Taylor's Method n=2

The only change from above is highlighted in red.

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=2;

data={{t[0],y[0]}};
Do[
t[i]=TaylorMethod[t[i-1],y[i-1],n][[1]];
y[i]=TaylorMethod[t[i-1],y[i-1],n][[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:

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

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]))|"}}]

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

Taylor's Method n=4

The only change from above is highlighted in red.

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=4;

data={{t[0],y[0]}};
Do[
t[i]=TaylorMethod[t[i-1],y[i-1],n][[1]];
y[i]=TaylorMethod[t[i-1],y[i-1],n][[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->{Darker[Green],PointSize[0.02]}];
Show[Taylorplot4,plotsol,BaseStyle->20]

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->{Darker[Green],PointSize[0.02]},BaseStyle->20]

Comparison

Show[Taylorplot1,Taylorplot2,Taylorplot4,plotsol,PlotRange->All]

Show[errorplot1,errorplot2,errorplot4]

Global Truncation Error Plots as a function of h

Let's plot the error as a function of h, to get a sense of how well our error analysis we did in class works. Notice we don’t have easy access to the bounds on the second derivative of the solution that enters into the analysis, so we will just look for a rough estimate of the order of the error by seeing how well the error fits linear or quadratic functions of h.

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

errorh={};
h=1/2;
While[h>(10^(-3)),

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

tfinal=1/2;
num=tfinal/h;

n=1;
k=1;

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

errorh=Append[errorh,{h,Abs[y[num]-sol[t[num]]]}];
h=h/2;
]

TableForm[errorh,TableHeadings->{None,{"h","Error in y(0.5)"}}]

errorhplot1 = ListPlot[errorh,Frame->True,FrameLabel->{"h","error in y(0.5)"},PlotStyle->{Red,PointSize[0.02]},BaseStyle->20]

The error is essentially linear in h (notice the fudge factor of 1/2.8 since we don’t know (y^((2)))((Subscript[c, i])) on the interval).

Plot[(h/2.7),{h,0,1/2},PlotStyle->Red];
Show[%,errorhplot1,BaseStyle->20]

errorhplot1log= ListLogLogPlot[errorh,Frame->True,FrameLabel->{"h","error in y(0.5)"},PlotStyle->{Red,PointSize[0.02]},BaseStyle->20]

Here is the result when k=2, so we expect the global trunctation error to be quadratic.

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

errorh={};
h=1/2;
While[h>(10^(-3)),

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

tfinal=1/2;
num=tfinal/h;

n=1;
k=2;

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

errorh=Append[errorh,{h,Abs[y[num]-sol[t[num]]]}];
h=h/2;
]

TableForm[errorh,TableHeadings->{None,{"h","Error in y(1/2)"}}]

errorhplot2 = ListPlot[errorh,Frame->True,FrameLabel->{"h","error in y(1/2)"},PlotStyle->{Blue,PointSize[0.02]},BaseStyle->20]

This error is essentially quadratic in h:

Plot[((h^2)/20),{h,0,1/2},PlotStyle->Blue];
Show[%,errorhplot2,BaseStyle->20]

errorhplot2log = ListLogLogPlot[errorh,Frame->True,FrameLabel->{"h","error in y(1/2)"},PlotStyle->{Blue,PointSize[0.02]},BaseStyle->20]

Show[errorhplot1log,errorhplot2log,PlotRange->All]