5 Linear Multistep Methods
Linear multistep methods (LMMs) use information from several previous time levels to advance the numerical solution of an initial value problem
\bm{z}'(t)=\bm{f}(\bm{z}(t),t), \qquad \bm{z}(t_0)=\bm{z}_0.
On a uniform grid t_n=t_0+nh, a general k-step linear multistep method has the form
\boxed{ \sum_{j=0}^{k}\alpha_j\bm{z}_{n+j} = h\sum_{j=0}^{k}\beta_j\bm{f}_{n+j}, }
where
\bm{f}_{n+j}=\bm{f}(\bm{z}_{n+j},t_{n+j}),
and the coefficients \alpha_j,\beta_j are fixed scalars. We assume \alpha_k\neq0 and normalize the method so that \alpha_k=1 whenever convenient.
The method is explicit if \beta_k=0 and implicit if \beta_k\neq0.
Unlike one-step Runge–Kutta methods, a k-step method requires k starting values. These may be generated by a suitable one-step method of sufficiently high order.
5.1 Foundations
Characteristic Polynomials
Associated with an LMM are the two polynomials
\boxed{ \rho(\xi)=\sum_{j=0}^{k}\alpha_j\xi^j, \qquad \sigma(\xi)=\sum_{j=0}^{k}\beta_j\xi^j. }
They compactly encode consistency, zero-stability, and absolute stability.
For the scalar test equation
z'=\lambda z, \qquad q=h\lambda,
the recurrence becomes
\sum_{j=0}^{k}(\alpha_j-q\beta_j)z_{n+j}=0,
whose stability polynomial is
\boxed{ \pi(\xi;q)=\rho(\xi)-q\sigma(\xi). }
For a fixed q, numerical stability requires the roots of \pi(\xi;q) to satisfy the root condition: every root must lie in |\xi|\le1, and roots on the unit circle must be simple.
Consistency and Order
Insert the exact smooth solution into the multistep formula and define the local defect
\bm{d}_{n+k} = \sum_{j=0}^{k}\alpha_j\bm{z}(t_{n+j}) - h\sum_{j=0}^{k}\beta_j\bm{z}'(t_{n+j}).
Expanding about t_n gives
\bm{z}(t_n+jh) = \sum_{m=0}^{\infty} \frac{j^m h^m}{m!}\bm{z}^{(m)}(t_n).
Therefore the coefficient of h^m\bm{z}^{(m)}(t_n) is
\frac1{m!} \left( \sum_{j=0}^{k}\alpha_jj^m - m\sum_{j=0}^{k}\beta_jj^{m-1} \right), \qquad m\ge1,
while for m=0 it is \sum_j\alpha_j.
Hence an LMM has order at least p if and only if
\boxed{ \sum_{j=0}^{k}\alpha_j=0, }
and
\boxed{ \sum_{j=0}^{k}\alpha_jj^m = m\sum_{j=0}^{k}\beta_jj^{m-1}, \qquad m=1,\ldots,p. }
In terms of the characteristic polynomials, consistency (order at least one) is equivalent to
\boxed{ \rho(1)=0, \qquad \rho'(1)=\sigma(1). }
The first identity states that constant exact solutions are reproduced; the second states that linear exact solutions are reproduced to first order.
5.2 Adams methods
Adams methods are obtained by integrating
\bm{z}(t_{n+1}) = \bm{z}(t_n) + \int_{t_n}^{t_{n+1}}\bm{f}(\bm{z}(t),t)\,\mathrm dt
and replacing the integrand by an interpolation polynomial through previously computed values of \bm{f}.
Adams–Bashforth methods
Adams–Bashforth methods use only already known values
\bm{f}_n,\bm{f}_{n-1},\ldots,
and are therefore explicit.
Using a constant interpolation polynomial gives explicit Euler:
\boxed{ \bm{z}_{n+1}=\bm{z}_n+h\bm{f}_n } \qquad\text{(AB1, order 1).}
Linear interpolation through \bm{f}_n and \bm{f}_{n-1} gives
\boxed{ \bm{z}_{n+1} = \bm{z}_n +h\left( \frac32\bm{f}_n- \frac12\bm{f}_{n-1} \right) }
(AB2, order 2).
The next two members are
\boxed{ \bm{z}_{n+1} = \bm{z}_n +h\left( \frac{23}{12}\bm{f}_n - \frac{16}{12}\bm{f}_{n-1} + \frac{5}{12}\bm{f}_{n-2} \right) }
(AB3, order 3), and
\boxed{ \bm{z}_{n+1} = \bm{z}_n +h\left( \frac{55}{24}\bm{f}_n - \frac{59}{24}\bm{f}_{n-1} + \frac{37}{24}\bm{f}_{n-2} - \frac{9}{24}\bm{f}_{n-3} \right) }
(AB4, order 4).
Derivation of AB2
Set
t=t_n+\theta h, \qquad 0\le\theta\le1.
The linear interpolation polynomial through \theta=0 and \theta=-1 is
\bm{p}(\theta) =(1+\theta)\bm{f}_n-\theta\bm{f}_{n-1}.
Therefore
\int_{t_n}^{t_{n+1}}\bm{f}(\bm{z}(t),t)\,\mathrm dt \approx h\int_0^1\bm{p}(\theta)\,\mathrm d\theta = h\left( \frac32\bm{f}_n- \frac12\bm{f}_{n-1} \right),
which is AB2.
Adams–Moulton methods
Adams–Moulton methods include the unknown endpoint value \bm{f}_{n+1} in the interpolation formula and are therefore implicit.
The first members are
\boxed{ \bm{z}_{n+1}=\bm{z}_n+h\bm{f}_{n+1} }
(AM1, implicit Euler, order 1),
\boxed{ \bm{z}_{n+1} = \bm{z}_n + \frac h2 \left(\bm{f}_{n+1}+\bm{f}_n\right) }
(AM2, trapezoidal rule, order 2),
\boxed{ \bm{z}_{n+1} = \bm{z}_n +h\left( \frac5{12}\bm{f}_{n+1} + \frac23\bm{f}_n - \frac1{12}\bm{f}_{n-1} \right) }
(AM3, order 3), and
\boxed{ \bm{z}_{n+1} = \bm{z}_n +h\left( \frac38\bm{f}_{n+1} + \frac{19}{24}\bm{f}_n - \frac5{24}\bm{f}_{n-1} + \frac1{24}\bm{f}_{n-2} \right) }
(AM4, order 4).
Because \bm{f}_{n+1}=\bm{f}(\bm{z}_{n+1},t_{n+1}), an Adams–Moulton step generally requires a nonlinear solve.
Predictor–corrector implementation
A classical inexpensive implementation combines an Adams–Bashforth predictor with an Adams–Moulton corrector. For example, AB2 can predict
\bm{z}_{n+1}^{(0)} = \bm{z}_n +h\left( \frac32\bm{f}_n- \frac12\bm{f}_{n-1} \right),
then the trapezoidal rule corrects it via
\bm{z}_{n+1}^{(m+1)} = \bm{z}_n + \frac h2 \left[ \bm{f}_n+ \bm{f}(\bm{z}_{n+1}^{(m)},t_{n+1}) \right].
Repeated correction is a fixed-point iteration. For stiff problems, Newton-type solves are generally preferable.
5.3 Backward Differentiation Formulas (BDF)
BDF methods interpolate the solution values and differentiate the interpolation polynomial at the newest time. Consequently, only the newest right-hand-side value appears:
\boxed{ \sum_{j=0}^{k}\alpha_j\bm{z}_{n+j} = h\beta_k\bm{f}_{n+k}. }
After normalization, the first six BDF formulas may be written as follows.
BDF1
\boxed{ \bm{z}_{n+1}-\bm{z}_n =h\bm{f}_{n+1} }
(order 1; implicit Euler).
BDF2
\boxed{ \frac32\bm{z}_{n+1} -2\bm{z}_n +\frac12\bm{z}_{n-1} =h\bm{f}_{n+1} }
(order 2).
BDF3
\boxed{ \frac{11}{6}\bm{z}_{n+1} -3\bm{z}_n +\frac32\bm{z}_{n-1} -\frac13\bm{z}_{n-2} =h\bm{f}_{n+1} }
(order 3).
BDF4
\boxed{ \frac{25}{12}\bm{z}_{n+1} -4\bm{z}_n +3\bm{z}_{n-1} -\frac43\bm{z}_{n-2} +\frac14\bm{z}_{n-3} =h\bm{f}_{n+1} }
(order 4).
BDF5
\boxed{ \frac{137}{60}\bm{z}_{n+1} -5\bm{z}_n +5\bm{z}_{n-1} -\frac{10}{3}\bm{z}_{n-2} +\frac54\bm{z}_{n-3} -\frac15\bm{z}_{n-4} =h\bm{f}_{n+1} }
(order 5).
BDF6
\boxed{ \frac{49}{20}\bm{z}_{n+1} -6\bm{z}_n +\frac{15}{2}\bm{z}_{n-1} -\frac{20}{3}\bm{z}_{n-2} +\frac{15}{4}\bm{z}_{n-3} -\frac65\bm{z}_{n-4} +\frac16\bm{z}_{n-5} =h\bm{f}_{n+1} }
(order 6).
BDF methods are especially important for stiff problems. BDF1 and BDF2 are A-stable. BDF methods of orders 3–6 are not A-stable, but they retain large sectors of the left half-plane in their stability regions and are widely used in variable-order stiff integrators.
A crucial structural fact is that the classical BDF family is zero-stable only for orders 1 through 6; the formulas with k\ge7 violate the root condition.
5.4 Stability and Convergence
Zero-Stability
Consistency controls the defect made in a single idealized step. It does not by itself guarantee convergence: perturbations introduced at successive steps must not be amplified without bound as h\to0.
To isolate this mechanism, set h=0 in the LMM. The homogeneous recurrence is
\sum_{j=0}^{k}\alpha_j\bm{v}_{n+j}=\bm{0}.
Its behavior is determined by the first characteristic polynomial \rho.
Definition 5.1 (Definition: Zero-Stability) A linear multistep method is zero-stable if solutions of its homogeneous recurrence remain uniformly bounded on every fixed finite time interval when the starting values are uniformly bounded, independently of sufficiently small h.
The algebraic criterion is the fundamental root condition.
Theorem 5.1 (Root Condition) A linear multistep method is zero-stable if and only if every root \xi_i of
\rho(\xi)=0
satisfies
|\xi_i|\le1,
and every root with |\xi_i|=1 is simple.
Why the Root Condition Is Necessary
If \rho has a root |\xi|>1, the homogeneous recurrence contains a mode proportional to \xi^n, which grows geometrically. Since on a fixed interval n\sim T/h\to\infty as h\to0, arbitrarily small perturbations can be amplified without bound.
If a root on the unit circle has multiplicity m>1, the recurrence contains generalized modes of the form
n^r\xi^n, \qquad r=1,\ldots,m-1,
which also become unbounded as n\to\infty. Thus all unit-modulus roots must be simple.
Conversely, if all roots satisfy the stated condition, the companion recurrence matrix has powers that remain uniformly bounded; hence perturbations of the starting data remain bounded on finite intervals. This proves sufficiency.
Adams Methods Are Zero-Stable
All Adams–Bashforth and Adams–Moulton methods have
\rho(\xi)=\xi^k-\xi^{k-1} =\xi^{k-1}(\xi-1).
Thus the roots are 0 and 1, with 1 simple, so every Adams method is zero-stable.
BDF and Zero-Stability
For BDF methods, the root condition holds for k\le6 and fails for k\ge7. This is why the classical BDF family stops at order six in practical ODE solvers.
Convergence and the Dahlquist Equivalence Theorem
For an LMM, define the global error
\bm{e}_n=\bm{z}(t_n)-\bm{z}_n.
A method is convergent if, for every sufficiently regular well-posed IVP and consistent starting values,
\max_{0\le n\le N}\|\bm{e}_n\|\to0 \qquad\text{as }h\to0,
while Nh=T remains fixed.
The fundamental theorem of linear multistep analysis is the following.
Theorem 5.2 (Dahlquist Equivalence Theorem) For a linear multistep method, the following are equivalent:
- the method is convergent;
- the method is consistent and zero-stable.
Equivalently,
\boxed{ \text{consistency}+\text{zero-stability} \iff \text{convergence}. }
Proof Idea: Sufficiency
Insert the exact solution into the method. The global error satisfies a perturbed recurrence of the form
\sum_{j=0}^{k}\alpha_j\bm{e}_{n+j} = h\sum_{j=0}^{k}\beta_j \left[ \bm{f}(\bm{z}(t_{n+j}),t_{n+j}) - \bm{f}(\bm{z}_{n+j},t_{n+j}) \right] + \bm{d}_{n+k},
where \bm{d}_{n+k} is the local defect. If the method has order p, then
\|\bm{d}_{n+k}\|\le Ch^{p+1}.
Zero-stability gives a uniform bound on the discrete resolvent of the homogeneous recurrence. If \bm{f} is Lipschitz in its state argument, the nonlinear difference satisfies
\|\bm{f}(\bm{u},t)-\bm{f}(\bm{v},t)\| \le L\|\bm{u}-\bm{v}\|.
Applying the zero-stability bound and then a discrete Grönwall inequality yields
\max_{0\le n\le N}\|\bm{e}_n\| \le C_T \left( \max_{0\le j<k}\|\bm{e}_j\| +h^p \right).
Therefore consistent starting values and consistency of order p imply convergence of order p on every fixed finite interval.
Why Both Hypotheses Are Necessary
If a method converges for every sufficiently smooth IVP, it must reproduce constant and linear solutions asymptotically; this forces
\rho(1)=0, \qquad \rho'(1)=\sigma(1),
so the method is consistent.
If it were not zero-stable, arbitrarily small perturbations in the starting values could be amplified without bound even for the trivial equation
\bm{z}'=\bm{0}.
That would contradict convergence. Hence zero-stability is necessary as well.
Absolute Stability of Multistep Methods
For the test equation z'=\lambda z, a multistep method leads to
\pi(\xi;q)=\rho(\xi)-q\sigma(\xi), \qquad q=h\lambda.
The absolute-stability region is the set of all q\in\mathbb C such that every root \xi_i(q) of \pi(\xi;q)=0 satisfies
|\xi_i(q)|\le1,
with all roots on the unit circle simple.
The boundary of the region can often be parametrized by setting \xi=e^{i\theta}:
\boxed{ q(\theta) = \frac{\rho(e^{i\theta})}{\sigma(e^{i\theta})}, \qquad 0\le\theta\le2\pi, }
where \sigma(e^{i\theta})\neq0.
This formula is the multistep analogue of the Runge–Kutta stability function.
Dahlquist Barriers
Two classical impossibility results explain important design limitations of linear multistep methods.
- First Dahlquist barrier. A zero-stable explicit k-step linear multistep method has order at most k. More generally, a zero-stable k-step LMM satisfies p\le k+1 when k is odd and p\le k+2 when k is even.
- Second Dahlquist barrier. An A-stable linear multistep method has order at most two.
Consequently, no high-order classical LMM can be both A-stable and of arbitrary order. The trapezoidal rule and BDF2 attain order two while remaining A-stable; higher-order BDF methods sacrifice full A-stability for sectorial stability.
5.5 Practical Comparison
The three main families have complementary roles.
| Family | Explicit/implicit | Typical strength | Main limitation |
|---|---|---|---|
| Adams–Bashforth | explicit | low cost per step for non-stiff problems | restricted stability region |
| Adams–Moulton | implicit | higher accuracy using endpoint information | nonlinear solve may be required |
| BDF | implicit | strong performance on stiff dissipative problems | classical family zero-stable only through order 6 |
Linear multistep methods reuse previously computed information and can therefore require fewer new evaluations of \bm{f} per step than Runge–Kutta methods. Their main complications are startup, step-size changes, storage of history, and the additional zero-stability analysis absent from one-step methods.
Modern variable-step, variable-order stiff solvers are often based on BDF or closely related formulas, while modern non-stiff solvers frequently use embedded Runge–Kutta methods because adaptive step-size control and dense output are easier to organize in a one-step framework.