Numerical Approximation to First Order Systems of Equations

Heun’s Method

This method is similar to Euler’s Method in that you use successive approximations to generate a plot which is “close” to the solution of a differential equation.  In Heun’s Method, there are more intermediate steps before you get your next approximation.  These intermediate steps improve the overall performance of the method.  (See the lab for a greater explanation of the derivation of Heun’s Method.)

Heun[funct_,t0_,y0_,tfinal_,numsteps_]:=
Module[{f=funct, beginPoint=t0, w0=y0, endPoint=tfinal, n=numsteps,t,w,h,m1,m2,m3,j},
t=beginPoint;
w=w0;
h=(endPoint-beginPoint)/n;
tList={t};  (*Initialize our list of t-values *)
wList ={w}; (*Initialize our list of w-vaules: the approximations to the solution *)

(* Here we loop from 1 to n performing Heun's method and appending the results to the lists *)
For[j=1,j<n+1,j++,
m1=h*f[t,w];
m2=h*f[t+h/3,w+m1/3];
m3=h*f[t+2*h/3,w+2*m2/3];
t=t+h;
w=w+(m1+3*m3)/4;
tList=Append[tList,t];
wList=Append[wList,w];
];
Return[{tList,wList}]
]

Using Heun’s Method

Consider the IVP:  y'=y-(t^2)+1 where y(0)=(1/2) over the interval 0≤t≤2 with 10 subintervals.  Note that this is a non-autonomous differential equation.

mySol=y[t]/.DSolve[{y'[t]==y[t]-t^2+1,y[0]==1/2}, y[t],t][[1]]

Now, let’s use Heun’s Method for the approximation:

myF[t_,y_]=y-t^2+1;
t0=0.0;  (* This is the beginning point *)
y0=0.5;
endP=2.0;
mySteps=10.0;

approx1=Heun[myF, t0, y0, endP, mySteps]

We need to turn the two returned lists (of size mySteps+1) into one list of ordered pairs.

myPairs=Table[{approx1[[1,k]],approx1[[2,k]]}, {k, 1, mySteps+1}]

myH=(endP-t0)/mySteps
myplot10=ListPlot[myPairs,
Joined->True,
Mesh->All, PlotMarkers->{Automatic,Medium},
PlotStyle->{Blue,Thick},
AxesOrigin->{0,0},
PlotLegends->LineLegend[{Blue},{StringForm["h = ``", NumberForm[myH,3]]}]
];
(* Next Show our list of ordered pairs *)
Show[myplot10,
AxesLabel->{t,y},
AxesStyle->Bold,
PlotLabel->Style[Text["Heun's Method"],16]
]

dField=VectorPlot[
{1, myF[t,y]}, {t, -0.5, 2.5}, {y, -1, 6},
Axes->True,
AxesLabel->{t,y},
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->Directive[Darker[Green],Thick,"Line"],
StreamScale->None
];  (* Note that if I don't use a semi-colon here, it will save the plot and still show the plot. *)
Show[dField,myplot10,
PlotLabel->Style["Direction Field,\!\(\*FractionBox[\(\(\\ \)\(dy\)\), \(dt\)]\)=y-(t^2)+1",Bold,16]
]

Now For a System

Now lets try a systems example.  Suppose we wish to numerically approximate the IVP

({{({{y''-y y'}})+(y^2)=0,     y(0)=-1,    y'(0)=2.}})

To convert this to a system of first-order equations we let (Subscript[v, 1])=y and (Subscript[v, 2])=y' to get the system:

({{({{((d (Subscript[v, 1]))/dt)=(Subscript[v, 2])}, {((d (Subscript[v, 2]))/dt)=(Subscript[v, 1]) (Subscript[v, 2])-(Subscript[v, 1]^2)}})}})

We will need to use a Module for our RHS function so that we can use a vector (list) for the dependent variables:

myg[t_,v_]:=Module[{v1=v[[1]],v2=v[[2]]},
{v2,v1 v2-v1^2}
]

Notice that this is an autonomous system, so the RHS function does not depend on t, but our method still needs to have a placeholder for that variable.  Now run Heun’s method.  Note we will approximate the solution on the interval 0⩽t⩽4.

mygapprox=Heun[myg,0.0,{-1.0,2.0},4,50];  (* Initial conditions are a vector *)

Here we suppress the output because of the number of steps (we don’t need to see the numbers).  A phase-plane graph of this approximation is

myggraphPhase=ListPlot[mygapprox[[2]],
Joined->True,PlotMarkers->{Automatic,Medium},
PlotStyle->{Blue,Thick},
AxesOrigin->{0,0},
AspectRatio->1,
AxesLabel->{"y","y'"},
PlotLegends->LineLegend[{Blue},{"n=50"}]
]

Let’s overlay this on a direction field...

dField2=VectorPlot[
{v2, v1 v2-v1^2}, {v1, -4, 3.5}, {v2, -6.5, 2.5},
Axes->True,
AxesLabel->{(Subscript[v, 1]), (Subscript[v, 2])},
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->Directive[Darker[Green],Thick,"Line"],
StreamScale->None
];  (* Note that if I don't use a semi-colon here, it will save the plot and still show the plot. *)
Show[dField2,myggraphPhase,
PlotLabel->Style["Direction Field for the System",Bold,16]
]

In order to graph the state space variables we will need to create lists of t vs (Subscript[v, 1]) and t vs (Subscript[v, 2]).  Note n was 50, so we have 51 data points.

typairs=Table[{mygapprox[[1,k]],mygapprox[[2,k,1]]},{k,1,51}];
typrimepairs=Table[{mygapprox[[1,k]],mygapprox[[2,k,2]]},{k,1,51}];
tyDoubPrimepairs=Table[{mygapprox[[1,k]],mygapprox[[2,k,1]] mygapprox[[2,k,2]] - mygapprox[[2,k,1]] mygapprox[[2,k,1]]},{k,1,51}];

Now the two state-space graphs:

myggraphState=ListPlot[{typairs,typrimepairs, tyDoubPrimepairs},
Joined->True,PlotMarkers->{Automatic,Medium},
PlotStyle->{Directive[Blue,Thick],Directive[Red,Thick], Directive[Darker[Green],Thick]},
AxesOrigin->{0,0},
AxesLabel->{"t","y, y', and y''"},
PlotLegends->LineLegend[{Blue,Red, Darker[Green]},{"y","y'", "y''"}]
]