Lab 03 Runge-Kutta Methods

In the Beginning

As we have seen, depending on the function f, it is not always possible to solve initial value problems (IVP),

y'(t)=f(t,y), y((Subscript[t, 0]))=(Subscript[y, 0])

in a closed form.  In this case we often turn to numerical methods to approximate the solution to the IVP on some interval.  Euler’s method is one such numerical method, but Euler’s method is often not accurate enough to do a good job approximating solutions.  In this notebook we explore a family of numerical methods that include both low-order methods, like Euler’s method, and higher order (more accurate) methods.

First some background material.  In calculus you have studied Taylor polynomials:

f(x)=(Subscript[P, n])(x)+(Subscript[R, n])(x)

where

(Subscript[P, n])(x)=f((Subscript[x, 0]))+f'((Subscript[x, 0]))(x-(Subscript[x, 0]))+...+(f^((n)))((Subscript[x, 0])) ((((x-(Subscript[x, 0])))^n)/(n!)) = (∑_(k=0)^n)(f^((k)))((Subscript[x, 0]))((((x-(Subscript[x, 0])))^k)/(k!))

is the (n^th)-degree Taylor polynomial centered at (Subscript[x, 0]) with remainder

(Subscript[R, n])(x)=(f^((n+1)))(ξ) ((((x-(Subscript[x, 0])))^(n+1))/((n+1)!)),

where ξ is between x and (Subscript[x, 0]).  Here is a multivariable version of Taylor polynomials:

Theorem: Suppose that f(t,y) and all of its partial derivatives of order less than or equal to n+1 are continuous on D={(t,y):a⩽t⩽b, c⩽y⩽d} and let ((Subscript[t, 0]),(Subscript[y, 0]))∈D.  For every (t,y)∈D, there exists τ between t and (Subscript[t, 0]) and γ between y and (Subscript[y, 0]) with

f(t,y)=(Subscript[P, n])(t,y)+(Subscript[R, n])(t,y)

where

(Subscript[P, n])(t,y)=(∑_(k=0)^n)(∑_(j=0)^k)((((∂^k)f)/(∂(t^j)∂(y^(k-j))))((Subscript[t, 0]),(Subscript[y, 0])) ((((t-(Subscript[t, 0])))^j)/(j!)) ((((y-(Subscript[y, 0])))^(k-j))/((k-j)!)))

and

(Subscript[R, n])(t,y)=(∑_(k=0)^(n+1))((((∂^(n+1))f)/(∂(t^k)∂(y^(n+1-k))))(τ,γ) ((((t-(Subscript[t, 0])))^k)/(k!)) ((((y-(Subscript[y, 0])))^(n+1-k))/((n+1-k)!))) .

(Subscript[P, n])(t,y) is called the (n^th)-degree Taylor polynomial in two variables for the function f about ((Subscript[t, 0]),(Subscript[y, 0])), and (Subscript[R, n])(t,y) is the remainder term associated with (Subscript[P, n])(t,y).

Basic Runge-Kutta Methods

The idea behind the Runge-Kutta Methods is to use Taylor Series to find parameter values in a method so that the method matches the Taylor polynomial form of the solution to the IVP.  For instance, suppose that we have a method of the form y((Subscript[t, 0])+h)=(Subscript[y, 0])+a·f((Subscript[t, 0])+α,(Subscript[y, 0])+β), where (Subscript[y, 0])=y((Subscript[t, 0])) and we have parameters a, α, and  β.  We then expand the right-hand-side of the method using the theorem and match terms in the Taylor polynomial

y((Subscript[t, 0])+h)=(Subscript[y, 0])+h·y'((Subscript[t, 0]))+((h^2)/2)·y''((Subscript[t, 0]))+...+((h^n)/(n!))·(y^((n)))((Subscript[t, 0]))

Since y'(t)=f(t,y), we have y''(t)=(∂/(∂t))[f(t,y)]=(Subscript[f, t])(t,y)+(Subscript[f, y])(t,y)·y'(t).  Thus, the first three terms in the above Taylor polynomial are

y((Subscript[t, 0])+h)=(Subscript[y, 0])+h·f((Subscript[t, 0]),(Subscript[y, 0]))+((h^2)/2)·{(Subscript[f, t])((Subscript[t, 0]),(Subscript[y, 0]))+(Subscript[f, y])((Subscript[t, 0]),(Subscript[y, 0]))·f((Subscript[t, 0]),(Subscript[y, 0]))}.

Now expand the method in a Taylor Series:

y((Subscript[t, 0])+h)=(Subscript[y, 0])+a·{f((Subscript[t, 0]),(Subscript[y, 0]))+α·(Subscript[f, t])((Subscript[t, 0]),(Subscript[y, 0]))+β·(Subscript[f, y])((Subscript[t, 0]),(Subscript[y, 0]))+HOT}

Matching like terms gives us: a=h, a α=((h^2)/2) and a β=((h^2)/2)·f((Subscript[t, 0]),(Subscript[y, 0])).  Since we have three parameter, we are only able to match three terms to get three conditions on the parameters.  Also note that with α and β  terms of order h, then the higher order terms (HOT) are of order (h^2).  This gives us the midpoint method:

(Subscript[w, 0])=(Subscript[y, 0]) and (Subscript[w, j+1])=(Subscript[w, j])+h·f((Subscript[t, j])+(h/2),(Subscript[w, j])+(h/2)·f((Subscript[t, j]),(Subscript[w, j]))),

where (Subscript[t, j])=(Subscript[t, j-1])+h=(Subscript[t, 0])+j·h, and y((Subscript[t, j])) is approximated by (Subscript[w, j]) with local truncation error O((h^2)).

If we include more parameters, then we might be able to match some of the HOT and get a higher order method.  For instance when the method is of the form y((Subscript[t, 0])+h)=(Subscript[y, 0])+(Subscript[a, 1])·f((Subscript[t, 0]),(Subscript[y, 0]))+(Subscript[a, 2])·f((Subscript[t, 0])+α,(Subscript[y, 0])+β·f((Subscript[t, 0]),(Subscript[y, 0]))), then there is more flexibility in choosing the parameters to get an O((h^2)) method.  One important choice is when (Subscript[a, 1])=(Subscript[a, 2])=(h/2), α=β=h, or Modified Euler’s method:

(Subscript[w, 0])=(Subscript[y, 0]) and (Subscript[w, j+1])=(Subscript[w, j])+(h/2)(f((Subscript[t, j]),(Subscript[w, j]))+f((Subscript[t, j])+h,(Subscript[w, j])+h·f((Subscript[t, j]),(Subscript[w, j])))),

which has better numerical properties than Euler’s or Midpoint methods as h→0.  Note that this method can be written as:

(Subscript[w, 0])=(Subscript[y, 0]) and then for j=1, 2, 3, ..., n-1,
(Subscript[m, 1])=h·f((Subscript[t, j]),(Subscript[w, j])),
(Subscript[m, 2])=h·f((Subscript[t, j])+h,(Subscript[w, j])+(Subscript[m, 1])),
(Subscript[w, j+1])=(Subscript[w, j])+(1/2)((Subscript[m, 1])+(Subscript[m, 2])).

For a third-order method we use a term that looks like f(t+(Subscript[α, 1]),y+(Subscript[β, 1])·f(t+(Subscript[α, 2]),y+(Subscript[β, 2])·f(t,y))).  The expansion using Taylor Series gets a bit scary for this one, but luckily some poor fool has already done it to get Heun’s method:

(Subscript[w, 0])=(Subscript[y, 0]) and then for j=1, 2, 3, ..., n-1,
(Subscript[m, 1])=h·f((Subscript[t, j]),(Subscript[w, j])),
(Subscript[m, 2])=h·f((Subscript[t, j])+(h/3),(Subscript[w, j])+(1/3)(Subscript[m, 1])),
(Subscript[m, 3])=h·f((Subscript[t, j])+((2h)/3),(Subscript[w, j])+(2/3)(Subscript[m, 2])),
(Subscript[w, j+1])=(Subscript[w, j])+(1/4)((Subscript[m, 1])+3(Subscript[m, 3]))

Implementation

The following module implements 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}]
]

The following is how you use the Module to approximate the solution to the IVP: y'=y-(t^2)+1, y(0)=0.5 on the interval 0⩽t⩽2 with 10 subintervals.

myf[t_,y_]=y-(t^2)+1

myapprox10=Heun[myf,0.0,0.5,2.0,10]

This generated a list with two elements: the list of t-values and the list of corresponding approximation values.  In order to graph the approximate solution we will need to create a list of ordered pairs.  Do this with the Table command, remembering that we had 10 steps (so 11 points).

mypairs=Table[{myapprox10[[1,k]],myapprox10[[2,k]]},{k,1,11}]

Now we can plot the approximation to the solution of the IVP:

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

Let’s compare this to the direction field for this DE.

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 lets try an example using a system of equations.  Suppose we wish to numerically approximate the IVP y''+y=0, with y(0)=1 and y'(0)=0.  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:

(Subscript[v, 1])'=(Subscript[v, 2]),   (Subscript[v, 1])(0)=1,
(Subscript[v, 2])'=-(Subscript[v, 1]), (Subscript[v, 2])(0)=0.

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

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 that the solution to this IVP is y(t)=cos(t) so we will approximate the solution on the interval 0⩽t⩽2π.

mygapprox=Heun[myg,0.0,{1.0,0.0},2*Pi,20];

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=20"}]
]

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 20, so we have 21 data points.

typairs=Table[{mygapprox[[1,k]],mygapprox[[2,k,1]]},{k,1,21}];
typrimepairs=Table[{mygapprox[[1,k]],mygapprox[[2,k,2]]},{k,1,21}];

Now the two state-space graphs:

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

One last Comment:  Since we know what the actual solution is for y, we can calculate the actual error in this approximation at each time-step.

mygerrorPairs=Table[
{mygapprox[[1,k]],Abs[Cos[mygapprox[[1,k]]]-mygapprox[[2,k,1]]]},
{k,1,21}];

Graphing this gives:

mygerrorPlot=ListPlot[mygerrorPairs,
Joined->True,PlotMarkers->{Automatic,Medium},
PlotStyle->{Darker[Green],Thick},
AxesOrigin->{0,0},
AxesLabel->{"t","abs error"},
PlotLegends->LineLegend[{Darker[Green]},{"error"}]
]

Assignment

The following is a forth-order Runge-Kutta method commonly used in mathematical software:

(Subscript[w, 0])=(Subscript[y, 0]) and for j=1, 2, 3, ..., n-1
(Subscript[m, 1])=h·f((Subscript[t, j]),(Subscript[w, j])),
(Subscript[m, 2])=h·f((Subscript[t, j])+(h/2),(Subscript[w, j])+(1/2)(Subscript[m, 1])),
(Subscript[m, 3])=h·f((Subscript[t, j])+(h/2),(Subscript[w, j])+(1/2)(Subscript[m, 2])),
(Subscript[m, 4])=h·f((Subscript[t, j])+h,(Subscript[w, j])+(Subscript[m, 3])),
(Subscript[w, j+1])=(Subscript[w, j])+(1/6)((Subscript[m, 1])+2(Subscript[m, 2])+2(Subscript[m, 3])+(Subscript[m, 4])).

Write a Module that implements this forth-order Runge-Kutta method.  Your Module should look similar to the above Heun’s method.  Put comments in the Module similar to those in Heun’s method.  You may discuss your work with others, but do not plagiarize others.

For the following IVP: y'=Exp[t]*Cos[y], y(0)=-(3/2), use the forth-order Runge-Kutta method to approximate the solution on the interval 0⩽t⩽4, using n=4, 8, 16, 32.  Graph each approximation separately along with the direction field for this IVP.  What can be said about this approximation method as the number of subintervals increases (step size decreases).  Use your own words, don’t plagiarize others!

How do the results in problem 2 compare with approximations using Euler’s method.  Use your own words, don’t plagiarize others!

For the following IVP: x''-(1-(x^2))x'+x=0, x(0)=0, x'(0)=1, convert to a system of equations and use the forth-order Runge-Kutta method to approximate a solution on the interval 0⩽t⩽50, with step sizes h=0.5, 0.1, 0.01 (Yes, you will need to figure out the corresponding number of subintervals to use to get these step sizes).  Generate both state-space and phase-space graphs (that look pretty) for these approximations (you may want to suppress the output of the lists when the number of steps is large).  Describe the behavior of the solution to this IVP.  Use your own words, don’t plagiarize others!!