Numerical Solutions to IVPs

Notebook Author/Course Instructor: Barry McQuarrie

Date: Sep 2015

Taylor's Method of Order n

We want to find a numerical approximation to the solution to the initial value problem:

y'(t)=f(t,y), a≤t, y(a)=c

Derivation of the Numerical Approximation (Done in Class)

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),i=0,1,2,3,...
(Subscript[y, i])=(Subscript[y, i-1])+(∑_(k=1)^n)(((f^((k-1)))((Subscript[t, i-1]),(Subscript[y, i-1])))/(k!))(h^k), i=1,2,3,4,...

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]),(Subscript[y, i]))

For n=2 you get:

(Subscript[y, i+1])=(Subscript[y, i])+h f((Subscript[t, i]),(Subscript[y, i]))+((h^2)/2)(d/(d t))f((Subscript[t, i]),(Subscript[y, i]))
=(Subscript[y, i])+h f((Subscript[t, i]),(Subscript[y, i]))+((h^2)/2)((∂/(∂ t))f((Subscript[t, i]),(Subscript[y, i]))+((∂/(∂ y))f((Subscript[t, i]),(Subscript[y, i])))(dy/dt))
=(Subscript[y, i])+h f((Subscript[t, i]),(Subscript[y, i]))+((h^2)/2)((∂/(∂ t))f((Subscript[t, i]),(Subscript[y, i]))+((∂/(∂ y))f((Subscript[t, i]),(Subscript[y, i])))f((Subscript[t, i]),(Subscript[y, 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!).

Example

We want to find a numerical approximation to the solution to the initial value problem:

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

Let's see if Mathematica can solve this:

Clear[t,y]
DSolve[{y'[t]==f[t,y[t]],y[0]==1/2},y[t],t]

Nope. So we definitely needed a numerical solution here.

Here is one way to implement this in Mathematica, using recursion relations:

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

h=0.1;
t[0]=-0.5;
y[0]=1.0;
t[i_]:=t[i]=t[i-1]+h
y[i_]:=y[i]=y[i-1]+h*f[t[i-1],y[i-1]]

euler=Table[{t[k],y[k]},{k,0,50}]
eulerplot=ListPlot[euler,PlotStyle->{Red},AxesLabel->{"t","y"}]

To get the Taylor Method of Order 2,  we need to do partial derivatives in Mathematica--here is a good way to do them:

Derivative[1,0][f][t,y]
Derivative[0,1][f][t,y]

Now, we can implement the method:

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

h=0.1;
t[0]=-0.5;
y[0]=1.0;
t[i_]:=t[i]=t[i-1]+h
y[i_]:=y[i]=y[i-1]+h*f[t[i-1],y[i-1]]+((h^2)/2)(Derivative[1,0][f][t[i-1],y[i-1]]+Derivative[0,1][f][t[i-1],y[i-1]]*f[t[i-1],y[i-1]])

taylor2=Table[{t[k],y[k]},{k,0,50}]
taylor2plot=ListPlot[taylor2,PlotStyle->{Blue},AxesLabel->{"t","y"}]

Show[eulerplot,taylor2plot,ImageSize->{600,Automatic},BaseStyle->20,PlotRange->All]

We can take a smaller value of h to improve the results in both cases:

Clear[f,t,y]
f[t_,y_]=y-(t^2)+Sin[t y] t;
h=0.05;
t[0]=-0.5;
y[0]=1.0;
t[i_]:=t[i]=t[i-1]+h
y[i_]:=y[i]=y[i-1]+h f[t[i-1],y[i-1]]
euler=Table[{t[k],y[k]},{k,0,100}];
eulerplot=ListPlot[euler,PlotStyle->{Red},AxesLabel->{"t","y"}];


Clear[f,t,y]
f[t_,y_]=y-(t^2)+Sin[t y] t;
h=0.05;
t[0]=-0.5;
y[0]=1.0;
t[i_]:=t[i]=t[i-1]+h
y[i_]:=y[i]=y[i-1]+h f[t[i-1],y[i-1]]+(1/2) (h^2) (Derivative[1,0][f][t[i-1],y[i-1]]+Derivative[0,1][f][t[i-1],y[i-1]] f[t[i-1],y[i-1]])
taylor2=Table[{t[k],y[k]},{k,0,100}];
taylor2plot=ListPlot[taylor2,PlotStyle->{Blue},AxesLabel->{"t","y"}];

Show[eulerplot,taylor2plot,BaseStyle->20]

It is a good idea to vary h to ensure you have chosen small enough to get accurate results

You might also wish to overlay the IVP solution on a direction field.

fieldeqn=VectorPlot[{1,f[t,y]},{t,-2,6},{y,0,10},VectorScale->{Tiny,Small,None}];
Show[fieldeqn,eulerplot,taylor2plot,BaseStyle->20]

Investigating Error Estimates

Here is an initial value problem that Mathematica can solve:

Clear[t,y]
f[t_,y_]=Sin[t]*Exp[-t]-Cos[4t];
DSolve[{y'[t]==f[t,y[t]],y[0]==1/2},y[t],t]

ysol[t_]=-(1/4) (E^(-t)) (-4 (E^t)+2 Cos[t]+2 Sin[t]+(E^t) Sin[4 t])

Let's see how the error in our numerical approximation look when compared to the exact solution. We want the error as a function of h, so let's look at y(2) for different values of h.

sol=ysol[2.0]

You can ignore all the Mathematica programming if you like--the important thing is the table at the end. This just changes h values and computes y(2) using both Euler’s method and Taylor’s method of Order 2.

Clear[h]
error={};
For[
h=1,h>(10^(-5)),h=h*(10^(-1)),

Clear[t,y,y2];

t[0]=0;
t[i_]:=t[i]=t[i-1]+h;

y[0]=0.5;
y[i_]:=y[i]=y[i-1]+h f[t[i-1],y[i-1]];

y2[0]=0.5;y2[i_]:=y2[i]=y2[i-1]+h*f[t[i-1],y2[i-1]]+((h^2)/2)(Derivative[1,0][f][t[i-1],y2[i-1]]+Derivative[0,1][f][t[i-1],y2[i-1]]*f[t[i-1],y2[i-1]]);


Table[{t[j],y[j],y2[j]},{j,0,(2/h)}];

error=Append[error,{h,Abs[((y[(2/h)]-sol)/sol)],Abs[((y2[(2/h)]-sol)/sol)]}]
]

Style[TableForm[N[error],TableHeadings->{None,{"h","n=1","n=2"}}],Large]

The n=1 gains roughy a decimal of accuracy with every decimal decrease in h.

The n=2 gains two decimals of accuracy!

This is why we might want to use n=3 or even n=4. Round off error will put a limit on your gains, however.

For Your Information: Implementing Taylor's Method of Order 4 Using Mathematica

In practice, Taylor's Method of Order 4 is commonly used, or other advanced numerical techniques like Runga-Kutta or Adams-Bashforth (a predictor-corrector technique). We could certainly work out Taylor's Method of Order 4 by hand, but here is how Mathematica can help you set this up. I use one of Mathematica's advanced programming techniques, the Module, to do this.

Determining the derivative in 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, which is done recursively.

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

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

Taylor's Method n=4

The following bit of code requires the earlier Module to be defined.

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

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

h = 0.05;
num=100;
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}]

There is no output from the above, but the Taylor's Method of Order 4 approximation to the solution is now stored in data. Notice this is much slower than Taylor’s Method of Order 2!

Here is a table of the numerical approximation to the solution:

taylor4plot=ListPlot[data,PlotStyle->{Orange},AxesLabel->{"t","y"}];
Show[eulerplot,taylor2plot,taylor4plot,BaseStyle->20]