SMSTC (2024/25)

Numerical Methods

www.smstc.ac.uk

Contents

1 Numerical Methods for ODEs 1 1–1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–2
1.1.1 The Standard ODE Initial Value Problem (IVP) . . . . . . . . . . . . . . . . . . . 1–2
1.1.2 Notation and general ideas about approximating ODEs . . . . . . . . . . . . . . . 1–3
1.1.3 Newton’s method for systems of algebraic equations . . . . . . . . . . . . . . . . . 1–3
1.2 Euler’s Method and variants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–4
1.2.1 Derivations of Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–5
1.2.2 Extensions of Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–5
1.2.3 Errors and Convergence of the Euler Scheme . . . . . . . . . . . . . . . . . . . . . 1–7
1.3 Runge-Kutta (RK) methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–8
1.3.1 Explicit RK methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–8
1.3.2 Implicit RK methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–9
1.4 Linear Multistep Methods — LMMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–10
1.4.1 General form of LMMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–10
1.4.2 Examples of LMM Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–11
1.5 Desired and essential properties of schemes . . . . . . . . . . . . . . . . . . . . . . . . . . 1–11
1.5.1 Zero-stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–12
1.5.2 Absolute stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–12
1.5.3 Local Truncation Error (LTE) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–13
1.5.4 Consistency of scheme and initial data . . . . . . . . . . . . . . . . . . . . . . . . . 1–15
1.5.5 Consistency of initial data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–15
1.5.6 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–15

2 Numerical Methods for ODEs 2 2–1

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2 Systems of ODEs and stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2.1 Linear systems of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2.2 Nonlinear systems of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
2.2.3 Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
2.2.4 A-stability and A(α)-stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–4
2.3 Computational effort and efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–5
2.4 Practical global error estimation for fixed step schemes . . . . . . . . . . . . . . . . . . . . 2–5
2.5 Error Estimation and Time Step Adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . 2–7
2.5.1 Local error estimation for RK methods . . . . . . . . . . . . . . . . . . . . . . . . 2–8
2.5.2 Local error estimation for implicit LMMs . . . . . . . . . . . . . . . . . . . . . . . 2–8
2.5.3 Time step control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–9
2.5.4 How is the global error related to the local error tolerance? . . . . . . . . . . . . . 2–10
2.6 What solvers are available and how do we use them? . . . . . . . . . . . . . . . . . . . . . 2–10

0–1
SMST C: Numerical Methods (i)

2.6.1 Explicit solvers for non-stiff problems . . . . . . . . . . . . . . . . . . . . . . . . . 2–11

2.6.2 Implicit LMM solvers for stiff problems . . . . . . . . . . . . . . . . . . . . . . . . 2–11
2.6.3 Other sources of ODE software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–13
2.7 The Method of Lines for simple, time-dependent PDEs . . . . . . . . . . . . . . . . . . . . 2–16
2.7.1 Output and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–18
2.8 What has not been covered? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–19
 SMSTC (2024/25)
Numerical Methods
Lecture 1: Numerical Methods for ODEs 1
Dugald B Duncan, Heriot-Watt Universitya

www.smstc.ac.uk

Contents
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–2
1.1.1 The Standard ODE Initial Value Problem (IVP) . . . . . . . . . . . . . . . . . 1–2
1.1.2 Notation and general ideas about approximating ODEs . . . . . . . . . . . . . 1–3
1.1.3 Newton’s method for systems of algebraic equations . . . . . . . . . . . . . . . 1–3
1.2 Euler’s Method and variants . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–4
1.2.1 Derivations of Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–5
1.2.2 Extensions of Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–5
1.2.3 Errors and Convergence of the Euler Scheme . . . . . . . . . . . . . . . . . . . 1–7
1.3 Runge-Kutta (RK) methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–8
1.3.1 Explicit RK methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–8
1.3.2 Implicit RK methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–9
1.4 Linear Multistep Methods — LMMs . . . . . . . . . . . . . . . . . . . . . . . 1–10
1.4.1 General form of LMMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–10
1.4.2 Examples of LMM Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–11
1.5 Desired and essential properties of schemes . . . . . . . . . . . . . . . . . . 1–11
1.5.1 Zero-stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–12
1.5.2 Absolute stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–12
1.5.3 Local Truncation Error (LTE) . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–13
1.5.4 Consistency of scheme and initial data . . . . . . . . . . . . . . . . . . . . . . . 1–15
1.5.5 Consistency of initial data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–15
1.5.6 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1–15

a D.B.Duncan@hw.ac.uk

1–1
SMST C: Numerical Methods 1–2

1.1 Introduction
These two lectures are about the accurate and efficient approximation of the solution of systems of
ordinary differential equations - hereafter abbreviated to “ODEs”. The aim is to introduce key ideas and
the language used to describe the attributes of numerical methods, as opposed to doing rigorous analysis.
There are numerous references describing methods for solving ODEs and related problems. The book [1]
gives an approachable introduction, in particular for stiff problems and DAEs, while [5] is also very
accessible, with interesting examples and final sections on geometric integration. The book [6] gives a
general introduction to numerical ODEs and PDEs, as do the more general books on numerical analysis
such as [2, 3]. For a modern view of numerical methods from a dynamical systems point of view the book
by [7] is a good source.
There is usually no need to write your own solver for ODEs. There are lots of good ODE software
packages available that many people have devoted years of their working life to. You should first see
what software is already available, and only consider writing a solver yourself if no software is suitable.
However it is useful to have a good idea of the theory and ideas that underlie solving ODEs and hence
understand the limitations and advantages of packages.

1.1.1 The Standard ODE Initial Value Problem (IVP)

Most numerical analysis texts and general purpose computer codes deal with the following, standard
ODE initial value problem (IVP). It consists of a system of M ≥ 1 first order ODEs for M solution
components, coupled with the initial condition on the M solution components: find u(t) ∈ RM satisfying

du
= f (t, u), given u(0) = u0 . (1.1)
dt
Usually we want the solution over a fixed, finite time interval 0 ≤ t ≤ T . The ODE above is linear if it
can be written in the form
du
= A(t)u + b(t)
dt
where A ∈ RM ×M and b ∈ RM are functions of t only. Otherwise it is nonlinear.
Different notations for derivatives, e.g.

du d2 u
ut ≡ u′ ≡ u̇ ≡ u(1) ≡ , utt ≡ u′′ ≡ ü ≡ u(2) ≡ ,
dt dt2
are used in different settings.
Note that we can write many higher order ODEs in this form — see the example below. Also, systems of
complex valued ODEs can be written as a standard real system of twice the dimension by equating real
and imaginary parts of each complex ODE. The following example illustrates the conversion of a real,
second order ODE into a two component system of first order ODEs.

Example 1.1 (General forced Van der Pol equation) Consider the general forced Van der Pol equa-
tion
x′′ + ϕ(x)x′ + x = p(t), x(0) = a, x′ (0) = b.
Here x(t) is the current in an electronic circuit, x′ its rate of change with time t, p(t) is a known forcing
function and ϕ(x) is a function determined by the properties of the circuit. Now define u1 (t) = x(t) and
u2 (t) = x′ (t), then this second order ODE is equivalent to the first order system

u′1 = u2 , u1 (0) = a
u′2 = p(t) − ϕ(u1 )u2 − u2 , u2 (0) = b.
T
Setting u = (u1 , u2 ) and  
u2
f (t, u) =
p(t) − ϕ(u1 )u2 − u1
we get our standard ODE IVP (1.1).
SMST C: Numerical Methods 1–3

When trying to solve any mathematical problem it is worth reviewing if there is any information to tell
you that the solution actually exists, if it is unique, if it has any special properties that might be useful
or troublesome etc.. We are not going to go into the theory of existence and uniqueness for ODEs, but
you can find that in standard texts. However, note that even for continuous f it is not true that (1.1)
has a unique solution – see the next example.

Example 1.2 Consider the simple ODE IVP u′ = u1/2 , u(0) = 0. This has two solutions u(t) ≡ 0 and
u(t) = t2 /4, and hence the solution is not unique.

1.1.2 Notation and general ideas about approximating ODEs

To construct time-stepping numerical approximation methods for a well posed system of ODEs for t ∈
[0, T ] we start by discretizing (splitting) the interval [0, T ] into small pieces, called a mesh or grid,
and delineated by time levels (or grid points or mesh points) labelled tn for n = 0, 1, 2, . . . , N . The
approximate solution at these time levels is

un ≈ u(tn ) n = 0, 1, 2, . . . , N

and we will examine methods to find these approximations. Interpolation may be used to fill in the gaps
in the approximate solution between these time levels, although often we simply choose the time levels
to give output at the times we are interested in.
To make life simple for now we take constant steps of size ∆t between time levels so that

tn = n∆t, tn+1 − tn = ∆t, ∆t = T /N > 0.

In the next lecture we consider how to change the step size as the calculation progresses from t = 0 to
t = T with the goal of making the overall process as efficient as possible.
The usual aim is to obtain an approximate solution un that is close enough to the exact solution. For
a given numerical method the only control we have is to choose the timestep size ∆t. For any sensible
method we want the difference between exact and approximate solutions (we will come on to how that is
measured) to → 0 as the step size ∆t → 0. Usually we want the approximate solution over a fixed time
interval t ∈ [0, T ], and so

smaller ∆t ⇒ more steps to get from t = 0 to t = T

and hence more work.

1.1.3 Newton’s method for systems of algebraic equations

Many of the ODE approximations we consider below are said to be implicit since in general they require
the solution of a nonlinear system of equations at each time step. There are many ways to do this, and
some are specifically designed for ODE solvers. Newton’s method is one of the most common, and it is
useful to see it before we proceed.
Given the system of M algebraic equations
g(x) = 0
in M unknowns x, Newton’s method consists of the following.

ˆ Make an initial guess x[0] at the solution.

ˆ For ν = 0, 1, . . . do
 −1
∂g [ν]
x[ν+1] = x[ν] − (x ) g(x[ν] ) (1.2)
∂x
and stop when ∥x[ν+1] − x[ν] ∥ is small enough.
ˆ Then x[ν+1] approximates the solution (root) of g(x) = 0.
SMST C: Numerical Methods 1–4

Exact
1 ∆t=1
∆ t = 0.5
∆ t = 0.25

Solution 0.8

0.6

0.4

0.2

0
0 1 2 3 4 5 6 7 8 9 10
time t

Figure 1.1: Explicit Euler method approximation of u′ = u(1 − u) with initial condition u(0) = 0.01 and
various stepsizes ∆t. The symbols show the computed values and linear interpolation (“join the dots”)
is used to fill in the detail between.

Here ∂g/∂x is the M × M Jacobian matrix with elements ∂gi /∂xj . Note that if g is linear in x, then
Newton’s method gives the exact solution in one iteration.
Of course one should usually not compute the matrix inverse, and the iteration (1.2) is often done in two
steps: solve the linear system of equations
 
∂g [ν]
(x ) δ [ν+1] = −g(x[ν] )
∂x
for δ [ν+1] ; and set
x[ν+1] = x[ν] + δ [ν+1] .

1.2 Euler’s Method and variants

The explicit Euler Method is the most basic of numerical ODE schemes. The scheme is
un+1 = un + ∆t f (tn , un ) , n = 0, 1, . . . , N − 1, (1.3)
0 n+1
given the initial data u = u0 . The method is said to be explicit since we get the solution u at the
new time level directly by evaluating a formula without solving any algebraic equations. It is generally
too inefficient and limited for massive calculations, but it is good for motivating other schemes and some
of the basic concepts.

Example 1.3 Consider the approximate solution of the logistic equation IVP u′ = u(1 − u) with u(0) =
0.2 by the explicit Euler method, using steps of size ∆t = 0.1.
We have
u1 = u0 + ∆t u0 (1 − u0 ) = 0.2 + 0.1 × 0.2 × 0.8 = 0.216
u2 = u1 + ∆t u1 (1 − u1 ) = 0.216 + 0.1 × 0.216 × 0.784 = 0.2329344
...
u10 = u9 + ∆t u9 (1 − u9 ) = ··· = 0.399696 . . .
SMST C: Numerical Methods 1–5

Hence u(1) ≈ u10 = 0.399696 . . .. (The exact solution is u(1) = 0.40460967519 . . ..)

1.2.1 Derivations of Euler’s method

We now derive (1.3) in two different ways: firstly using Taylor’s theorem; and secondly by approximating
an integral of the ODE.

Derivation of Euler using Taylor’s theorem

Suppose the solution u(t) has 2 continuous derivatives on [0, T ] and then apply Taylor’s Theorem so that

(tn+1 − tn )2 ′′ n
u(tn+1 ) = u(tn ) + tn+1 − tn u′ (tn ) +


u (ξ )
2
for some ξ n in the interval (tn , tn+1 ). Since ∆t = tn+1 − tn and u′ (t) satisfies the differential equation
(1.1) we get
∆t2 ′′ n
u(tn+1 ) = u(tn ) + ∆t f (tn , u(tn )) + u (ξ ). (1.4)
2
Euler’s method is obtained by neglecting the remainder term and replacing u(tn ) and u(tn+1 ) by their
approximations un and un+1 .

Derivation of Euler by approximate integration

An alternative derivation of the explicit Euler method is found by integrating the ODE between tn and
tn+1 to get
Z tn+1 Z tn+1 Z tn+1
′ n+1 n
u (s)ds = f (s, u(s))ds ⇔ u(t ) = u(t ) + f (s, u(s)) ds (1.5)
tn tn tn

and then applying the quadrature rule (approximate integration)

Z b
g(s)ds ≈ (b − a)g(a). (1.6)
a

This gives
Z tn+1
n+1 n
u(t ) = u(t ) + f (s, u(s)) ds ≈ u(tn ) + (tn+1 − tn )f (tn , u(tn ))
tn

and Euler’s method follows as for the last stage of the Taylor approach.
The approach of approximating the integral is illustrated further below and is particularly useful in
constructing schemes for stochastic differential equations later in this stream.

1.2.2 Extensions of Euler

Let us take the integral approximation idea above a little further. For example, if we replace (1.6) by the
one parameter family of integration rules
Z tn+1
g(s) ds ≈ ∆t (1 − θ) g(tn ) + θ g(tn+1 ) ,


tn

where constant θ ∈ [0, 1] we get

u(tn+1 ) ≈ u(tn ) + ∆t (1 − θ)f (tn , u(tn )) + θf (tn+1 , u(tn+1 )) ,

n = 0, 1, , N − 1.

This leads to the θ-method

un+1 = un + ∆t (1 − θ)f (tn , un ) + θf (tn+1 , un+1 ) ,

n = 0, 1, , N − 1. (1.7)

With θ = 0 we recover the explicit Euler method and with θ = 1 we get the Implicit Euler method

un+1 = un + ∆t f (tn+1 , un+1 ) (1.8)

SMST C: Numerical Methods 1–6

which is widely used in approximating the ODE systems that arise in the numerical solution of parabolic
PDEs. Another commonly used scheme in PDEs and elsewhere has θ = 1/2, giving

∆t
un+1 = un + f (tn , un ) + f (tn+1 , un+1 ) ,


(1.9)
2
and is called the Trapezoidal rule.
The Trapezoidal rule and implicit Euler method are examples of implicit schemes. In order to get un+1
we need to solve a system of algebraic equations

Implicit Euler: un+1 − ∆t f (tn+1 , un+1 ) − |{z}

un = 0,
known

∆t n+1 n+1 ∆t n n
Trapezoidal Rule: un+1 − f (t ,u ) − un − f (t , u ) = 0
2 | 2{z }
all known
n+1 n+1
for the unknown u . Both have the form g(u ) = 0 ready to fit directly into Newton’s method (1.2).
The initial guess un+1,[0] may be obtained in various ways, e.g.

un+1,[0] = un or un+1,[0] = un + ∆t f (tn , un ),

using explicit Euler as a predictor of the value in the second case. These algebraic equations are linear
when f is linear in u and nonlinear in the general case that f is nonlinear in u. Implicit schemes clearly
need much more work for each step than explicit schemes, but we will see that they often allow larger
time steps to be used and for a wide class of ODEs are more efficient.
A final Euler-related method is obtained by approximating an approximation! Instead of solving the
Trapezoidal rule algebraic equations exactly for un+1 , use Euler to approximate un+1 on the right hand
side of (1.9) to get
 
∆t  n n
un+1 = un + n n n n 
f (t , u ) + f (t + ∆t, u + ∆tf (t , u )) . (1.10)
2 | {z }
≈un+1

This is called the Improved Euler scheme and is a Runge-Kutta method. It is usually implemented
and written in the equivalent Runge-Kutta format (1.14) that we discuss later. Notice that it involves
a term of the form f (t + ∆t, u + ∆t f (t, u)) – i.e. composition of f with itself rather than just a linear
combination of f ’s.
There are various natural ways to generalise the explicit Euler method beyond what we have seen above.
One approach is to use information from more than one previous time level to get the approximation un+1
at tn+1 . This results in Linear Multistep Methods (LMMs) which we examine in Section 1.4. The
other main approach is like the Improved Euler scheme above, essentially to to use function evaluations
between time levels tn and tn+1 as described in Section 1.3 on Runge–Kutta Methods.
SMST C: Numerical Methods 1–7

0
Max Error for t [0,10]
Error at t=10
-1 Slope 1

log10(Absolute Error) -2

-3

-4

-5

-6

-7
-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0
log10(Step size t)
Figure 1.2: Error in the explicit Euler method approximation of u′ = u(1 − u) with initial condition
u(0) = 0.01 and various stepsizes ∆t. Both the maximum error over the time interval and the error at
t = 10 appear to be proportional to ∆t as ∆t → 0.

1.2.3 Errors and Convergence of the Euler Scheme

We want to evaluate the error in approximations of the ODE (1.1). To do this we need to think about
how to compare the continuous solution from the ODE with the approximations given at a set of discrete
points in time that we get from approximations like Euler. We either need to interpolate the discrete
points or simply compare at the grid points. The standard choice is to compare the ODE and difference
scheme at the grid points.
We define the global error at time tn by

en = u(tn ) − un , n = 0, 1, 2, . . . , N.

In other words the difference between exact and approximate solutions. Ideally we want to show that
as ∆t is reduced the global error is also reduced. This certainly appears to be happening in the Euler
example in Figure 1.1, where the approximation gets better as ∆t is reduced. Figure 1.2 evaluates this
more carefully, and we see that the global error in the Euler approximation on the fixed time interval
appears to behave like
Global Error in explicit Euler = O(∆t)
as ∆t → 0 in this example. To verify this, work out the slope of the lines as ∆t → 0 (i.e. moving to the
left on the plot).
A useful step to proving convergence is to look at the local truncation error (LTE) described in Section
1.5.3. For explicit Euler this is

u(tn+1 ) − u(tn )
LTE(tn ) = − f (tn , u(tn )) = O(∆t)
∆t
where u(t) is an exact solution of the ODE u′ = f (t, u), and the end result is found by Taylor expansion
about t = tn . The LTE gives the order of the method. For explicit Euler it is order 1. If other conditions
SMST C: Numerical Methods 1–8

are satisfied (and they are for Euler applied to this problem), the order of the LTE gives the convergence
rate observed in the experiment.
Watch out! Different definitions of the LTE are used in various reference books. Sometimes this involves
an extra power of ∆t.

1.3 Runge-Kutta (RK) methods

In this and the next section we introduce the two families of approximations schemes that are used
most often. This section introduces Runge-Kutta (RK) methods and the next section introduces Linear
Multistep Methods (LMMs). After that, we examine their properties in detail in Section 1.5 where we
see that they can give much more accurate and efficient results than explicit Euler.
Runge-Kutta methods are one-step methods. Given the ODE, its approximate solution un ≈ u(tn ) and
the timestep size ∆t, RK nethods give the approximate solution un+1 through a sequence of intermediate
calculations called stages. The approximate solution at time levels before tn is not used in the formula.
Notation varies for RK schemes, but the version below is commonly used.

1.3.1 Explicit RK methods

Perhaps the best-known RK method is the so-called classical 4th order RK method. Sometimes it
is simply (but wrongly) called “the RK method”. It is:
k1 = ∆t f (tn , un )
 
2 1 1 1
k = ∆t f tn + ∆t, un + k ,
2 2
 
1 1 2
k3 = ∆t f tn + ∆t, un + k , (1.11)
2 2
k4 ∆t f tn + ∆t, un + k 3


=
 
1 1 1 2 1 3 1 4
un+1 n
= u + k + k + k + k .
6 3 3 6
Here, armed with the function f (t, u) for our ODE problem, we use ∆t and un to compute k1 , then ∆t,
un and k1 to compute k2 , and so on building sequentially until we have all the parts to compute the
approximate solution at the new time level un+1 . There are 4 intermediate stages required to complete
one step of this scheme. The ki values are simply intermediate calculations which are discarded at the
end of each step and computed afresh on the next step. This RK4 scheme is explicit since we do not
need to solve any algebraic equations to compute the intermediate stages or the final part.
The general s stage RK family of explicit methods can be written as
k1 = ∆t f (tn , un )
2
∆t f tn + c2 ∆t, un + a2,1 k1 ,


k =
k3 ∆t f tn + c3 ∆t, un + a3,1 k1 + a3,2 k2 ,


= (1.12)
... ...
s
∆t f tn + cs ∆t, un + as,1 k1 + · · · + as,s−1 ks−1


k =
Xs
un+1 = un + bi ki .
i=1

The key thing is that every stage of this process involves only an explicit evaluation of a formula involving
previously computed information — un gives k1 gives k2 gives · · · gives ks , and finally all are combined to
give un+1 . These schemes are often written in condensed form in terms of their coefficients in a Butcher
tableau
0 0 0 ··· 0
c2 a21 0 · · · 0
c A . .. .. .. .
T
:= .. . . . ..
b c a a ··· 0
s s1 s2

b1 b2 ··· bs
SMST C: Numerical Methods 1–9

Ps
where ci = j=1 aij (so c1 = 0) and the diagonal and upper portion of coefficient matrix A is zero. The
zero entries in A have not been left as gaps here, despite the common convention in textbooks.

Example 1.4 The Butcher tableau for the Classical RK4 scheme (1.11) is

0 0 0 0 0
1/2 1/2 0 0 0
1/2 0 1/2 0 0 .
1 0 0 1 0
1/6 1/3 1/3 1/6

The general two-stage second order explicit Runge-Kutta is given by the Butcher table

0 0 0
1/(2α) 1/(2α) 0 , (1.13)
1−α α

where α is a parameter. With α = 1/2 we recover the Improved Euler method first seen as (1.10)

k1 = ∆t f (tn , un )
k2 ∆t f tn + ∆t, un + k1 ,


=
1 1
un+1 un + k + k2 .


= (1.14)
2

1.3.2 Implicit RK methods

There are more general s stage implicit RK schemes defined in the same notation by
 
Xs
ki = ∆t f tn + ci ∆t, un + ai,j ki  , i = 1, 2, . . . , s.
j=1
s
X
un+1 = un + bi k i (1.15)
i=1

and in Butcher tableau form as

c1 a11 a12 ··· a1s
c2 a21 a22 ··· a2s
c A . .. .. .. ..
:= .. . . . . (1.16)
bT cs as1 as2 ··· ass
b1 b2 ··· bs
Ps
where ci = j=1 aij . The numbers ai,s , bi , ci may be chosen from a list of schemes in books and tables.
A horrendous coupled system of nonlinear algebraic equations must be solved to find the ki , and so these
fully-implicit RK schemes are reserved for special problems.

Example 1.5 The Implicit Euler scheme (1.8) is an implicit, 1-stage RK method:

k1 = ∆t f tn + ∆t, un + k1

un+1 = un + k1

with Butcher table

1 1
.
1
We have to solve a system of algebraic equations to get k1 — hence the scheme is implicit.

There are lots of different Runge-Kutta schemes available apart from those we have already met above.
SMST C: Numerical Methods 1–10

1.4 Linear Multistep Methods — LMMs

The second family we schemes that is widely used are LMMs. Here values un and f (tn , un ) from more
than one previously computed step are used to form the approximation. To see how such schemes can
arise, consider the following example.

Example 1.6 (Derivation of a multistep method) Integrate the differential equation u′ = f (t, u)
from tn to tn+2 to get
Z tn+2 Z tn+2 Z tn+2
u′ (s) ds = f (s, u(s)) ds ⇔ u(tn+2 ) − u(tn ) = f (s, u(s)) ds.
tn tn tn

Approximate the integral by Simpson’s rule

Z b
b−a
g(x)dx ≈ (g(a) + 4g(c) + g(b))
a 6
where midpoint c = (a + b)/2, to get
∆t
u(tn+2 ) − u(tn ) ≈ f (tn , u(tn )) + 4f (tn+1 , u(tn+1 )) + f (tn+2 , u(tn+2 )) .


3
This leads to the method
∆t
un+2 = un + f (tn , un ) + 4f (tn+1 , un+1 ) + f (tn+2 , un+2 ) .


(1.17)
3
In the example we see that in contrast to the one-step RK methods of Section 1.3, we require un and
un+1 to get un+2 and so (1.17) is a two-step method rather than a one-step method.

1.4.1 General form of LMMs

These schemes are called “linear” because they involve linear combinations of u’s and f ’s, and “multistep”
because (usually) more than one step is involved.
Given a sequence of equally spaced time levels tn = t0 + n∆t with step size ∆t, the general k-step
LMM can be written as
k
X k
X
n+j
αj u = ∆t βj f n+j , |α0 | + |β0 | =
̸ 0, αk ̸= 0 (1.18)
j=0 j=0

and we use the shorthand

f n ≡ f (tn , un ).
The method is defined through the parameters αj , βj , j = 0, 1, 2, . . . , k. Given the approximate solution
up to time tn+k−1 , we obtain the approximate solution un+k at the new time level tn+k from (1.18) as
k−1
X k−1
X
g(un+k ) ≡ αk un+k − ∆t βk f (tn+k , un+k ) + αj un+j − ∆t βj f n+j =0
j=0 j=0
| {z }
all previously computed

written in the form for the Newton iteration (1.2). If βk = 0, then the scheme is explicit since un+k can
be evaluated directly without the need to solve an algebraic system of equations. However, if βk ̸= 0,
then we do have to solve an algebraic system of equations at each time step for un+k and the scheme is
implicit.

The first k steps

Note that to get started, the k-step LMM needs the first k time levels of the approximate solution,
u0 , u1 , . . . , uk−1 , to be specified. The ODE initial conditions only give u0 , so something extra has to
be done. Standard approaches include using a one-step method to get u1 , . . . , uk−1 or using a one-step
method to get u1 , then a two-step method for u2 ,. . . , then a (k − 1)-step method to get uk−1 and then
continue with the k-step method.
SMST C: Numerical Methods 1–11

1.4.2 Examples of LMM Methods

Adams-Bashforth and Adams-Moulton methods have the general form
k
X
un+k = un+k−1 + ∆t βj f n+k
j=0

where given k, the βj are chosen to make the schemes as accurate as possible (Section 1.5.3.) Bashforth
has βk = 0 and Moulton has βk ̸= 0. A typical Adams-Bashforth method (AB4) is
∆t
un+4 = un+3 + 55f n+3 − 59f n+2 + 37f n+1 − 9f n


24
and it is explicit since the approximate solution can be found at each step without solving any algebraic
equations. In contrast, Adams-Moulton 4 is
∆t
un+3 = un+2 + 9f n+3 + 19f n+2 − 5f n+1 − 9f n


24
and is implicit since we do need to solve an algebraic equation to get un+3 .
Another common class of methods are the Backward Difference Formulae (BDF) with general form
k
X 1
∇i∆t un+k = ∆t f n+k (1.19)
i=1
i

where the backward difference operator ∇∆t and its iterates are defined by

∇∆t ur = ur − ur−1 ,
∇2∆t ur = ∇∆t (∇∆t ur )
= ∇∆t ur − ∇∆t ur−1
= (un − un−1 ) − (un−1 − un−2 )

and so on. BDFs are implicit methods. A typical example, written in standard LMM format (1.18)
normalised so that α2 = 1, is BDF2
4 1 2∆t n+2
un+2 − un+1 + un = f . (1.20)
3 3 3
We discuss Backward Difference Formulae some more in Section 2.2, since they have very useful stability
properties.
Schemes with special features are used for special types of problems. For example, the Leap-Frog method
(also known as the explicit midpoint rule) is a two step method

un+2 = un + 2∆tf n+1 .

It has particularly useful properties for Hamiltonian systems, like N -body gravitational problems (such
as the motion of the solar system or star clusters) and molecular dynamics. It is one of the most used
numerical methods, although better known in applications areas under names like the Verlet method.
For completeness, let us note that some, but not all, of the early examples we met can be cast as one-step
LMMs. e.g. Euler and Backward Euler are simple examples of one-step LMMs as well as of RK methods.

1.5 Desired and essential properties of schemes

Preliminary notation

Definition 1.1 For the general LMM (1.18), we define the first and second characteristic polyno-
mials
k
X Xk
ρ(z) = αj z j and σ(z) = βj z j .
j=0 j=0
SMST C: Numerical Methods 1–12

1.5.1 Zero-stability
A starting point for establishing if a numerical method for approximating ODEs is any good or not is to
see if it can solve u′ = 0. The solution of this ODE is clearly

u(t) = u(0) ∀ t ≥ 0.

This sounds too trivial a test problem, but it can actually give vital information about the scheme. For
example the Explicit Euler method (1.3) applied to u′ = 0 gives

un+1 = un n = 0, 1, . . . .

So does the Implicit Euler method (1.8). In fact this is the case for all RK methods.
The property related to solving u′ = 0 that is required for k-step LMMs is actually less demanding than
getting the right answer. It is called Zero Stability.

Theorem 1.1 (Zero Stability and the Root Condition) A linear multi-step method is zero-stable
if and only if all the roots of the characteristic polynomial ρ(z) (Definition 1.1) satisfy |z| ≤ 1, and any
root with |z| = 1 is simple.

The characteristic polynomial ρ(z) is obtained by applying the general LMM (1.18) to u′ = 0 to get
k
X
αj un+j = 0.
j=0

Substituting the ansatz

uj = z j
(z to the power j, not just a superscript label) and dividing by z n gives ρ(z) on the left hand side.

1.5.2 Absolute stability

The next level of test problem is to see if the scheme produces a sensible approximation of

u′ = λu λ ∈ C−

i.e. Re(λ) < 0. This ODE has solution

u(t) = exp(λt)u(0)

and, since Re(λ) < 0,

|u(t)| → 0 as t → ∞.
The approximation scheme is said to be absolutely stable if the analogous result holds for the approx-
imation:
|un | → 0 as n → ∞.
This property often depends on the choice of ∆t in the approximation scheme, and, almost always, it is
actually determined by the value of λ∆t.

Definition 1.2 The range of values of ∆tλ in the complex plane for which the scheme is absolutely
stable is called the Region of Absolute Stability while the section of the negative real axis starting
from hλ = 0 and intersecting the RAS is the Interval of Absolute Stability.

The direct approach often works.

Example 1.7 The explicit Euler scheme applied to u′ = λu (for λ ∈ C) gives

un+1 = un + ∆tλun = (1 + ∆tλ)un

for n = 0, 1, . . .. So
un = (1 + ∆tλ)n u0
SMST C: Numerical Methods 1–13

and un → 0 in general, if and only if |1 + ∆tλ| < 1. So the RAS is the interior of the disk centre −1
radius 1 in the complex ∆tλ plane. Restricting λ to be real gives the IAS to be ∆tλ ∈ (0, −2) where
λ ∈ R. So, when λ < 0 we have to choose
2
0 < ∆t <
−λ
for absolute stability.

Example 1.8 The explicit RK2 scheme (1.14) applied to u′ = λu gives

k1 = ∆tλun
2
∆tλ(un + k 1 ) = ∆tλ + (∆tλ)2 un


k =
 
1 1 1
un+1 n
k + k = 1 + ∆tλ + (∆tλ) un = R(∆tλ)un
2 2


= u +
2 2
for n = 0, 1, . . ., so that
un = R(∆tλ)n u0
and un → 0 in general, if and only if |R(∆tλ)| < 1. So the RAS is the interior of a complicated region
in the complex ∆tλ plane. Restricting λ to be real again gives the IAS to be ∆tλ ∈ (0, −2) where λ ∈ R.
So, when λ < 0 we have to choose
2
0 < ∆t <
−λ
for absolute stability.

In general, for RK methods applied to the test problem u′ = λu we get

un+1 = R(λ∆t)un

where we can get the stability function R directly as above, or by using the Butcher table notation
(1.16)
R(z) = 1 + zbT (I − zA)−1 1,
and 1 := (1, 1, . . . , 1)T . It is relatively easy to show that R(z) is a polynomial in z for explicit RK
methods, and that it is a rational polynomial function for implicit RK methods. We have
n
un = (R(λ∆t)) u0

and clearly we need |R(λ∆t)| < 1 to ensure that un → 0 as n → ∞ for absolute stability.
For k-step LMMs the requirement for absolute stability is that the roots z = zi , i = 1, . . . , k, of the
second characteristic equation
ρ(z) − ∆tλ σ(z) = 0 (1.21)
all satisfy |zi | < 1.

Example 1.9 The explicit Euler scheme has

ρ(z) = z − 1, σ(z) = 1

and (1.21) has root z = 1 + ∆tλ, which is less than 1 in modulus for ∆tλ inside the disk described in
Example 1.7.

1.5.3 Local Truncation Error (LTE)

Warning: There are two main ways to define this quantity given in different tetxbooks. They are directly
related, but differ by a factor of ∆t. I choose to use this one because it matches similar definitions for
PDEs.
We define the Local Truncation Error (LTE) for RK methods in terms of (1.12) to be
Ps
n u(tn+1 ) − u(tn ) − i=1 bi k i
LTE(t ) =
∆t
SMST C: Numerical Methods 1–14

where u(t) is an exact solution of the ODE u′ = f (t, u) and we generate the k i starting with un = u(tn ),
the exact solution.
For general LMMs (1.18) it is defined similarly as
 
k k
1 X X
LTE(tn ) =  αj u(tn+j ) − ∆t βj f (tn+j , u(tn+j ))
αk ∆t j=0 j=0

where again u(t) is an exact solution of the ODE u′ = f (t, u). This is neater if the standard normalisation
αk = 1 has been used to specify the scheme.
What we want to know is the behaviour of the LTE as ∆t → 0 since this tells us something about the
convergence properties of the scheme. The order of the LTE is the highest number p satisfying

LTE = O(∆tp ).

Very detailed multivariable Taylor expansion is needed for RK methods, and the following results are
available. For s stage implicit RK, methods of order p = 2s have been derived. However, for s stage
explicit RK, the maximum order is at best p = s, and often p < s as shown below.

Explicit RK stages s= 1 2 3 4 5 6 7 8 9 10
p
Best LTE is O(∆t ) with p= 1 2 3 4 4 5 6 6 7 7

For LMMs the Taylor expansion is much more straightforward. We use the ODE to replace the f by u′
(written as u(1) below) and get
k
X k
X
(αk ∆t) LTE(tn ) = αj u(tn+j ) − ∆t βj f (tn+j , u(tn+j ))
j=0 j=0
k
X k
X
= αj u(tn + j∆t) − ∆t βj u(1) (tn + j∆t)
j=0 j=0
k
(j∆t)2 (2) n
X  
= αj u(tn ) + (j∆t)u(1) (tn ) + u (t ) + · · ·
j=0
2!
k
(j∆t)2 (3) n
X  
−∆t βj u(1) (tn ) + (j∆t)u(2) (tn ) + u (t ) + · · · .
j=0
2!

Careful bookkeeping (or the use of Maple or Mathematica) gives the required results, and we get

LTE(tn ) = Cp+1 ∆tp u(p+1) (tn ) + O(∆tp+1 ) = O(∆tp )

where Cp+1 is simply a number that depends on the details of the scheme. If you really need to know
what it is (and sometimes that is useful), it is given by
 
k k
1  1 X 1 X
Cp+1 = j p+1 αj − j p βj  .
αk (p + 1)! j=0 p! j=0

An alternative and equivalent way to establish the order of accuracy for LMMs is given by the following.
It uses the characteristic polynomials defined in Definition 1.1.

Theorem 1.2 (eg [6]) The multistep method is of order p ≥ 1 if and only if there exists c ̸= 0 such that

ρ(z) − σ(z) ln(z) = c(z − 1)p+1 + O(|z − 1|)p+2 .

Once again, careful bookkeeping of Taylor expansions (this time expanding about z = 1) is required.
SMST C: Numerical Methods 1–15

1.5.4 Consistency of scheme and initial data

The approximation scheme is consistent with the ODE if

LTE → 0 as ∆t → 0.

So the order of the LTE has to be at least 1.

RK methods (1.12) and (1.15) are consistent if
s
X
bi = 1.
i=1

LMM (1.18) is consistent if both

k
X k
X k
X
αj = 0 and jαj − βj = 0,
j=0 j=0 j=0

or equivalently if both ρ(1) = 0 and ρ′ (1) − σ(1) = 0.

1.5.5 Consistency of initial data

As noted in Section 1.4.1, a k-step LMM needs the first k time levels of the approximate solution
u0 , u1 , . . . , uk−1 to get started. These data are said to be consistent if

us → u(0) as ∆t → 0

for each s = 0, . . . , k − 1.
Initial data for RK schemes are consistent if they satisfy the above with k = 1. Of course this is trivial if
one does the obvious thing and uses u0 = u(0), the initial data supplied with the original ODE problem.

1.5.6 Convergence
This is the real goal. What we want is that the error in the approximate solution everywhere over the
fixed time interval t ∈ [0, T ] should tend to zero as the step size ∆t → 0. One could write this for a scalar
problem as
||uapprox (t) − u(t)||L∞ [0,T ] → 0 as ∆t → 0.
Fortunately for most well-behaved ODE systems, LMMs satisfy the Dahlquist Equivalence Theorem.

Zero Stability + Consistency of scheme and initial data ⇔ Convergence.

If the initial data are in addition accurate enough, then the rate of convergence is the order of the LTE.
See e.g. [1, Thm. 5.1] and [8, Thm. 12.5] for some technical details.
So checking two relatively “simple” properties of the approximation scheme is enough to establish con-
vergence, and it does not have to be tackled directly.
Similar results hold for RK methods (and one-step methods in general). See e.g. [1, Thm. 3.1] or the
comprehensive coverage in [4].

References
[1] U. Ascher and L. Petzold. Computer Methods for Ordinary Differential Equations and Differential
Algebraic Equations. SIAM, 1998.
[2] K. E. Atkinson. An Introduction to Numerical Analysis. Wiley, 2 edition, 1989.

[3] R. L. Burden and J. D. Faires. Numerical Analysis. Brooks Cole, 7 edition, 2001.
[4] J. C. Butcher. The Numerical Analysis of ODEs. Wiley-Interscience, 1987.
SMST C: Numerical Methods 1–16

[5] D. F. Griffiths and D. J. Higham. Numerical Methods for Ordinary Differential Equations. Springer,
2010.
[6] A. Iserles. A first course in the Numerical Analysis of Differential Equations. Cambridge texts in
applied mathematics. CUP, Cambridge, UK, 1996.
[7] A. M. Stuart and A. R. Humphries. Dynamical Systems and Numerical Analysis. Cambridge Mono-
graphs on Applied and Computational Mathematics. CUP, Cambridge, UK, 1998.
[8] E. Süli and D. Mayers. An Introduction to Numerical Analysis. CUP, 2003.
 SMSTC (2024/25)
Numerical Methods
Lecture 2: Numerical Methods for ODEs 2
Dugald B Duncan, Heriot-Watt Universitya

www.smstc.ac.uk

Contents
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2 Systems of ODEs and stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2.1 Linear systems of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–2
2.2.2 Nonlinear systems of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
2.2.3 Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–3
2.2.4 A-stability and A(α)-stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–4
2.3 Computational effort and efficiency . . . . . . . . . . . . . . . . . . . . . . . 2–5
2.4 Practical global error estimation for fixed step schemes . . . . . . . . . . . 2–5
2.5 Error Estimation and Time Step Adaptivity . . . . . . . . . . . . . . . . . . 2–7
2.5.1 Local error estimation for RK methods . . . . . . . . . . . . . . . . . . . . . . 2–8
2.5.2 Local error estimation for implicit LMMs . . . . . . . . . . . . . . . . . . . . . 2–8
2.5.3 Time step control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–9
2.5.4 How is the global error related to the local error tolerance? . . . . . . . . . . . 2–10
2.6 What solvers are available and how do we use them? . . . . . . . . . . . . . 2–10
2.6.1 Explicit solvers for non-stiff problems . . . . . . . . . . . . . . . . . . . . . . . 2–11
2.6.2 Implicit LMM solvers for stiff problems . . . . . . . . . . . . . . . . . . . . . . 2–11
2.6.3 Other sources of ODE software . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–13
2.7 The Method of Lines for simple, time-dependent PDEs . . . . . . . . . . . 2–16
2.7.1 Output and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–18
2.8 What has not been covered? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2–19

a D.B.Duncan@hw.ac.uk

2–1
SMST C: Numerical Methods 2–2

2.1 Introduction
These two lectures are about the accurate and efficient approximation of the solution of systems of
ordinary differential equations - hereafter abbreviated to “ODEs”. The aim is to introduce key ideas and
the language used to describe the attributes of numerical methods, as opposed to doing rigorous analysis.
There is usually no need to write your own solver for ODEs. There are lots of good packages available
that many people have devoted years of their life to. However it is useful to have a good idea of the
theory and ideas that underlie solving ODEs and hence understand the limitations and advantages of
packages.
This lecture aims to give an overview of practical issues related to using methods and software packages.

2.2 Systems of ODEs and stiffness

All of the schemes we have considered apply to systems of ODEs, and much of the analysis of their
properties does as well. However the absolute stability analysis used so far only deals with the scalar
problem u′ = λu, and we need to go into more detail to reveal an important aspect of the numerical
solution of ODEs called stiffness.

2.2.1 Linear systems of ODEs

Consider the M × M linear system of ODEs

u′ = Au (2.1)

where for simplicity we assume that matrix A ∈ RM ×M can be written as

A = SDS −1

where S, D are matrices of the same dimension as A, and D is diagonal. If such S, D exist, then A and
D are said to be similar matrices and they have the same eigenvalues λµ (A), µ = 1, . . . , M . Hence the
diagonal elements of D are the eigenvalues of A.
Now define v(t) = S −1 u(t) for all t ≥ 0. It can be shown that applying a standard RK or LMM to (2.1)
is equivalent to applying the same scheme to the uncoupled system of ODEs

v′ = Dv, (2.2)

in the sense that the approximations generated satisfy vn = S −1 un , for all n ≥ 0. The scalar component
equations of (2.2) are
vµ′ = λµ (A)vµ , µ = 1, . . . , M (2.3)
where λµ (A) are the eigenvalues of A. Note that the eigenvalues are either real or in complex conjugate
pairs since A is real.
So the overall absolute stability of an approximation of the linear system (2.1) is determined by the
absolute stability of it applied to the eigencomponent ODEs (2.3). For the full system to be stable, all
of the eigencomponent problems must be stable and hence absolute stability is determined by the worst
cases in (2.3).

Example 2.1 Consider the linear ODE

u′ = Au = SDS −1 u

where    
−1 0 1 1
D= , S= .
0 −100 −0.8 0.2
The solution of the ODE is:
   
1 1
u(t) = c1 exp(−t) + c2 exp(−100t) (2.4)
−0.8 0.2
| {z } | {z }
slow decay fast decay
SMST C: Numerical Methods 2–3

1.5

U0 and U1
0.5

-0.5

-1
0 0.5 1 1.5 2 2.5 3 3.5 4
time t

Figure 2.1: Solution (u1 , u2 ) given by (2.4). Over a very short timescale we see decay from the fast
eigenvalue, then decay is dominated by the faster eigenvalue. For an explicit fixed step scheme the step
size is determined by the fast time scale. The fast transient at the start is like a boundary layer from the
earlier part of the course.

where c1 , c2 must be chosen to satisfy the initial condition u(0) = u0 . The solution is shown in Figure 2.1.
Stability for explicit Euler is restricted by the “fast” or big eigenvalue λ2 = −100 and we have
100∆t < 2 → ∆t < 2/100 = 0.02,
but the actual solution is completely dominated by the slow component after a short initial period. The
slow component would be reasonably well resolved with time steps
∆t = O(1/|λ1 |) = O(1)
but stability requires something much smaller. So the acceptable step size ∆t and overall effort is de-
termined by something that only contributes to the solution for a short time. After an initial period,
accuracy is determined only by the slow component.
Stiffness is showing up here in the large ratio of the timescales in the problem.

2.2.2 Nonlinear systems of ODEs

The general nonlinear problem u′ = f (t, u) is much harder to analyse for stability, but some information
can be gained by linearising it about a representative solution value (say about u∗ ) and apply the linear
system analysis to the associated problem
 
∂f
w′ = Aw where A = (t, u∗ ) .
∂u
Here A is the Jacobian matrix of the vector function f with respect to its argument u, evaluated at (t, u∗ ).

2.2.3 Stiffness
It is hard to get a clean definition of stiffness that covers all cases, because it can arise in different ways.

ˆ Stiff problems have timescales in the solution which differ by a large factor.
ˆ For the linear ODE u′ = Au where A has eigenvalues λj < 0, the traditional definition of stiffness
is that
def maxj |λj |
Stiffness Ratio = >> 1 .
minj |λj |
SMST C: Numerical Methods 2–4

Large ratio of timescales.

ˆ Higham and Trefethen:

It is generally agreed that the essence of stiffness is a simple idea: stability is more of a constraint
than accuracy.
ˆ Hairer and Wanner:
Stiff equations are problems for which explicit methods don’t work.
ˆ The step-size for an explicit scheme needs to be chosen to maintain stability rather than for the
accuracy of the scheme.

Example of extreme stiffness.

The Becker-Döring equation is an innocuous looking system of ODEs

u′ = A(u1 )u ,

where the matrix A is essentially tridiagonal and depends only on the first component of u. Typical
solutions evolve on time scales of O(1020 ) and the eigenvalues of the linearised problem are real and
negative with maxj |λj | ≈ 10.
The Euler stability limit is then ∆t < 0.2 implying that O(1021 ) steps would be required to compute
the solution. Unfortunately the CPU time per step ≈ 0.001 seconds on a powerful workstation, giving a
total CPU time of 1018 seconds or O(1010 ) years. One could wait until machines are 106 faster (in about
20 years) and do it on a 1000 processor parallel computer, but the calculation using explicit Euler would
still take about 10 years!
A good implicit package based on BDF methods takes only a few minutes.

2.2.4 A-stability and A(α)-stability

Schemes with no stability restrictions would be ideal for stiff problems. Following the definition of
absolute stability in Section 1.5.2, a scheme is said to be A-stable if the region of absolute stability
(RAS) contains the whole of the left half complex plane. i.e. all of Re(∆tλ) < 0.
Notes:

ˆ The simplest A-stable methods are the implicit Euler method (1.8) (with R(z) = 1/(1 − z)) and
the Trapezoidal Rule (1.9) (with R(z) = (1 + z/2)/(1 − z/2)).
ˆ Explicit RK methods are not A-stable. This is because R(z) is a polynomial in z for an explicit
method, and polynomials of degree 1 or more are unbounded as |z| → ∞.

ˆ Explicit LMMs are not A-stable.

Even the BDF methods (1.19) which are designed for stiff problems cannot combine high order LTEs,
zero stability and A-stability, but they do satisfy a slightly less restrictive stability requirement that their
RAS contains the sector in the left half complex plane with angle [180◦ − α◦ , 180◦ + α◦ ]. This is called
A(α)-stability. When α > 0 this includes all of the negative real axis, which is useful for stiff problems
with pure real eigenvalues. Note that only BDF1—BDF6 are zero-stable and hence convergent/useful.

BDF # 1 2 3 4 5 6 ≥7
√ √ √ √ √ √
Zero-stable No
LTE order 1 2 3 4 5 6 N.A.
◦ ◦ ◦ ◦ ◦ ◦ ◦
α 90 90 88 73 51 18 N.A.
SMST C: Numerical Methods 2–5

2.3 Computational effort and efficiency

For systems of ODEs the computational cost is usually measured in terms of the computer time (CPU
time) needed to do the calculation, and the amount of computer memory used. These both depend on
the details of the problem and on the scheme used, but some comparisons can be made by considering
the contributing factors to both of these. We assume below that we are dealing with an ODE system
where the evaluation of f (t, u) involves significant computational cost compared to e.g. simply adding
two vectors of the same size together.
For explicit methods, the CPU time is roughly proportional to the number of times the function f (t, u)
is evaluated. For RK methods that is the number of steps multiplied by the number of stages in each step.
For LMMs, with careful storage of previously computed f values, we only need one new f evaluation per
step, so it simply depends on the number of steps. The storage required is proportional to the number
of stages s for explicit RK and to k for k-step explicit LMMs.
For implicit methods, the CPU time depends on the contributions above with significant additional
cost for solving the algebraic equations on each step — e.g. by Newton (1.2). The storage costs also
increase since the Jacobian matrix must be computed, stored and space is required for solution of the
linear systems.
To cut CPU time, one can use higher accuracy or more stable implicit schemes with bigger steps so
long as the higher cost per step is outweighed by their number being reduced. Also, we should consider
adapting the time steps size to be small where errors are expected to be large and large when errors are
small. That way we concentrate work where it is required rather than uniformly over the whole time
interval.
We consider this and other ways of dealing with efficiency further in the next lecture.

2.4 Practical global error estimation for fixed step schemes

Summary: Compare the approximate solutions obtained with two different step sizes and use the differ-
ence between them as an estimate of the global error.
We need some new notation to make this discussion easier to deal with.

Definition 2.1 We rewrite the approximate solution obtained with fixed step size ∆t to show the depen-
dence on ∆t as u∆t (tn ) ≡ un , where tn = t0 + n∆t for n = 0, 1, 2, . . ..

If we have a scheme with LTE = O(∆tp ) for some integer p ≥ 1, and the ODE it is applied to is
well-enough behaved, then the convergence result in Section 1.5.6 tells us that the Global Error
u(t) − u∆t (t) = O(∆tp ) as ∆t → 0
for all t = tn ∈ [0, T ], the finite time interval of interest. Another way to write this when ∆t is small
enough is
u∆t (t) = u(t) + O(∆tp ) ≈ u(t) + c(t)∆tp
where the exact solution u and c do not depend on ∆t. If we compare approximate solutions obtained
with different stepsizes, usually ∆t and ∆t/2, then
∆tp
 
p
u∆t (t) − u∆t/2 (t) ≈ (u(t) + c(t) ∆t ) − u(t) + c(t) p
2
 
1
= c(t) ∆tp 1 − p
2
 
1
≈ 1 − p (u(t) − u∆t (t)) .
2
So the difference between approximate solutions with two different step sizes is approximately proportional
to the Global Error:  p 
2
∥u(t) − u∆t (t)∥ ≈ ∥u∆t (t) − u∆t/2 (t)∥ (2.5)
2p − 1
using an appropriate norm for vector-valued solutions. Table 2.1 gives an example of this for the logistic
ODE.
SMST C: Numerical Methods 2–6

Relative Error (measured) vs CPU time (measured)

0
10

−2
10

−4
10

−6
10

−8
10

Euler
Huen
−10
10 RK4
AB4
−4 −3 −2 −1 0
10 10 10 10 10
Computer Time

Figure 2.2: Maximum relative error in the approximate solution of the logistic equation u′ = u(1 − u)
with u(0) = 0.02 for t ∈ [0, 10] by various explicit methods plotted against the CPU time required on
my old laptop. The step size decreases from ∆t = 10/4 in multiples of 2 to ∆t = 10/215 — result
shown by a marker symbol for each scheme. The most efficient scheme is the one that gives the required
error (say 10−4 ) with the least CPU time. For all errors less than about 0.001 in this example, it is
Adams-Bashforth 4 (AB4). We see AB4 and the Classical RK method (RK4) decreasing in parallel since
both are 4th order accurate methods. Huen decreases with about 1/2 the slope since it is 2nd order, and
Euler decreases with 1/2 the slope again since it is first order.
SMST C: Numerical Methods 2–7

Number Step Size Approximate 4

of steps N ∆t = 5/N u(5) estimated error
1 5 0.32919 2 actual error
2 2.5 0.560787747852

log10(error)
4 1.25 0.670943530629 0
8 0.625 0.725665614442
16 0.3125 0.744405428442 −2
32 0.15625 0.749814551585
64 0.078125 0.751277482656
−4
128 0.039062 0.751659099945
256 0.019531 0.751756651338
512 0.0097656 0.751781318874 −6
−2 −1 0 1
1024 0.0048828 0.751787521467 log10(∆ t)

Table 2.1: The approximate solution at t = 5 for the well-behaved ODE u′ = u(1 − u) with u(0) = 0.02
using the Improved Euler method (1.14) with fixed step size ∆t = 5/N . The Figure shows the error
estimate (4/3)|u∆t (5) − u∆t/2 (5)| obtained from (2.5) with p = 2 as well as the actual error. There is
little difference except when the step size is large and the estimation process that gives (2.5) is inaccurate.

The slope of the graph in Table 2.1 asymptotes to 2 as ∆t → 0 because the scheme (1.14) is 2nd order
accurate. This sort of graph can act as a check that a fixed step size scheme is working correctly. It can
do this by showing that it converges with the correct rate (or not), by providing a global error estimate
and by indicating where the step sizes are small enough to enter the convergence zone so that we can
trust the error estimate and hence the approximate solution.

2.5 Error Estimation and Time Step Adaptivity

Summary: Compare the approximate solutions obtained over each step using two different schemes, use
the difference between the approximations as an estimate of the local error, and use that to select the time
step size to use locally.
So far we have assumed that the time step size will be fixed throughout a calculation. However it is often
much more efficient to adapt the time step size as the solution evolves in time to make it small when
errors are expected to be large and big when errors are small. That way the work is concentrated where
required, and we don’t waste effort where it is not. We need notation for variable timestep sizes.

Definition 2.2 Given the time levels t0 , t1 , . . ., not necessarily equally spaced, the local time step size is
defined by
∆tn−1/2 = tn − tn−1
for n = 1, 2, . . .. Looking at this the other way, given the sequence of timestep sizes ∆t1/2 , ∆t3/2 , . . ., the
time levels tn are defined by
n
X
tn = tn−1 + ∆tn−1/2 = t0 + ∆tj−1/2 .
j=1

Variable time step sizes are used in virtually all modern codes for ODEs. The usual approach in ODE
codes is to use an estimate of the local error (LE) defined next to determine what step size to use.

Definition 2.3 The local error is defined to be the error over the step from t = tn to t = tn+1 assuming
that no errors have been made up to time t = tn . In other words, apply the scheme to the local exact
problem û′ = f (t, û) with initial condition û(tn ) = un and calculate the error in solving this local problem
over one step.

A consequence of this definition is that there is a simple relationship between local error (LE) and the
local truncation error (LTE) defined in Section 1.5.3:

LE = ∆t LTE,
SMST C: Numerical Methods 2–8

and this helps when dealing with error estimates.

We start with ways to estimate the local error in RK methods and LMMs.

2.5.1 Local error estimation for RK methods

Assume that we have reached time t = tn , we have the approximation un ≈ u(tn ) and we are to use time
step size ∆t for the next step. We compute approximate solutions at tn+1 = tn + ∆t by two different
schemes, giving results un+1 by our main scheme and un+1
aux by an auxiliary scheme. If the main scheme
has LTE = O(∆tp ) and the auxiliary scheme has LTEaux = O(∆tq ) with q > p, then
..
un+1 = û(tn+1 ) + LE . un+1
aux = û(tn+1 ) + LEaux
..
= û(tn+1 ) + ∆t LTE . = û(tn+1 ) + ∆t LTEaux
..
= û(tn+1 ) + O(∆tp+1 ) . = û(tn+1 ) + O(∆tq+1 )
where the exact local solution is given by û′ = f (t, û) with û(tn ) = un . So there exists a constant C
which depends on tn but not ∆t such that
∥LE∥ ≈ C ∆tp+1
and
un+1 − un+1
aux = (û(tn+1 ) − û(tn+1 )) + LE − LEaux
= LE + O(∆tq+1 )
≈ LE
when ∆t is small enough and q > p. Clearly
∥LEest (∆t)∥ = ∥un+1 − un+1
aux ∥ (2.6)
is a practical approximation of local error when ∆t is small enough. The auxiliary solution un+1
aux is used
for this estimate and then discarded.

Embedded pairs

To cut down the computational effort, a pairs of RK methods of order p and p + 1 are constructed sharing
as many of the same stage calculations (the k in (1.12)) as possible. The simplest example is the forward
Euler method embedded in a modified trapezoidal rule:
k1 = ∆t f (tn , un )
k2 = ∆t f (tn + ∆t, un + k1 )
n+1
u = un + k1
1 1
un+1 = un + k + k2


aux
2
LE estimate = un+1 − un+1
aux .

Probably the most famous embedded formula is the explicit Fehlberg 4(5) pair which uses a method of
order 4 with an auxiliary scheme of order 5. It has 6 different stages (the k evaluations) instead of at least
9 for two arbitrarily selected explicit RK methods of these orders. The details are not very instructive but
can be found for example in [1]. The methods of Dormand and Prince are now probably more commonly
used. The 4(5) pair which is the basis of the Matlab code has 7 stages.
In practice embedded RK pairs are usually based on explicit schemes and so are not often used for stiff
problems.

2.5.2 Local error estimation for implicit LMMs

Implicit LMMs are in widespread use to deal with stiff ODEs. Let us rewrite the general implicit LMM
(1.18) in an equivalent form for this Section. We have
k−1
X k−1
X
αk un+1 − βk ∆tn+1/2 f (tn+1 , un+1 ) + αj un+j−k+1 − ∆tn+1/2 βj f n+j−k+1 ≡ g(un+1 ) = 0,
j=0 j=0
SMST C: Numerical Methods 2–9

and we need to solve the nonlinear system of algebraic equations, g(un+1 ) = 0 for un+1 at each timestep.
This is usually done using Newton’s method (Section 1.1.3) or a variant of it. The initial guess for un+1
may be obtained by using an explicit LMM called the predictor in the literature. The implicit method
we actually want to use is called the corrector. If we match the explicit and implicit LMMs to have the
same order of accuracy (say p) we get the following way to estimate the local error.
The predictor method can be chosen to play the part of the auxiliary scheme in the previous subsection
about RK methods. Typically a predictor of the same order of accuracy as the corrector is used. From
the LTE definition of Section (1.5.3), and from the relationship LE = ∆t LTE, we have

LE = Cp+1 ∆tp+1 û(p+1) (tn ) + O(∆tp+2 ), LEaux = C̃p+1 ∆tp+1 û(p+1) (tn ) + O(∆tp+2 )

for the main scheme (the corrector) and auxiliary scheme (the predictor) respectively where the exact
local solution is given by û′ = f (t, û) with û(tn ) = un . Then assuming that Cp+1 ̸= C̃p+1 , we follow the
approach of the previous subsection to get

un+1 − un+1
aux = (û(tn+1 ) − û(tn+1 )) + LE − LEaux
= (Cp+1 − C̃p+1 ) ∆tp+1 û(p+1) (tn ) + O(∆tp+2 )
!
Cp+1 − C̃p+1
= LE + O(∆tp+2 )
Cp+1
!
Cp+1 − C̃p+1
≈ LE
Cp+1

when ∆t is small enough. So

!
Cp+1
LEest = (un+1 − un+1
aux )
Cp+1 − C̃p+1

is an estimate of the local error and there exists constant C, independent of ∆t such that

Cp+1 p+1
∥LEest ∥ = ∥un+1 − un+1
aux ∥ ≈ C ∆t . (2.7)
Cp+1 − C̃p+1

Again the auxiliary solution un+1

aux is used to provide the estimate and then discarded. This estimate is
often called the Milne Device in textbooks.

2.5.3 Time step control

Now we use the local error estimate to control the time step size. Often an efficient strategy is to pick a
tolerance called TOL for the maximum acceptable local error and to choose the local step size ∆tn+1/2
such that
∥LE(∆tn+1/2 )∥ ≤ TOL with ∥LE(∆tn+1/2 )∥ ≈ TOL.
Roughly speaking this attempts to make the local error the same size — uniformly good or uniformly bad
depending on your viewpoint – for each timestep over the whole time interval of the calculation. Most
computer codes define TOL by specifying a relative and absolute tolerance in combination like this,

TOL = reltol∥un+1 ∥ + abstol, (2.8)

or by some similar approach.

Now we derive a formula for the step size that achieves this local error control. If we take a step of size
∆t from t = tn and obtain the local error estimate LEest (∆t), then we can estimate the optimal step size
∆t∗ we should have taken to make ∥LE(∆t∗ )∥ ≈ TOL from

∥LEest (∆t)∥
∥LEest (∆t)∥ ≈ C∆tp+1 ⇔ C≈
∆tp+1
and so
TOL TOL
∥LE(∆t∗ )∥ ≈ C∆tp+1
∗ ≈ TOL ⇔ ∆tp+1
∗ ≈ ≈ ∆tp+1 .
C ∥LEest (∆t)∥
SMST C: Numerical Methods 2–10

Hence the optimal step size is given by

 1/p+1
TOL
∆t∗ ≈ ∆t
∥LEest (∆t)∥

and the local error estimate is given by (2.6) for RK and (2.7) for LMMs. To be on the safe side of all
these estimates, in most computer codes the optimal step size is given by
 1/p+1
TOL
∆t∗ = Fsafe ∆t (2.9)
∥LEest (∆t)∥

where the safety factor Fsafe ∈ (0, 1).

A basic algorithm to use error estimation and time step control, written in the notation of Definition 2.2
is:

2–1. Set n = 0. Guess the initial step size ∆t1/2 .

2–2. Take a step of size ∆tn+1/2 from t = tn to tn+1 .
2–3. Estimate the Local Error (LEest ) in this step.

2–4. Set
!1/p+1
n+1/2 TOL
∆t∗ = Fsafe ∆t .
∥LEest (∆tn+1/2 )∥

2–5. If LEest > TOL, then the local error is too big.
Reset the step size to be ∆tn+1/2 := ∆t∗ .
Goto Algorithm Step 2–2 and try again.
2–6. If LEest ≤ TOL, then the local error is small enough.
Reset n := n + 1.
Set the new step size to be ∆tn+1/2 = ∆t∗ .
Goto Algorithm Step 2–2 and continue.

To make this algorithm work in real life one needs to add a stopping condition so that we stop the
calculation once we have approximated the solution all the way from t = 0 to t = T . Additionally one
might also add a mechanism to choose the step size to give outputs at pre-selected values of t of particular
interest.
Using a one-step scheme with a different size for each step is straightforward. However, you might have
noticed that k-step LMMs are based on having past information over k equal length steps and we lose
that if we have non-uniform step sizes. A workaround is required. We are not going to cover it here, but
is used routinely in ODE codes.

2.5.4 How is the global error related to the local error tolerance?
This was a hot topic in the 1980’s, but seems to be overlooked when discussing ODE solvers these days.
It is important, since specifying reltol and abstol in (2.8) does not give a guaranteed error bound for the
global error. It only directly controls the size local error estimate. The relationship between the local
error tolerance, and the global (or actual) error is not straightforward.
Fortunately the global error is usually more or less proportional to the local error tolerance. See Figure
2.3 for a simple example that works well. There we see that when reltol is small, the global error in each
case is roughly proportional to reltol, and is between 5 and 10 times reltol.

2.6 What solvers are available and how do we use them?

We will concentrate on the Matlab solvers. See the Matlab on-line help system or the book [3] for
full details. Note that Matlab function M-files are essentially the same as Fortran subroutines, and
SMST C: Numerical Methods 2–11

Global Error vs ... Global Error vs ...

0 0
10 10
ode23
ode45

−5 −5
10 10

−10 −10
10 10

−10 −5 0 −2 −1 0 1
10 10 10 10 10 10 10
Relative Tolerance Computer Time

Figure 2.3: Global error in the solution of the logistic ODE u′ = u(1 − u) with u(0) = 0.02 for t ∈ [0, 10]
using the Matlabembedded pair explicit RK codes ode23 and ode45 with abstol = 10−30 . The left plot
shows the relationship between the specified relative tolerance reltol in (2.8) and the global error. The
right hand plot shows the computational cost. The 4th order code ode45 is more efficient than the second
order code ode23.

similar things in Python, C etc.. Matlab has a range of very reliable stiff and non-stiff solvers with
good documentation. They are easily interchangeable and deal automatically with sparse linear algebra.
They can deal with combinations of ODEs and algebraic equations (called DAEs) as well.

2.6.1 Explicit solvers for non-stiff problems

The standard solvers available for non-stiff problems are embedded RK pair methods, like ode45 in
Matlab. The main thing you need to supply (in Matlab) is a function M-file that specifies the function
f (t, u) for the problem under consideration. You also need to specify the initial data, the output times
for the results and (optionally) the relative and absolute tolerances reltol and abstol in the style of (2.8).

Example 2.2 Approximate the solution of the forced oscillator equation

ϕ′′ + (ϕ2 − 4 − sin t)ϕ′ + ϕ3 = 0
with ϕ(0) = 0.5, ϕ′ (0) = −0.5 using a nonstiff solver for t ∈ [0, 10]. Display both ϕ and ϕ′ .
First rewrite the problem in first order form defining u1 (t) = ϕ(t), u2 (t) = ϕ′ (t):
u′1 = u2
u′2 = −(u21 − 4 − sin t)u2 − u31
with u(0) = (u1 (0), u2 (0))T = (0.5, −0.5)T .
Next encode this in a function M-file to define the function f (t, u) as oscode.m in Table 2.2. Finally,
write a script M-file runosc.m to run ode45 and plot the results (see Figure 2.4). This choice of tout in
runosc.m forces output at times 0, 0.1, 0.2, . . . , 10.

2.6.2 Implicit LMM solvers for stiff problems

Matlab uses a variant of the standard LMM backward difference formulae (1.19) in the code ode15s to
deal with stiff ODEs. These schemes are necessarily implicit and so require the solution of a nonlinear
system of algebraic equations at each time step. The full implementation details are complicated and are
described in [5].
We can see the main points by considering the general implicit LMM
k−1
X k−1
X
αk un+1 − βk ∆tn+1/2 f (tn+1 , un+1 ) + αj un+j−k+1 − ∆tn+1/2 βj f n+j−k+1 ≡ g(un+1 ) = 0
j=0 j=0
SMST C: Numerical Methods 2–12

Filename: oscode.m

function dudt = oscode(t,u)

%OSCODE f(t,u) for phi’’ + (phi^2 - 4 - sin t)phi’ + phi^3 = 0

dudt = [u(2);
-(u(1)^2-4-sin(t))*u(2) - u(1)^3];
Filename: runosc.m

%RUNOSC run the oscillator example

tout = 0:0.1:10; %%% output times

u0 = [0.5;-0.5]; %%% initial data

options = odeset(’reltol’,1.e-6,’abstol’,1.e-8); %%% solver parameters

[t,u] = ode45(@oscode,tout,u0,options); %%% run the solver

phi = u(:,1); %%% extract the solution

phiprime = u(:,2); %%% components phi,phi’

subplot(2,1,1)
plot(t,phi,’b-’,t,phiprime,’r--’) %%% plot the results
legend(’\phi’,’d\phi/dt’)
xlabel(’time t’,’fontsize’,14)

print -depsc oscfig %%% save graph as oscfig.eps

Table 2.2: M-files for the oscillator Example 2.2.

30
φ
20 dφ/dt

10

−10

−20
0 1 2 3 4 5 6 7 8 9 10
time t

Figure 2.4: Results of Example 2.2.

SMST C: Numerical Methods 2–13

which is obtained directly from (1.18) by substituting n := n − k + 1. Using Newton’s method 1.1.3 to
solve g(un+1 ) = 0 for un+1 on each time step, we get the iteration
 −1
n+1,[ν+1] n+1,[ν] ∂g n+1,[ν]
u =u − (u ) g(un+1,[ν] ) for ν = 0, 1, 2 . . .
∂u

where the initial guess un+1,[0] is provided by a predictor method and we stop when ∥un+1,[ν+1] −un+1,[ν] ∥
is small enough. Note that the Jacobian matrix is
∂g ∂f
= αk I − βk ∆tn+1/2 ,
∂u ∂u
where I is the M × M identity matrix, all commensurate with the dimension of u′ = f (t, u).
Clearly need to specify the function f (t, u) to make any ODE solver work (otherwise how would it know
which ODE it was to solve?). Additionally for implicit methods we need a way to calculate the Jacobian
matrix, and that is usually best done by providing another M-file or subroutine to calculate all the entries
in ∂f /∂u.
Further, for large systems, the linear algebra required in the Newton iteration can dominate the computer
time and it makes sense to exploit all possible ways to speed this up, like using sparse solvers. Fortunately
Matlab deals with sparse matrices automatically if you make sure that the output from the routine you
supply to calculate the Jacobian produces a matrix of sparse type.

Example 2.3 Approximate the solution of the Robertson chemical reaction system [3, Ch. 12.2.3]

u′1 = −0.04u1 + 104 u2 u3

u′2 = 0.04u1 − 104 u2 u3 − 3 × 107 u22
u′3 = 3 × 107 u22

with u(0) = (u1 (0), u2 (0), u3 (0))T = (1, 0, 0)T using a stiff solver for t ∈ [0, 5].
First encode this in a function M-file to define the function f (t, u) in chemode.m of Table 2.3. Next create
function M-file chemjac.m to calculate the Jacobian matrix. Finally, create a script M-file runchem.m to
run ode15s and plot the results (see Figure 2.5). This choice of tout in runchem.m gives output at the
time levels used by the solver.

2.6.3 Other sources of ODE software

This list is far from comprehensive.

ˆ Python, or more specifically the SciPy library, has the code solve ivp which covers some of the
same methods as Matlab.
ˆ The Nag library has many F77/C callable routines for solving ODEs making use of a wide range
of special structures in the linear algebra.
ˆ The Netlib WWW repository www.netlib.org contains a mountain of software (usually F77
source) for solving all sorts of things. The most appropriate packages for ODEs/DAEs are under
ode and I suggest vode, vodpk for ODEs and DASSL and DASPK are well respected LMM codes
for differential algebraic equations.
ˆ RK Stiff codes publicly available include Radau5 [2] and Stride by Burrage, Butcher and Chipman.

ˆ The Octave programming environment uses lsode and DASSL.

ˆ The SciLab programming environment also uses lsode.

SMST C: Numerical Methods 2–14

Filename: chemode.m

function dudt = chemode(t,u)

%CHEMODE f(t,u) for chemical reaction example

dudt = [-0.04*u(1) + 1.e4*u(2)*u(3);

0.04*u(1) - 1.e4*u(2)*u(3) - 3.e7*u(2)^2;
3.e7*u(2)^2];
Filename: chemjac.m

function dfdu = chemjac(t,u)

%CHEMJAC Jacobian matrix df/du for chemical reaction example

dfdu = [-0.04, 1.e4*u(3) , 1.e4*u(2);

0.04, - 1.e4*u(3) - 6.e7*u(2),-1.e4*u(2);
0 , 6.e7*u(2) , 0];
Filename: runchem.m

%RUNCHEM run the chemical reaction example

tout = [0,5]; %%% output times

u0 = [1;0;0]; %%% initial data

options = odeset(’reltol’,1.e-5,... %%% solver parameters

’abstol’,1.e-8,...
’Jacobian’,@chemjac);

[t,u] = ode15s(@chemode,tout,u0,options); %%% run the solver

c2 = u(:,2); %%% concentration c2

subplot(1,2,1) %%% plot the results

plot(t,c2,’r-’,t,c2,’bx’)
axis([-0.1,5,0,4e-5])
ylabel(’Concentration of chemical 2’,’fontsize’,14)
xlabel(’time t’,’fontsize’,14)

subplot(2,2,2) %%% focus on small t

plot(t,c2,’r-’,t,c2,’bx’)
axis([-0.01,0.05,3e-5,4e-5])

subplot(2,2,4) %%% focus on small t

plot(t,c2,’r-’,t,c2,’bx’)
axis([0,0.01,3.5e-5,3.7e-5])

xlabel(’time t’,’fontsize’,14)

Table 2.3: M-files for Example 2.3

.
SMST C: Numerical Methods 2–15

−5 −5
x 10 x 10
4 4

3.8
3.5
3.6

3 3.4
Concentration of chemical 2

3.2
2.5

0 0.02 0.04

2
−5
x 10
3.7
1.5

3.65

1
3.6

0.5
3.55

0 3.5
0 1 2 3 4 5 0 0.002 0.004 0.006 0.008 0.01
time t time t

Figure 2.5: Results of Example 2.3. The markers show the time levels used by the adaptive time stepping
in the solver, and the smaller graphs focus on the region of small time steps. There is a fast transient at
the start. The problem is fairly stiff, and explicit methods do not work well for it.
SMST C: Numerical Methods 2–16

2.7 The Method of Lines for simple, time-dependent PDEs

The method of lines (MOL) is a simple approach to approximating PDEs that depend on both space
and time variables. Usually a finite difference, finite element or spectral approximation is made in the
space direction(s), and this converts the PDE into a system of ODEs suitable for attack by our methods.
Although popular, it can go badly wrong for all but the most cooperative PDEs and we leave the difficult
ones to later in this stream.
One thing to bear in mind is that the tolerances reltol and abstol used in ODE solvers only apply to the
accuracy of the solution of the ODE system. They say nothing about the accuracy of approximating the
PDE by a system of ODEs, and that is much bigger topic covered later in the stream. If the ODEs do
not provide a good approximation of the PDE, then there is no point in trying to approximate them very
well. A good approximation of a bad approximation is a bad approximation overall.
Let us look at simple reaction-diffusion problems which can often be tackled by the MOL. A common
example problem is
∂u ∂2u
= ε2 2 + g(u) , g(u) = u − u2
∂t ∂x
u(0, t) = 1, u(1, t) = 0, t > 0 , u(x, 0) = α(x), x ∈ (0, 1) .
The solution is approximated by

uj (t) ≈ u(xj , t) where xj = jh for j = 0, . . . , N

with space mesh size

interval length 1
h= = .
N N
The second space derivative uxx (x, t) is approximated using a finite difference. Taylor expansion for small
h gives

h2 h3
 
u(xj − h, t) − 2u(xj , t) + u(xj + h, t) = u − hux + uxx − uxxx + . . .
2! 3! (xj ,t)
−2u(xj , t)
h2 h3
 
+ u + hux + uxx + uxxx + . . .
2! 3! (xj ,t)

2h4
= h2 uxx (xj , t) + uxxxx (xj , t) + O(h6 )
4!
and so
u(xj − h, t) − 2u(xj , t) + u(xj + h, t) uj−1 (t) − 2uj (t) + uj+1 (t)
uxx (xj , t) ≈ ≈ .
h2 h2
The other parts of the PDE are easily approximated:
duj ∂u
(t) ≈ (xj , t) and g(uj (t)) ≈ g(u(xj , t))
dt ∂t
are obvious (but not unique) choices.
We put these three approximations together to define the MOL approximation
 
duj 2 uj−1 − 2uj + uj+1
=ε + g(uj )
dt h2

at the interior points j = 1, . . . , N−1, with u0 (t) = 1, uN (t) = 0. Now we write this in the standard form

u′ = f (t, u) ≡ Au + g(u) + b

where

u = (u1 , u2 , . . . , uN −1 )T , g(u) = (g(u1 ), g(u2 ), . . . , g(uN −1 ))T , b = (ε2 /h2 , 0, . . . , 0)T

SMST C: Numerical Methods 2–17

are all vectors of length N −1, and the N −1 × N −1 matrix

 
−2 1
. .
ε2  1 .. ..
 

A= 2  
h  .. .
. ..

1 
1 −2

is tridiagonal.

Specify the ODEs — function M-file rdodes.m

This specifies the ODE function f (t, u) = Au + g(u) + b, including the system size parameter N , matrix
A and boundary condition vector b as arguments of the function. This saves on recomputing them each
time the function is called.

function dudt = rdodes(t,u,N,A,b)

%RDODES example function for reaction-diffusion in Matlab
% ODE system is du/dt = Au + u-u^2 + b
% A is a N-1 x N-1 tridiagonal matrix

dudt = A*u + u - u.^2 + b; %%% ODE system

Specify the Jacobian — function M-file rdjac.m

Similar to the previous function, but this time specifying the Jacobian matrix
∂f
= A + diag(1 − 2uj : j = 1, . . . , N − 1),
∂u
which is a sparse N − 1 × N − 1 tridiagonal matrix.

function dfdu = rdjac(t,u,N,A,b)

%RDJAC Jacobian matrix for reaction-diffusion in Matlab
% ODE system is du/dt = Au + u-u^2 + b
% A is a tridiagonal matrix

dfdu = A + spdiags(1-2*u,0,N-1,N-1); %%% sparse matrix

Set up and solve — script M-file runrdodes.m

%RUNRDODES run the reaction-diffusion example

% u_t = epsilon^2 u_xx + u-u^2
% u(0,t) = 1, u(1,t) = 0, u(x,0) = 0

%%%---------------------------------------------------------
%%% Problem settings
%%%---------------------------------------------------------

epsilon = 0.05; %%% physical property

tspan = [0:0.1:10]; %%% output times

N = 400; %%% space approximation settings

dx = 1/N; %%% space mesh size

RTOL = 1.e-4; %%% ODE solver settings

ATOL = 1.e-8;

%%%---------------------------------------------------------
SMST C: Numerical Methods 2–18

%%% Set up matrix A = tridiag(1,-2,1) * epsilon^2/dx^2,

%%% the boundary information b and the initial data
%%%---------------------------------------------------------

A = sparse(N-1,N-1); %%% make A sparse

A(1,1:2) = [-2,1]; %%% define it row by row

for j=2:N-2
A(j,j-1:j+1) = [1,-2,1];
end
A(N-1,N-2:N-1) = [1,-2];

A = (epsilon^2/dx^2)*A;

b = [(epsilon^2/dx^2);zeros(N-2,1)]; %%% boundary info

u = zeros(N-1,1); %%% initial data

%%%---------------------------------------------------------
%%% Run the ODE solvers ode45 and ode15s
%%%---------------------------------------------------------

disp(’-------------------------------------’);
disp(’TEST ODE45’)
disp(’-------------------------------------’);

options=odeset(’stats’,’on’,’AbsTol’,ATOL,’RelTol’,RTOL);

tic %%% start timer

[tout,yout] = ode45(@rdodes,tspan,u,options,N,A,b);
t1 = toc; %%% stop timer

disp(’ ’)
disp([’ODE solve by ode45 took ’,num2str(t1),’ seconds’]);
disp(’ ’)

%%%----------------------------------------------------------

disp(’-------------------------------------’);
disp(’TEST ODE15S’)
disp(’-------------------------------------’);

options=odeset(’stats’,’on’,’AbsTol’,ATOL,’RelTol’,RTOL,’Jacobian’, @rdjac);

tic %%% start timer

[tout,yout] = ode15s(@rdodes,tspan,u,options,N,A,b);
t2 = toc; %%% stop timer

disp(’ ’)
disp([’ODE solve by ode15s took ’,num2str(t2),’ seconds’]);

2.7.1 Output and discussion

The output from runrdodes is below:

>> runrdodes
-------------------------------------
TEST ODE45
SMST C: Numerical Methods 2–19

-------------------------------------
4839 successful steps
314 failed attempts
30919 function evaluations

ODE solve by ode45 took 4.7071 seconds

-------------------------------------
TEST ODE15S
-------------------------------------
246 successful steps
0 failed attempts
291 function evaluations
1 partial derivatives
38 LU decompositions
289 solutions of linear systems

ODE solve by ode15s took 0.65984 seconds

We see that the implicit code ode15s was significantly faster than the explicit code ode45. The rest of
the diagnostic output above is produced as a result of setting stats to on in the options statement. We
see that ode45 did no linear algebra (since it is explicit), but used a large number of function evaluations.
On the other hand, ode15s did some linear algebra (LU decompositions and solutions of linear systems)
and a small number of function evaluations and the net effect is that it took much less time overall.
Both codes use schemes of roughly the same order of accuracy (order 4 and order up to 5 respectively), so
there is no obvious difference there. The main reason that ode15s is more efficient is that this problem is
stiff, and ode15s is able to take fewer, bigger steps than ode45, which is hampered by its linear stability
restriction. Indeed we can see that ode45 shows quite a few “failed attempts”. These are steps which
have to be done again because the step size predicted by the local time step control turned out to be too
big — a common indication that the problem is stiff.

2.8 What has not been covered?

A lot!

ˆ Differential Algebraic Equations (DAEs) — these are ODEs combined with algebraic constraints
(see e.g. [1]). Matlab and other ODE software deals with this reasonably easily.

ˆ Delay Differential Equations — ODEs with a delay term. Some software (Matlab dde23 and
ddesd), but not so easy.
ˆ Event detection — at the simplest level, work out when the solution passes through a target value.
Matlab and other ODE software deals with this reasonably easily.

ˆ Hamiltonian systems of any kind (e.g. molecular dynamics, N-body gravtitional problems) and
systems with special invariants.
ˆ Highly oscillatory systems.

ˆ Boundary value problems — similar to elliptic PDE problems coming in later lectures.

ˆ Stochastic differential equations — coming in later lectures. Take care – don’t use standard ODE
solvers with stochastic forcing of any kind.

References
[1] U M Ascher and L R Petzold, Computer Methods for Ordinary Differential Equations and Differential-
Algebraic Equations, SIAM (1998).
SMST C: Numerical Methods 2–20

[2] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II, Stiff and Differential Algebraic
Equations. Number 14 in Springer Computational Mathematics. Springer, Berlin, 1991.
[3] D J Higham and N J Higham, Matlab Guide, SIAM (2005).
[4] K Segeth, A posteriori error estimation with the finite element method of lines for a nonlinear parabolic
equation in one space dimension, Numer. Math. 83, 455-475 (1999).

[5] L F Shampine and M W Reichelt, The MATLAB ODE suite, SIAM J. Sci. Comput., 18, 1–22 (1997).