The Taylor Series Method For ODE

This notebook gives some examples for the Taylor method for solving ordinary differential equations. We start with Euler's method.

The Euler Method for Solving ODE numerically

The idea of Euler is that 

x'(t) = f(t,x)

so 

x(t+h) = x(t) + h x'(t) + O(h^2)

           = x(t) + h f(t,x) + O(h^2)
           
We can use this to define the Euler step from x(t) to x(t+h) by 

x -> x + h f(t,x)

Euler[f_,x0_,t0_,t1_,n_] := Module[{EulerStep,h},

h = (t1 - t0)/(n-1);
EulerStep[{t_,x_,h_}] := {t+h,x +h*f[t,x],h};
NestList[EulerStep,{ t0,x0,h},n-1][[All,1;;2]]

];

?NestList

Now we try applying this to the ODE:

	x'(t) = x(t),  t0 = 0, x(0) = 1. 
	
We already know the solution, but we try it with Mathematica just to see how it works:

ExactSoln = (x /. DSolve[{x'[t] == x[t],x[0] == 1},x,t])[[1]]

We now need to use our Euler code. For this equation, we write

f[t_,x_] := x;
Euler[f,1.0,0.0,1.0,10]

This is easier to understand in a plot:

Manipulate[Show[Plot[ExactSoln[t],{t,0,1}],ListLinePlot[Euler[f,1.0,0.0,1.0,n],PlotMarkers->Automatic],

AxesLabel->{Style["t",Larger,FontFamily->"Helvetica",Bold],Style["x",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Exact and Euler Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")],

{{n,10},2,100,1}]

We can go ahead and look at just the result of the calculation as well:

EulerVal[f_,x0_,t0_,t1_,n_] := Module[{EulerStep,h},

h = (t1 - t0)/(n-1);
EulerStep[{t_,x_,h_}] := {t+h,x +h*f[t,x],h};
Nest[EulerStep,{ t0,x0,h},n-1][[2]]

];

EulerVal[f,1.0,0.0,1.0,100]

ExactSoln[1.0]

As you can see, even for this comparatively trivial case, there aren't a lot of correct digits in the solution here. We go ahead and make the plot of correct digits just to see how fast this method is converging to a correct solution.

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];
ListLinePlot[Table[{n,Digits[EulerVal[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct Digits in Euler Solution of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

The Taylor Method in General

The idea of the higher order Taylor methods is that 

x'(t) = f(t,x)

so 

x(t+h) = x(t) + h x'(t) + (h^2/2) x''(t) + O(h^3)

           = x(t) + h f(t,x) + ?
           
We know if 

x'(t) = f(t,x)

we can differentiate both sides with respect to t to get

x''(t) = f_t(t,x) + f_x(t,x) * x'(t) = f_t(t,x) + f_x(t,x) * f(t,x).
           
We can use this to define the Taylor step from x(t) to x(t+h) by 

x -> x + h f(t,x) + (h^2/2) * (f_t(t,x) + f_x(t,x) * f(t,x)).

Taylor2[f_,x0_,t0_,t1_,n_] := Module[{Taylor2Step,h,ft,fx,t,x,u,v},

ft[t_,x_]:= D[f[u,v],u] /. {u->t,v->x};
fx[t_,x_]:= D[f[u,v],v] /. {u->t,v->x};

h = (t1 - t0)/(n-1);
Taylor2Step[{t_,x_,h_}] := {t+h,x +h*f[t,x] + (h^2/2)*(ft[t,x] + fx[t,x]*f[t,x]),h};
NestList[Taylor2Step,{ t0,x0,h},n-1][[All,1;;2]]

];

Taylor2[f,1.0,0.0,1.0,10]

As before, we'll show this in a plot

Manipulate[Show[Plot[ExactSoln[t],{t,0,1}],ListLinePlot[Euler[f,1.0,0.0,1.0,n],PlotMarkers->Automatic],
ListLinePlot[Taylor2[f,1.0,0.0,1.0,n],PlotMarkers->Automatic,PlotStyle->{Red}],
AxesLabel->{Style["t",Larger,FontFamily->"Helvetica",Bold],Style["x",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Exact, Euler, and Taylor2 Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")],
{{n,10},2,100,1}]

We can also check the values for the solver:

Taylor2Val[f_,x0_,t0_,t1_,n_] := Module[{Taylor2Step,h,ft,fx,t,x,u,v},

ft[t_,x_]:= D[f[u,v],u] /. {u->t,v->x};
fx[t_,x_]:= D[f[u,v],v] /. {u->t,v->x};

h = (t1 - t0)/(n-1);
Taylor2Step[{t_,x_,h_}] := {t+h,x +h*f[t,x] + (h^2/2)*(ft[t,x] + fx[t,x]*f[t,x]),h};
Nest[Taylor2Step,{ t0,x0,h},n-1][[2]]

];

Taylor2Val[f,1.0,0.0,1.0,10]

ExactSoln[1.0]

We'll now make a plot of the correct digits in the Euler and Taylor2 methods:

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];
Show[ListLinePlot[{Table[{n,Digits[EulerVal[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}],
Table[{n,Digits[Taylor2Val[f,1.0,0.0,1.0,n],ExactSoln[1.0]]},{n,10,1000,10}]},
Filling->Axis,PlotMarkers->Automatic,
AxesLabel->{Style["Steps (n)",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct Digits in Euler and Taylor2 Solutions of x' = x",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]]

Now we can try some fun stuff! Keep in mind that the Taylor method was comparatively easy to code because Mathematica is happy to evaluate the partial derivatives of f for us. In general (for instance, in something like a C++ environment, or where the function f(t,x) has no easy symbolic expression), this dependence on being able to symbolically compute partials will become a real problem.