Linear multistep method
From Scholarpedia
| This article is undergoing 1 initial review (1 completed); It may contain inaccuracies and unapproved changes made by anonymous reviewers. | |||||||||||||||||||||
Author: Dr. Ernst Hairer, Department of Mathematics, University of Geneva, Switzerland
Author: Dr. Gerhard Wanner, Department of Mathematics, University of Geneva, Switzerland
Dr. Ernst Hairer accepted the invitation on 3 October 2008 (self-imposed deadline: 3 April 2009).
Linear multistep methods constitute an important class of numerical integrators for ordinary differential equations, and particular methods are well suited for solving non-stiff and stiff equations as well as Hamiltonian systems over long time intervals. In contrast to one-step methods, multistep methods increase efficiency by using the information from several previous solution approximations.
Contents |
Methods for nonstiff problems
Consider an initial value problem
- (1)
where the solution
can be scalar- or vector-valued. Let
be a sequence of grid
points, which for simplicity is supposed to be equidistant
with step size
. A numerical method is an algorithm that
yields approximations
to the solution
at the grid points.
Explicit Adams methods
These methods are introduced by J.C. Adams (1883) for solving practical problems
of capillary action.
They are based on the following idea:
suppose that
values
approximating
are known. Compute the
derivatives
- (2)
and replace in the integrated form of (1)
- (3)
the integrand
by the polynomial
interpolating the values (2).
Then evaluate the integral analytically and obtain the next approximation to the
solution,
. After advancing the scheme by one step, this procedure can be repeated to obtain
, and so on. See Figure 1 for an illustration using the
logistic growth equation
.
Newton's formula for polynomial interpolation can be written as
where
is the backward difference operator
,
, etc.
This
leads to the following equations (for more details see Hairer, Nørsett & Wanner (1993), p. 357-358):
- (5)
Implicit Adams methods
If
were known, the interpolation
polynomial could be stretched over the entire interval
and thereby
an improved precision would be obtained. The corresponding formulas have the form
- (6)
and represent an implicit equation for computing
.
The special cases
and
are the implicit Euler
method and the trapezoidal rule, respectively. They
are actually one-step methods.
Predictor-Corrector schemes
An elegant way of solving the implicit equations in (6) is
to predict a value
by the explicit
method (5) with the same value of
, and replace
in the corrector (6) the missing value
by
. This
is the commonly used implementation
of the implicit Adams methods for non-stiff problems.
The complete algorithm created by this idea is illustrated in Figure 2.
General formulation and consistency
For theoretical investigations it is convenient to use a general framework that comprises
all explicit and implicit Adams methods as well as the BDF schemes considered below.
Following the seminal work of Dahlquist (1956),
a linear multistep method (linear
-step method) is written in the form
- (7)
(notice that the indices are shifted when compared to the formulas (5) and (6)). The generating polynomials for its coefficients,
- (8)
play an important role in the theory for these methods.
Adams methods are characterized by the fact that
.
A linear multistep method is called
consistent if it reproduces the exact solution for the differential equation
, when exact
starting approximations are used.
This is equivalent to
and
.
If it is exact for the differential equations
with
,
the method is said to have order
. In terms of the generating polynomials, this means that
The accuracy of methods having the same order can be distinguished by their error constant
Stability and convergence
The basic theory for linear multistep methods has been developed by Dahlquist (1959) and is
presented in the classical book of Henrici (1962). The main difference to one-step methods
is that consistency alone does not guarantee convergence. A linear
multistep method is called stable (or "zero-stable"), if for the differential equation
all solutions
of the difference equation (7) remain bounded. This is equivalent to the
property that all zeros of
satisfy
with those lying on the unit circle being simple. The highest
attainable order of linear
-step methods is
. However, for
these
methods are of no practical use because of the following result.
First Dahlquist barrier. The order
of a stable linear
-step
method satisfies

if
is even,

if
is odd,



if
, in particular if the method is explicit.
Convergence theorem. If a linear multistep method (7) is stable and of order
, then
it is convergent of order
.
Convergence of order
means that for sufficiently accurate starting approximations
the global error
satisfies
on compact intervals
, where
. The constant symbolized by the
notation is bounded by
, where
is a Lipschitz constant of the vector field
.
Methods for stiff problems
Continuing the computation of the above logistic growth equation
, using
, with the explicit Adams process
of order 3 leads to a bad surprise.
Instead of approaching the steady-state solution
, the numerical solution exhibits
a violent instability phenomenon (see Figure 3).
Stiff differential equations are characterized by the fact that the use of explicit methods requires very small step sizes, and certain implicit methods (like the implicit Euler method) perform much better. Stiff differential equations arise frequently as singularly perturbed problems in chemical reaction systems and in electrical circuitry, and as space discretizations of parabolic partial differential equations.
Stiff stability analysis
The foregoing stability analysis for nonstiff problems, which is uniquely based on the polynomial
, is not sufficient for stiff differential equations. There are many
attempts for studying the stability of methods for stiff problems. The most popular approach,
initiated by Dahlquist (1963), is motivated by linearization of
in the neighborhood
of a stationary solution,
translation and, in the case of systems, diagonalization. It is based on
- (9)
where
represents an eigenvalue of the Jacobian of
and
can be a complex number. In the above case, we have
.
Method (7) applied to this equation gives
,
whose stability depends on the roots of
- (10)
The stability domain
is the set of complex numbers
, for which all roots of (10)
satisfy
with those lying on the unit circle being simple.
For the explicit Adams method with
, equation (10) becomes
- (11)
or
- (12)
An elegant way of avoiding the high degree equation in (11)
is to insert the limiting values
, for
real, into (12).
This leads to the root locus curve for
, parts of which must
form the boundary
of the stability domain
. These domains
are drawn for four explicit Adams methods in Figure 4 and shrink rapidly with increasing
.
In the example from the beginning of this section, the value of
was
, far beyond
the stability bound
.
Implicit Adams methods for
have finite stability domains as well,ƒ
but they are considerably larger than those for explicit Adams methods. However, if they are solved
by a predictor-corrector approach, a good deal of the stability is lost, as seen
in Figure 5 for the case
.
Backward differentiation formulas (BDF)
This class of methods was
discovered, together with the phenomenon of stiffness, by Curtiss & Hirschfelder (1952) in calculations for chemistry, and their extreme importance for stiff problems has been recognized since the work of Gear (1971). In contrast to the
Adams methods, where polynomial interpolations were used for the slopes
,
we here interpolate the
-values
and require for the interpolating
polynomial
the collocation condition
. The differentiation of Newton's interpolation formula (4)
then leads to the methods
- (13)
(
denotes again the backward difference operator)
or, when expanded,
- (14)
These are implicit methods and work, for
, significantly better for stiff
problems, as illustrated in the movie of Figure 6. For
they
become unstable.
Their stability domains, always the outside of the root locus curve,
extend to
and can be seen
in Figure 7.
For stiff problems, their implementation requires some Newton-like
procedures for the implicit equations, because a predictor-corrector
approach would destroy this good stability property.
For more details we refer to the entry on BDF methods.
These methods are furthermore a very important tool for the
treatment of differential-algebraic equations.
A-stability
The analytical solution of Dahlquist's test equation (9) is stable
if and only if
. A very desirable property of a numerical method
is that, under this condition, the numerical solution remains stable as well,
i.e., that
- (15)
This property is called A-stability (Dahlquist 1963). Inspecting the
above stability domains, we observe that the implicit Adams methods for
(which are the implicit Euler and trapezoidal method) as well as the
BDF methods for
(the first one is implicit Euler again) have this
property. A close inspection of the stability domains for the higher order
BDF methods shows that some parts close to the imaginary axis, which
correspond to oscillatory problems, do not produce stable solutions.
This observation is part of the famous second Dahlquist barrier (Dahlquist 1963):
Second Dahlquist barrier. An A-stable multistep method must be
of order
. Among all multistep methods of order 2, the trapezoidal rule
has the smallest error constant.
This severe order bound was the origin of numerous efforts for breaking
it by generalizing multistep methods in various directions: Enright's second
derivative methods, second derivative BDF methods, blended
multistep methods (Skeel & Kong), super future points (Cash 1980),
multistep collocation methods, multistep methods of Radau type, which all fall into the class of
general linear methods.
For more details, consult Section V.3 of Hairer & Wanner (1996).
Also these methods, if A-stable, have their order barrier, which was for many
years an open conjecture (Daniel-Moore 1968): An A-stable generalized
multistep method, with
implicit stages per step, must be of
order
. The smallest error constant for such methods is that of the corresponding Gauss-Runge-Kutta method.
Order stars.
An elegant idea for proving the second Dahlquist barrier and the Daniel-Moore conjecture, together
with many other stability results, is the consideration of the corresponding
order star
. Instead of the stability domain
,
where the complex roots
of (10)
were compared, in absolute values, to 1, we now compare them to
,
the absolute value of the exact solution, which they tend to approximate, i.e.,
. Since
is multi-valued, this set has to be considered on the associated Riemann surface.
The order star gives much insight into the structure of the function
, because
it transforms analytic properties into geometric
properties as follows: the star-shaped structure with
sectors at the origin indicates
that the method is of order
, and A-stability is characterized by the fact that poles
are in the right half-plane and the intersection of
with the imaginary axis is empty.
The order star for the BDF3 method is drawn in Figure 8, and explains that such a method of order 3
cannot be A-stable.
A
-stability.
This is a concept, weaker than A-stability, but nevertheless very useful. It requires that the stability
domain contains a sector with angle
in the negative half plane;
more precisely, a linear multistep method (7) is called A
-stable if
. BDF methods are A
-stable
with angle given by
Convergence
The interest of applying BDF methods (or other A-stable or A
-stable multistep methods)
relies in the fact that they can be applied to stiff differential equations with step sizes that are much larger than the inverse of a Lipschitz
constant of
. In this situation, classical convergence theory (for non-stiff problems) does not provide any information.
It is nevertheless possible to prove estimates
of the global error
(for sufficiently accurate starting approximations and for A
-stable methods of order
) in the following situations:
- singularly perturbed problems
: if the eigenvalues of
are in the negative half plane, then the constant in the error estimation
can be chosen independently of
.
- discretized parabolic equations and stiff differential equations satisfying a one-sided Lipschitz condition.
Precise statements and proofs can be found, for example, in Sections V.7, V.8, and VI.2 of the monograph Hairer & Wanner (1996).
Methods for Hamiltonian systems
For long-time integration of conservative differential equations special integrators are usually more appropriate. The study of such integrators is the main topic of geometric numerical integration (Hairer, Lubich & Wanner 2006), and there are interesting geometric integrators among multistep methods.
Second order Hamiltonian systems
An important class of conservative systems are second order differential equations
- (16)
with a constant, symmetric, positive definite mass matrix
, and a smooth
potential
.
This equation is Hamiltonian with energy
,
where
is the momentum. The exact flow is known to be symplectic,
and the energy is preserved along solutions.
Typical examples are
-body
systems, e.g., planetary motion in astronomy or simulations in molecular
dynamics.
Symmetric linear multistep methods
It is possible to write (16) as a first order system and to apply any multistep method considered above. Eliminating the velocity leads to a formula
- (17)
with coefficients
that are different from those above.
We again introduce generating polynomials
A linear multistep method
(17) has order
for differential equations (16) if
and it is stable if all zeros of
satisfy
with those lying on the unit circle being at most double.
Order
and stability guarantee convergence, and the global error
satisfies
on compact intervals
.
For the integration of Hamiltonian systems, we are mainly interested
in long-time integration. For example, the solution of the
harmonic oscillator
satisfies
and remains
bounded for all times. Can we have
a similar behavior for the numerical solution, when
is fixed and
? When applied to the harmonic oscillator, method (17)
becomes a linear difference relation with characteristic equation
. It turns out that near
energy conservation for the harmonic oscillator can be achieved only
if all roots have modulus one. This in turn implies that
the method has to be symmetric, i.e.,
- (18)
Lambert & Watson (1976) initiated the study of the
long-time behavior of symmetric multistep methods.
Quinlan & Tremaine (1990) constructed high order
methods (17) and applied them successfully to long-time
integrations of the outer solar system. One of their methods
of order
with
is given by
- (19)
Notice that the method is explicit (i.e.,
) and therefore
its implementation with constant steps is straight-forward. Due to possible
numerical resonances, symmetric multistep methods are recommended
only for computations with high accuracy.
Long-time behavior
For the study of qualitative and quantitative features of
the numerical solution on intervals that are far beyond the
applicability of classical convergence theory, backward error
analysis is the appropriate tool. It turns out that a symmetric method
has to be
-stable (i.e., apart
from the double root at
,
all zeros of
are simple and of modulus one).
It has been shown in Hairer & Lubich (2004) that symmetric,
-stable
multistep methods (17) of order
have a long-time behavior that it
very similar to that of symplectic one-step methods. In particular,
their numerical solutions satisfy:
- for a Hamiltonian system (16) the total energy
is conserved up to
on intervals of length
;
- quadratic first integrals of the form
(such as the angular momentum in
-body problems) are conserved up to
on intervals of length
(here, numerical approximations
for the derivative are computed by symmetric finite differences);
- for completely integrable Hamiltonian systems, action variables are in general conserved up to
and angle variables are bounded by
on intervals of length
.
Implementation and software
The efficient implementation of linear multistep methods is a nontrivial task. A code has to be able to automatically select the step size (which can be done either by formulas with grid-dependent coefficients or by using a Nordsieck vector formulation) and to choose an optimal order. Among the available software (most of it can be downloaded from http://www.netlib.org/ode/) let us mention the following:
- ODE/STEP is a widely used code for the solution of non-stiff differential equations. It is based on Adams methods and documented in the monograph of Shampine & Gordon (1975).
- DIVA is an Adams integrator by Fred Krogh. It can directly integrate second order differential equations.
- LSODE is the Livermore Solver of A.~Hindmarsh. It solves stiff or nonstiff differential equations of first order.
- DASSL is a differential-algebraic equation solver, based on BDF schemes (see Brenan, Campbell & Petzold 1989). The most recent version and various modifications can be downloaded from the web page of L.~Petzold at http://www.cs.ucsb.edu/~cse/software.html
- MEBDF is an implementation of modified extended BDF methods (Cash & Considine 1992). Various modifications can be obtained from http://www2.imperial.ac.uk/~jcash/
References
- Adams J.C. (1883) Appendix in: F. Bashforth (1883) An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid. With an explanation of the method of integration employed in constructing the tables which give the theoretical form of such drops". Cambridge Univ. Press.
- Brenan K.E., Campbell S.L. and Petzold L.R. (1989 and 1996) Numerical solution of initial-value problems in differential-algebraic equations. North-Holland, New York, and SIAM, Philadelphia.
- Cash J.R. and Considine S. (1992) An MEBDF code for stiff initial value problems. ACM Trans. Math. Software 18:142-158.
- Curtiss C.F. and Hirschfelder J.O. (1952) Integration of stiff equations. Proc. Nat. Acad. Sci. 38:235-243.
- Dahlquist G. (1956) Convergence and stability in the numerical integration of ordinary differential equations. Math. Scand. 4:33-53.
- Dahlquist G. (1959) Stability and error bounds in the numerical integration of ordinary differential equations. Trans. of the Royal Inst. of Techn., Nr. 130, Stockholm, Sweden.
- Dahlquist G. (1963) A special stability problem for linear multistep methods. BIT 3:27-43.
- Gear C.W. (1971) Numerical initial value problems in ordinary differential equations. Prentice Hall.
- Hairer E and Lubich C (2004) Symmetric multistep methods over long times. Numer. Math. 97:699-723.
- Hairer E, Lubich C. and Wanner, G. (2006) Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Second edition. Springer-Verlag, Berlin. [1]
- Hairer E, Nørsett S.P. and Wanner, G. (1993) Solving Ordinary Differential Equations I: Nonstiff Problems. Second edition. Springer-Verlag, Berlin. [2]
- Hairer E. and Wanner, G. (1996) Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Second edition. Springer-Verlag, Berlin. [3]
- Henrici P. (1962) Discrete Variable Methods in Ordinary Differential Equations. Volume 1, John Wiley & Sons, New York.
- Lambert J.D. and Watson I.A. (1976) Symmetric multistep methods for periodic initial value problems. J. Inst. Maths. Applics. 18:189-202.
- Quinlan G.D. and Tremaine S. (1990) Symmetric multistep methods for the numerical integration of planetary orbits. Astron. J. 100:1694-1700.
- Shampine L. and Gordon M. (1975) "Computer Solution of Ordinary Differential Equations: The Initial Value Problem." Freeman.
