Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2015
Lecture topic: Numerical Solutions and Autonomous Systems

Taylor's Method of Order n for Systems

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

(Subscript[x, j])'(t)=(Subscript[f, j])(t,(Subscript[x, 1]),(Subscript[x, 2]),...) , j=1,2,3,4,...,N

(Subscript[x, j])(0)is given for all j.

Derivation of the Taylor’s Method of Order 2

Our approximation will be a set of points, ((Subscript[t, i]),(Subscript[(Subscript[x, 1]), i]),(Subscript[x, 2i])), i=0,1,2,3,...,N which is initialized by ((Subscript[t, 0]),(Subscript[(Subscript[x, 1]), 0]),(Subscript[x, 20]))=(a,b,c). The derivation is a simple extension of the result for Taylor's method of order n for one unknown y(t).

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

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

Now, from the differential equation, we have:

((Subscript[x, j])^((1)))(t)=(Subscript[f, j])(t,(Subscript[x, j])(t))
 ((Subscript[x, j])^((2)))(t)=((Subscript[f, j])^((1)))(t,(Subscript[x, j])(t))
 ((Subscript[x, j])^((3)))(t)=((Subscript[f, j])^((2)))(t,(Subscript[x, j])(t))
⋮
 ((Subscript[x, j])^((k)))(t)=((Subscript[f, j])^((k-1)))(t,(Subscript[x, j])(t))

 ((Subscript[x, j])^((k)))((Subscript[t, i]))=((Subscript[f, j])^((k-1)))((Subscript[t, i]),(Subscript[x, j])((Subscript[t, i]))) (derivative ((Subscript[f, j])^((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[(Subscript[x, 1]), 0])=b
(Subscript[(Subscript[x, 2]), 0])=c

(Subscript[t, i+1])=(Subscript[t, i])+h
(Subscript[((Subscript[x, j])), i+1])=(Subscript[((Subscript[x, j])), i])+(∑_(k=1)^n)((((Subscript[f, j])^((k-1)))((Subscript[t, i]),(Subscript[x, j])((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. If we choose n=2, we get Taylor's method of order 2. I am going to surpress the functional dependence for some clarity here, but everything on the right hand sides should be evaluated at the previous point, i. For example. ((Subscript[f, j])^((k-1)))=((Subscript[f, j])^((k-1)))((Subscript[t, i]),(Subscript[((Subscript[x, 1])), i]),(Subscript[((Subscript[x, 2])), i])).

(Subscript[((Subscript[x, j])), i+1])=(Subscript[((Subscript[x, j])), i])+(∑_(k=1)^2)(((Subscript[f, j])^((k-1)))/(k!))(h^k)
=(Subscript[((Subscript[x, j])), i])+h (Subscript[f, j])+((h^2)/2)(((∂(Subscript[f, j]))/(∂ t))+((∂(Subscript[f, j]))/(∂ (Subscript[x, 1])))·((Subscript[dx, 1])/dt)+((∂(Subscript[f, j]))/(∂ (Subscript[x, 2])))·((Subscript[dx, 2])/dt))
=(Subscript[((Subscript[x, j])), i])+h (Subscript[f, j])+((h^2)/2)(((∂(Subscript[f, j]))/(∂ t))+((∂(Subscript[f, j]))/(∂ (Subscript[x, 1])))·(Subscript[f, 1])+((∂(Subscript[f, j]))/(∂ (Subscript[x, 2])))·(Subscript[f, 2]))

So Taylor's Method of order 2 for a 2×2 system is:

h  (user defined)
(Subscript[t, 0])=a
(Subscript[(Subscript[x, 1]), 0])=b
(Subscript[(Subscript[x, 2]), 0])=c

(Subscript[t, i])=(Subscript[t, i-1])+h
(Subscript[((Subscript[x, 1])), i])=(Subscript[((Subscript[x, 1])), i-1])+h (Subscript[f, 1])+((h^2)/2)(((∂(Subscript[f, 1]))/(∂ t))+((∂(Subscript[f, 1]))/(∂ (Subscript[x, 1])))·(Subscript[f, 1])+((∂(Subscript[f, 1]))/(∂ (Subscript[x, 2])))·(Subscript[f, 2]))
(Subscript[((Subscript[x, 2])), i])=(Subscript[((Subscript[x, 2])), i-1])+h (Subscript[f, 2])+((h^2)/2)(((∂(Subscript[f, 2]))/(∂ t))+((∂(Subscript[f, 2]))/(∂ (Subscript[x, 1])))·(Subscript[f, 1])+((∂(Subscript[f, 2]))/(∂ (Subscript[x, 2])))·(Subscript[f, 2]))

Example for a Linear Nonhomogeneous System

(Subscript[x, 1])'(t)=(Subscript[f, 1])(t,(Subscript[x, 1]),(Subscript[x, 2]))=((Subscript[x, 1])^2)-(Subscript[x, 2])-t
 (Subscript[x, 2])'(t)=(Subscript[f, 2])(t,(Subscript[x, 1]),(Subscript[x, 2]))=(Subscript[x, 1])-((Subscript[x, 2])^2)+Sin[t]

(Subscript[x, 1])(0)=-1
(Subscript[x, 2])(0)=1

Here is one way to implement Taylor's method for systems in Mathematica, using recursion relations:

Clear[f1,f2,t,x1,x2,h]
f1[t_,x1_,x2_]=(x1^2)+x2-t;
f2[t_,x1_,x2_]=x1-(x2^2)+Sin[t];
h=0.02;
t[0]=0;
x1[0]=-1.;
x2[0]=1.;
t[i_]:=t[i]=t[i-1]+h
x1[i_]:=x1[i]=x1[i-1]+h f1[t[i-1],x1[i-1],x2[i-1]]+(1/2) (h^2) (Derivative[1,0,0][f1][t[i-1],x1[i-1],x2[i-1]]+Derivative[0,1,0][f1][t[i-1],x1[i-1],x2[i-1]] f1[t[i-1],x1[i-1],x2[i-1]]+Derivative[0,0,1][f1][t[i-1],x1[i-1],x2[i-1]] f2[t[i-1],x1[i-1],x2[i-1]])
x2[i_]:=x2[i]=x2[i-1]+h f2[t[i-1],x1[i-1],x2[i-1]]+(1/2) (h^2) (Derivative[1,0,0][f2][t[i-1],x1[i-1],x2[i-1]]+Derivative[0,1,0][f2][t[i-1],x1[i-1],x2[i-1]] f1[t[i-1],x1[i-1],x2[i-1]]+Derivative[0,0,1][f2][t[i-1],x1[i-1],x2[i-1]] f2[t[i-1],x1[i-1],x2[i-1]])
taylor2=Table[{x1[k],x2[k]},{k,0,100}];
taylor2plot=ListPlot[taylor2,PlotStyle->{PointSize[0.02],Red},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{All,All},BaseStyle->20]

Using NDSolve for Systems

Clear[x,y]
eqns={x'[t]==((x[t])^2)+y[t]-t,y '[t]==x[t]-((y[t])^2)+Sin[t],x[0]==-1.0,y[0]==1.0}
sol=First[{x,y}/.NDSolve[eqns,{x,y},{t,0,2}]]

ndsolveplot=ParametricPlot[{sol[[1]][t],sol[[2]][t]},{t,0,2},PlotRange->All,AxesLabel->{"x","y"},PlotRange->{All,All},BaseStyle->20]

Show[taylor2plot,ndsolveplot,BaseStyle->20]

Autonomous Systems

Example 1

x'=y-1
 y'=(x^2)-1

family=((y^2)/2)-y==((x^3)/3)-x+c

familylist=family/.c->Range[-10,10,2];

lines=ContourPlot[Evaluate[familylist],{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->Purple,BaseStyle->20];

separatrix=ContourPlot[((y^2)/2)-y==((x^3)/3)-x+(1/6),{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->{Thick,Red},BaseStyle->20];

nullclines=ContourPlot[{y-1==0,(x^2)-1==0},{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->Blue,BaseStyle->20];

CP1=Graphics[{Red,Disk[{1,1},0.15]}];
CP2=Graphics[{Red,Disk[{-1,1},0.15]}];

field=VectorPlot[{y-1,(x^2)-1},{x,-5,5},{y,-5,5},VectorScale->{Small,Small,None}];

Show[lines,separatrix,nullclines,CP1,CP2,field]

The blue lines are the nullclines, which is where one of the derivatives is zero (notice the arrows on the direction field are horizontal or vertical on the nullclines).

Where the nullclines intersect (the red dots) we have equilibrium points (critical points).

The equilibrium at (1,1) is unstable, the one at (-1,1) is stable.

The red line forms a separatrix, which bounds a basin of attraction for the points that remain close to (-1,1).

Example 2 Van der Pol Oscillator

x'=y
 y'=-x+(1-(x^2))y

Verify that we couldn’t solve the autonomous DE when we eliminated t:

DSolve[y'[x]==((-x+(1-(x^2))y[x])/(y[x])),y[x],x]

Let’s use NDSolve to generate a solution.

Clear[x,y]
eqns={x'[t]==y[t],y '[t]==-x[t]+(1-((x[t])^2))y[t],x[0]==-0.2,y[0]==0.2};
sol=First[{x,y}/.NDSolve[eqns,{x,y},{t,0,50}]];

ndsolveplot=ParametricPlot[{sol[[1]][t],sol[[2]][t]},{t,0,50},PlotRange->All,AxesLabel->{"x","y"},PlotRange->{All,All},BaseStyle->20,PlotStyle->Purple];

nullclines=ContourPlot[{y==0,-x+(1-(x^2))y==0},{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->Blue,BaseStyle->20];

CP1=Graphics[{Red,Disk[{0,0},0.15]}];

field=VectorPlot[{y,-x+(1-(x^2))y},{x,-5,5},{y,-5,5},VectorScale->{Small,Small,None}];

Show[nullclines,ndsolveplot,CP1,field]

The blue lines are the nullclines, which is where one of the derivatives is zero (notice the arrows on the direction field are horizontal or vertical on the nullclines).

Where the nullclines intersect (the red dot) we have equilibrium points (critical points).

The critical point is unstable.

The solution moves away from (0,0), then appears to approach a limit cycle. If we start at a different iniital point, we still aproach the same set of points in the limit cycle.

plot={};
Do[
Clear[x,y,sol];
eqns={x'[t]==y[t],y '[t]==-x[t]+(1-((x[t])^2))y[t],x[0]==RandomReal[{-5,5}],y[0]==RandomReal[{-5,5}]};
sol=First[{x,y}/.NDSolve[eqns,{x,y},{t,0,50}]];

ndsolveplot=ParametricPlot[{sol[[1]][t],sol[[2]][t]},{t,0,50},PlotRange->All,AxesLabel->{"x","y"},PlotRange->{All,All},BaseStyle->20,PlotStyle->Purple];

AppendTo[plot,ndsolveplot];

,{i,1,50}]

nullclines=ContourPlot[{y==0,-x+(1-(x^2))y==0},{x,-5,5},{y,-5,5},PlotRange->{{-5,5},{-5,5}},FrameLabel->{"x","y"},ContourStyle->Blue,BaseStyle->20];

CP1=Graphics[{Red,Disk[{0,0},0.15]}];

field=VectorPlot[{y,-x+(1-(x^2))y},{x,-5,5},{y,-5,5},VectorScale->{Tiny,Small,None}];

Show[nullclines,plot,CP1,field]