Chapter 6: Methods for solving ordinary differential equations

6.1 General

Problems involving ordinary differential equations (ODEs) are common in Physics. Take for example Newton's second law that can be written in the following differential form (using standard notation)

(((d^2)x)/(dt^2))=((∑F)/m)		(1)

This is a second-order differential equation that we may need to solve for the position x given the forces F and the initial conditions. We already saw examples of this such as a falling body without air resistance. In principle Equation (1) can be solved. However, Equation (1) is a second-order ODE. The numerical methods we will learn in this chapter are for solving first-order ODEs. There is a way around this problem: an ODE of any order can always be reduced to a set of first-order ODEs. For example, (1) can be rewritten as the following set of two first-order ODEs

(dv/dt)=((∑F)/m)		(2)	
	v=(dy/dt)			(3)

This principle can be generalized to higher order. In this course we will mostly deal with first and second-order ODEs. When dealing with a second-order ODE we will reduce it to a set of two first-order ODEs before applying the methods for solving first-order ODEs numerically. 
In general terms the objective of this chapter is to develop methods for solving the following general first-order ODE problem

(dy/dt)=f(t,y(t))	(4)

where f(t,y) can be any functions, from a simple constant to a more complicated nonlinear function. We will see examples applied to math and physics in this chapter and the following ones.

6.2 Euler Method

The simplest method for solving numerically Equation (4) is by discretizing it using a forward difference scheme for the derivative, i.e.

(((Subscript[y, n+1])-(Subscript[y, n]))/(Δt))=f((Subscript[t, n]),(Subscript[y, n]))	(5)

where the subscript n is an integer index that refers to each step of the solution. For example, the initial condition corresponds to the index n=1.
Re-arranging (5) we can get an expression for the solution at the future step n+1 as a function of all quantities known from the previous step n i.e.,

(Subscript[y, n+1])=Δt f((Subscript[t, n]),(Subscript[y, n]))+(Subscript[y, n])
where
	n=1,2,3,…,(Subscript[n, max])
and
	(Subscript[n, max])=((Duration of the simulation)/(Δt)).
The tabulated time is,
	(Subscript[t, n+1])=n Δt.

With a computer we can code this algorithm within a loop that automatically increment the index n=1,2,3,…,(Subscript[n, max]) and thus calculate discrete values of the solution of the problem, i.e. (Subscript[y, n]). Solving an equation like this with a computer is sometime referred to as carrying out a simulation. In practice, we cannot request the computer to compute the solution for n=1,2,3,…,∞. We need to stop at some point otherwise the calculation will never end. Therefore, when performing a simulation like this we need a priori to specify the duration of the simulation. The maximum number of iterations, i.e. (Subscript[n, max]), can then be calculated from the timestep Δt.

The choice of the time step Δt depends on the accuracy needed for a particular problem and the computational power available to do the computation. Some problems in physics are so complicated that only large timestep can be used for practical application. Ideally though you want the time step to be as small as possible for greater accuracy.
The Euler scheme is easy to understand and code. It is however not very accurate. The error of the scheme is of order Δt. This means that decreasing the timestep Δt by a factor of two decreases the error by a factor of two.

6.3 Euler step by step

Remove["Global`*"]

f[y_,t_]:=((t+1)/(y+1))

dt=0.15;
vp=VectorPlot[{dt,f[y,t]dt },{t,0,1.5},{y,0,4},VectorStyle->{"Segment"}]

sol=DSolve[{y'[t]==f[y[t],t],y[0]==y0},y[t],t]

Y[y0_,t_]=y[t]/.sol[[2]]

Show[vp,Plot[{Y[0.29,t],Y[1.425,t],Y[2.575,t]},{t,0,1.5},PlotStyle->{{Red,Thick},{Red,Thick},{Red,Thick}},PlotRange->All]]

y0=1;
t0=0;

f[y0,t0]//N

lin1[t_]:=y0+f[y0,t0]*(t-t0)

dt=0.3;
psol=Plot[Y[y0,t],{t,0,1.5}];
p1=Plot[lin1[t],{t,t0,t0+dt},PlotStyle->Red];

Show[psol,p1]

y1=lin1[t0+dt]
t1=t0+dt

lin2[t_]:=y1+f[y1,t1](t-t1)

p2=Plot[lin2[t],{t,t1,t1+dt},PlotStyle->Black];
Show[psol,p1,p2]

y2=lin2[t1+dt]
t2=t1+dt

lin3[t_]:=y2+f[y2,t2](t-t2)

p3=Plot[lin3[t],{t,t2,t2+dt},PlotStyle->Green];
Show[psol,p1,p2,p3]

Better to automate this.

t[1]=t0;
y[1]=y0;
dt=0.3;
tinit=0;
tfinal=1.5;

nmax=Round[((tfinal-tinit)/dt)]

Do[
t[i+1]=t[i]+dt;
y[i+1]=y[i]+f[y[i],t[i]]dt,{i,1,nmax}]

?t

?y

ytab=Table[{t[i],y[i]},{i,1,nmax+1}]

Show[psol,ListPlot[ytab,PlotStyle->{Red,PointSize[0.02]}],PlotRange->All]

6.4 Midpoint or Second-Order Runge-Kutta Method

The Euler scheme can be improved by using the midpoint scheme, also called the second-order Runge-Kutta scheme. The scheme is as follows (see your notes and comments given in class for more details on this scheme).

(Subscript[s, 1])=f((Subscript[t, n]),(Subscript[y, n]))
	(Subscript[y, n+1/2])=((Δt)/2) (Subscript[s, 1])+(Subscript[y, n])		(This is really just the Euler step)
	(Subscript[s, 2])= f((Subscript[t, n+1/2]),(Subscript[y, n+1/2]))	
	(Subscript[y, n+1])=Δt (Subscript[s, 2])+(Subscript[y, n])	
where
	n=1,2,3,…,(Subscript[n, max])
and
	(Subscript[n, max])=((Simulation Time)/(Δt)).
The tabulated time is,
	(Subscript[t, n+1])=n Δt.

The error of the scheme is of order ((Δt)^2). This means that decreasing the timestep Δt by a factor of  2 decreases the error by a factor of 4. It also means that compare to the Euler scheme and for the same time step used the accuracy of the scheme is better.

6.5 Fourth-Order Runge-Kutta Method

The most common scheme is the classical fourth-order Runge-Kutta formula (again see notes provided in class for more details).

(Subscript[s, 1])=f((Subscript[t, n]),(Subscript[y, n]))
	(Subscript[y, n+1/2])=((Δt)/2) (Subscript[s, 1])+(Subscript[y, n])		
	(Subscript[s, 2])= f((Subscript[t, n+1/2]),(Subscript[y, n+1/2]))
	(Subscript[y, n+1/2]^( * ))=((Δt)/2) (Subscript[s, 2])+(Subscript[y, n])	
	(Subscript[s, 3])=f((Subscript[t, n+1/2]),(Subscript[y, n+1/2]^( * )))
	(Subscript[y, n+1]^( * ))=Δt (Subscript[s, 3])+(Subscript[y, n])	
	(Subscript[s, 4])= f((Subscript[t, n+1]),(Subscript[y, n+1]^( * )))
	(Subscript[y, n+1])=((Δt)/6) ((Subscript[s, 1])+2(Subscript[s, 2])+2(Subscript[s, 3])+(Subscript[s, 4]))+(Subscript[y, n])	
where
	n=1,2,3,…,(Subscript[n, max])
and
	(Subscript[n, max])=((Simulation Time)/(Δt)).
The tabulated time is,
	(Subscript[t, n+1])=n Δt.

The error of the scheme is of order ((Δt)^4). This means that decreasing the timestep Δt by a factor of 2 decreases the error by a factor of 16. It also means that compare to the Euler or midpoint schemes the accuracy is typically much better.

6.6 Examples

Example 1

Solve numerically the following ODE

(dy/dt)=-y

with initial condition

y(0)=1

using the Euler, midpoint and fourth-order Runge-Kutta method, and compare graphically the numerical solutions to the exact solution. Take Δt=0.5 and carry out the simulation up to t=4.

Remove["Global`*"]

(* Set the parameters of the problem *)
T=4;
Δt=0.5;
nmax=Round[T/Δt];
yi=1; (* Initial condition *)

(* Find the exact solution *)
Exact[t_]=y[t]/.DSolve[{y'[t]==-y[t],y[0]==yi},y[t],t]//First;

(* EULER SCHEME *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Euler algorithm *)
Do[t[n+1]=n Δt;
y[n+1]=- Δt y[n]+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolEuler=Table[{t[n],y[n]},{n,1,nmax+1}];

SolEuler

(* MIDPOINT SCHEME *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Midpoint algorithm *)
Do[t[n+1]=n Δt;
s1=-y[n];
ytmp=((Δt)/2)s1+y[n];
s2=-ytmp;
y[n+1]=Δt s2+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolMidPoint=Table[{t[n],y[n]},{n,1,nmax+1}];

(* RUNGE-KUTTA *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Runge-Kutta algorithm *)
Do[t[n+1]=n Δt;
s1=-y[n];
ytmp=((Δt)/2)s1+y[n];
s2=-ytmp;
ytmp=((Δt)/2)s2+y[n];
s3=-ytmp;
ytmp=Δt s3+y[n];
s4=-ytmp;
y[n+1]=((Δt)/6)(s1+2 s2+2s3+s4)+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolRK4=Table[{t[n],y[n]},{n,1,nmax+1}];

(* Plot and compare the results of the three methods graphically *)

plot1=Plot[Exact[t],{t,0,T},AxesLabel->{"x","f(x)"}];
plot2=ListPlot[SolEuler,PlotStyle->RGBColor[0,0,1]];
plot3=ListPlot[SolMidPoint,PlotStyle->RGBColor[0,1,0]];
plot4=ListPlot[SolRK4,PlotStyle->RGBColor[1,0,0]];
Show[plot1,plot2,plot3,plot4,DisplayFunction->$DisplayFunction]

Blue dots: Euler
Green dots: Midpoint
Red dots: Runge-Kutta
Black line: Exact solution
It is clear that the Runge-Kutta algorithm is the most accurate.

Example 2

Solve numerically the following ODE

(dy/dt)=-2y+3(e^(-4t))

with initial condition

y(0)=1

using the Euler, midpoint and fourth-order Runge-Kutta method, and compare graphically the numerical solutions to the exact solution. Take Δt=0.2 and carry out the simulation up to t=3.

Remove["Global`*"]

(* Set the parameters of the problem *)
T=3;
Δt=0.2;
nmax=Round[T/Δt];
yi=1; (* initial condition *)

(* EXACT SOLUTION *)
Exact[t_]=y[t]/.DSolve[{y'[t]==-2y[t]+3(E^(-4t)),y[0]==yi},y[t],t]//First;

(* EULER SCHEME *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Euler algorithm *)
Do[t[n+1]=n Δt;
y[n+1]=Δt(-2y[n]+3(E^(-4t[n])))+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolEuler=Table[{t[n],y[n]},{n,1,nmax+1}];

(* MIDPOINT SCHEME *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Midpoint algorithm *)
Do[t[n+1]=n Δt;
s1=-2y[n]+3(E^(-4t[n]));
ytmp=((Δt)/2)s1+y[n];
s2=-2ytmp+3(E^(-4(t[n]+t[n+1])/2));
y[n+1]=Δt s2+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolMidPoint=Table[{t[n],y[n]},{n,1,nmax+1}];

(* RUNGE-KUTTA *)

(* Initial condition *)
y[1]=yi;
t[1]=0;

(* Runge-Kutta algorithm *)
Do[t[n+1]=n Δt;
s1=-2y[n]+3(E^(-4t[n]));
ytmp=((Δt)/2)s1+y[n];
s2=-2ytmp+3(E^(-4(t[n]+t[n+1])/2));
ytmp=((Δt)/2)s2+y[n];
s3=-2ytmp+3(E^(-4(t[n]+t[n+1])/2));
ytmp=Δt s3+y[n];
s4=-2ytmp+3(E^(-4t[n+1]));
y[n+1]=((Δt)/6)(s1+2 s2+2s3+s4)+y[n],
{n,1,nmax}]; 

(* Store the solution into a list for later analysis *)
SolRK4=Table[{t[n],y[n]},{n,1,nmax+1}];

(* Plot and compare the results of the three methods graphically *)

plot1=Plot[Exact[t],{t,0,T},AxesLabel->{"x","f(x)"}];
plot2=ListPlot[SolEuler,PlotStyle->RGBColor[0,0,1]];
plot3=ListPlot[SolMidPoint,PlotStyle->RGBColor[0,1,0]];
plot4=ListPlot[SolRK4,PlotStyle->RGBColor[1,0,0]];
Show[plot1,plot2,plot3,plot4,DisplayFunction->$DisplayFunction]

Blue dots: Euler
Green dots: Midpoint
Red dots: Runge-Kutta
Black line: Exact solution
It is clear that the Runge-Kutta algorithm is the most accurate.

Example 3 Writing reusable routines

So far, whenever we've implemented the midpoint scheme or 4th order RK, we've taken a template and rewritten the parts that contain the right-hand-sides of our ODE's.  Wouldn't it be nice if we could use the same code over and over again?

Let's solve the following ODE:

(dy/dt)+3y=cos(t)+5

with

y(0)=3,

from t=0 to t=10 using a timestep Δt=0.1.

Remove["Global`*"]

Define our right-hand-side function

f[y_,t_]:=Cos[t]+5-3 y

Define parameters and initial conditions

T=10;
Δt=0.1;
t[1]=0;
y[1]=3;

nmax=Round[T/Δt];

(* Midpoint scheme *)
Do[t[n+1]=n Δt+t[1];
s1=f[y[n],t[n]];
ytmp=((Δt)/2)s1+y[n];
s2=f[ytmp,t[n]+Δt/2];
y[n+1]=Δt s2+y[n],
{n,1,nmax}]

Interpolate our solution to get a smooth curve.

yM=Interpolation[Table[{t[n],y[n]},{n,1,nmax+1}]];

We can now get a function value for any value of t within initial and final times.

yM[2.3]

(* 4th order Runge-Kutta scheme *)
Do[t[n+1]=n Δt+t[1];
s1=f[y[n],t[n]];
ytmp=((Δt)/2)s1+y[n];
s2=f[ytmp,t[n]+Δt/2];
ytmp=((Δt)/2)s2+y[n];
s3=f[ytmp,t[n]+Δt/2];
ytmp=Δt s3+y[n];
s4=f[ytmp,t[n+1]];
y[n+1]=((Δt)/6) (s1 + 2 s2 + 2 s3 +s4)+y[n],

{n,1,nmax}]

Interpolate our solution to get a smooth curve.

yNR=Interpolation[Table[{t[n],y[n]},{n,1,nmax+1}]];

Compare on a single graph the exact solution with your numerical solutions.

yE[t_]=y[t]/.DSolve[{y'[t]==f[y[t],t],y[0]==3},y[t],t]//First

Plot[{yM[t],yNR[t],yE[t]},{t,0,T},PlotStyle->{Blue,Red}]

Plot the difference between your solutions and the exact solution as a function of time.

Plot[{yNR[t]-yE[t],yM[t]-yE[t]},{t,0,T},PlotStyle->{Blue,Red}]

The error for the 4th order R-K method is much smaller.

Example 4 A mechanics problem

Let's solve the following falling body problem:

m a=-m g - b v

with

y(0)=3m, and v(0)=20 m/s

from t=0 to t=4 using a timestep Δt=0.1, where b=0.1 and m=3.  First, break up the force equation into to ODE’s

	(dv/dt)=-g - (b/m)v
(dy/dt)=v

Remove["Global`*"]

Define our right-hand-side functions

fv[y_,v_,t_]:=-g -b/m v
fy[y_,v_,t_]:=v

Define parameters and initial conditions

T=4;
Δt=0.1;
t[1]=0;
y[1]=3;
v[1]=20;
g=9.81;
m=3;
b=0.1;

nmax=Round[T/Δt];

(* Midpoint scheme *)
Do[t[n+1]=n Δt+t[1];
s1v=fv[y[n],v[n],t[n]];
s1y=fy[y[n],v[n],t[n]];
vtmp=((Δt)/2)s1v+v[n];
ytmp=((Δt)/2)s1y+y[n];
s2v=fv[ytmp,vtmp,t[n]+Δt/2];
s2y=fy[ytmp,vtmp,t[n]+Δt/2];
v[n+1]=Δt s2v+v[n];
y[n+1]=Δt s2y+y[n],
{n,1,nmax}]

Interpolate our solution to get a smooth curve.

yM=Interpolation[Table[{t[n],y[n]},{n,1,nmax+1}]];

We can now get a function value for any value of t within initial and final times.

Plot[yM[t],{t,0,T}]