Chapter 2
Numerical methods for solving
ODEs
We will study two methods for finding approximate solutions of ODEs. Such methods may be used for
(at least) two reasons
the ODE does not have an exact solution or, the solution is difficult to find
even if we know or could find the exact solution we are only interested in obtaining a approximate
solution (for example to estimate the solutions at a given time or to draw a plot of the solution.)
The methods to be considered are both easy to use but are not as fast or as reliably accurate as
the methods most often used on a computer to solve real problems. Nevertheless they illustrate the
important aspects of the techniques involved in numerical integration of ODEs.
2.1
Eulers method
Consider the initial value problem
dy
= f (t, y), y(a) = A.
(2.1)
dt
We will assume that this has a unique solution y(t). The value of this solution at t = a + h may be
approximated using Taylors theorem to obtain
dy
(a),
dt
and then using the ODE (2.1) to express the derivative of y in terms of y and t,
y(a + h) y(a) + h
y(a + h) y(a) + hf (a, y(a)).
If we define t0 := a, y0 := y(a), t1 := a + h, we can express the above result in the following way;
y1 := y0 + hf (t0 , y0 )
gives an approximation to y(t1 ). This method of approximating the solution to an initial value problem
is usually called Eulers method. The absolute error (or sometimes just error ) and the relative error
r in this approximation are
|y(t1 ) y1 |
= |y(t1 ) y1 |, r =
.
y(t1 )
Figure 2.1 gives a geometrical interpretation of Eulers method in which the solution between t = t0
and t = t1 is approximated by the tangent to the solution curve at (t0 , y0 ).
We wish to repeat this process through a number of steps in time and to describe this it is useful to
extend the above notation;
t2 := t1 + h,
t3 := t2 + h,
..
.
y2 := y1 + hf (t1 , y1 )
y3 := y2 + hf (t2 , y2 )
..
.
18
Figure 2.1: Geometrical interpretation of Eulers method
that is,
t0 := a, ti+1 := ti + h for i = 0, 1, 2, . . .
y0 := y(t0 ), yi+1 := yi + hf (ti , yi ) for i = 0, 1, 2, . . . .
(2.2)
(2.3)
Here yi gives an approximation to y(ti ) and the lines joining (t0 , y0 ) to (t1 , y1 ) to to (tn , yn ) gives an
approximation to the solution curve y as a function of t on the time interval [t0 , tn ], that is [a, a + nh],
as shown in Figure 2.2.
Figure 2.2: Approximate solution found by Eulers method
Eulers method algorithm To approximate the solution to the initial value problem (2.1) on the
interval [a, b] Eulers method expressed in algorithmic form as follows
Pseudo code
Set t=a and y=A
Store t and y
Repeat while t < b
{
Set F=f(t,y), t=t+h and y=y+hF
Store t and y
}
Maple code
t:=a; y:=A;
ans:=[t,y];
while t<b
do
F:=f(t,y); t:=t+h; y:=y+h*F;
ans:=ans,[t,y];
od;
19
This algorithm is illustrated in the following example.
Example 2.1 Use Eulers method with h = 1.0, h = 0.5 and h = 0.25 to approximate y(1) where y(t)
is the solution of the initial value problem
dy
= y + t,
dt
y(0) = 1.
Show that the exact solution of this IVP is y(t) = 2et t 1 and find the relative error in y(1) for each
step length.
Solution We use the algorithm described above with F = f (t, y) = y + t and work to six significant
figures throughout.
h = 1.0
F
1.0
t
0.0
1.0
y
1.0
2.0
Therefore h = 1 give the approximation y1 = 2.0 to y(1).
h = 0.5
F
1.0
2.0
t
0.0
0.5
1.0
y
1.0
1.5
2.5
Therefore h = 0.5 give the approximation y2 = 2.5.
h = 0.25
F
1.0
1.5
2.125
2.90625
t
0.0
0.25
0.5
0.75
1.0
y
1.0
1.25
1.625
2.15625
2.88281
Therefore h = 0.25 give the approximation y4 = 2.88281.
For y(t) = 2et t 1,
dy
= 2et 1 = y + t,
dt
y(0) = 2e0 0 1 = 1,
therefore y(t) = 2et t1 is the solution of the IVP. The exact answer is therefore y(1) = 2e2 = 3.43656.
The relative errors for the different values of h are
h
1.0
0.5
0.25
r
0.418023
0.272528
0.161135
Observe in the above table that as the step size h is halved, the relative error is also halved. For
smaller h this phenomenon persists
h
0.1
0.01
0.001
r
0.0724795
0.00783632
0.000788579
20
and we see that the relative error is of the same order as the step size. In fact this relationship between
error size and step length is true in general and we write
r = O(h1 ),
and say that Eulers method is an order 1 method.
The Maple code to carry out Eulers method in the case h = 0.25 is shown below;
h:=0.25;
f:=(t,y)->t+y;
t:=0; y:=1;
ans:=[t,y];
while t<1
do
F:=f(t,y); t:=t+h; y:=y+h*F;
ans:=ans,[t,y];
od:
After this, ans stores the values of t and y as the sequence
[0,1],[0.25,1.25],[0.5,1.6250],[0.75,2.156250],[1,2.88281250].
To plot the corresponding approximate solution curve, the following Maple code may be used
with(plots);
pointplot([ans],style=line);
and
with(plots);
display({plot(2*exp(t)-t-1,t=0..1),pointplot([ans],style=line),pointplot([ans])});
which plots the exact solution, this curve and the points (ti , yi ) on the same axes. The resulting graph
is shown in Figure 2.3.
2.2
The Euler predictor-corrector method
This method is also sometimes called the improved Eulers method. The aim of this method is the same
as Eulers method, to approximate y(t1 ) = y(a + h) for the solution of the IVP (2.1). First, as in Eulers
method the tangent to the solution through (t0 , y0 ), of gradient f (t0 , y0 ), is used predict an approximation
z1 to y(t1 ), that is
t1 := t0 + h, z1 := y0 + hf (t0 , y0 ).
The gradient of the solution curve through (t1 , z1 ) is then used to correct the first approximation by
using the line through (t0 , y0 ) whose gradient is the average of f (t0 , y0 ) and f (t1 , z1 )
y1 := y0 +
h
(f (t0 , y0 ) + f (t1 , z1 )).
2
This is illustrated in Figure 2.4
Euler predictor-corrector method algorithm To approximate the solution to the initial value
problem (2.1) on the interval [a, b] the Euler predictor corrector method expressed in algorithmic form
as follows
Pseudo code
Set t=a and y=A
Store t and y
Repeat while t < b
{
Set F=f(t,y), t=t+h and z=y+hF
Set Z=f(t,z), y=y+h(F+Z)/2
Store t and y
}
Maple code
t:=a; y:=A;
ans:=[t,y];
while t<b
do
F:=f(t,y); t:=t+h; z:=y+h*F;
Z:=f(t,z); y:=y+h*(F+Z)/2;
ans:=ans,[t,y];
od;
21
Euler's method and exact solution
2.5
1.5
1
0
0.2
0.4
0.6
0.8
Figure 2.3: Eulers method and exact solution in Maple
Example 2.2 Repeat the numerical parts of Example 2.1 using the Euler predictor-corrector method
rather than Eulers method. Compare the relative errors for the two methods for the different step
lengths.
Solution Here F = t + y and Z = t + z. In evaluating F and Z we always use the most recent values
of the arguments; for F the values of t and y used are on the previous row but for Z they are on the
same row.
h = 1.0
F
1.0
t
0.0
1.0
z
2.0
h = 0.5
22
Z
3.0
y
1.0
3.0
Figure 2.4: The Euler predictor-corrector method
F
1.0
2.25
t
0.0
0.5
1.0
z
1.5
2.875
Z
2.0
3.875
y
1.0
1.75
3.28125
h = 0.25
F
1.0
1.5625
2.28320
3.20660
t
0.0
0.25
0.50
0.75
1.00
z
1.25
1.70312
2.35400
3.25826
Z
1.50
2.20312
3.10400
4.25826
y
1.0
1.3125
1.78320
2.45660
3.38971
Here we tabulate the relative errors for the two methods and different step lengths
h
1.0
0.5
0.25
Euler r
0.418023
0.272528
0.161135
Euler p-c r
0.127034
0.0451934
0.0136328
Using Maple we can also calculate approximations for smaller step lengths, together with the corresponding relative errors
h
0.1
0.01
0.001
Euler r
0.072479
0.007836
0.000789
Euler p-c r
0.00244430
0.00002616
0.00000026
From the above tables it is clear that the Euler predictor-corrector relative error is proportional to
the square of the step length, in the example at least. In fact, it may be shown that this is true generally
for this method and so
r = O(h2 ),
i.e. Euler predictor-corrector is an order 2 method.
23