Prev
Unit 9 · Ordinary Differential Equations: Initial Value Problems
Multistep Methods
55m100 XP
Next
🛰️Ordinary Differential Equations: Initial Value Problems

Multistep Methods

Adams–Bashforth and Adams–Moulton, predictor–corrector, starting values, zero-stability and the root condition

Cambridge, 1883 — Aberdeen, 1918
In 1883 Francis Bashforth published An Attempt to Test the Theories of Capillary Action, comparing measured shapes of liquid drops with those predicted by the equations of surface tension. The drop profile is the solution of a differential equation, and Bashforth needed it to many digits, by hand. John Couch Adams — who had predicted the position of Neptune from irregularities in the orbit of Uranus — gave him a method. Instead of re-sampling the slope several times in each step, it fitted a polynomial through the slopes already computed at earlier steps and integrated that polynomial across the next step. Thirty-five years later, during the First World War, the American astronomer Forest Ray Moulton was asked to improve the range tables of the US Army's artillery at the Aberdeen Proving Ground. Shells experience drag that depends on speed, air density changing with height, and the rotation of the Earth; the old approximate formulas were not good enough. Moulton adapted Adams's methods, added an implicit corrector, and organised teams of human computers to produce trajectories. The predictor–corrector schemes he described in New Methods in Exterior Ballistics (1926) carry his name with Adams's.
A linear multistep method uses the last kk values to take the next step. The Adams methods integrate an interpolating polynomial through recent values of ff:
yn+1=yn+∫tntn+1p(t) dt.y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} p(t)\,dt.
If pp interpolates fn,fn−1,…f_n, f_{n-1}, \dots (known values), the method is explicit — Adams–Bashforth. If it also uses fn+1=f(tn+1,yn+1)f_{n+1} = f(t_{n+1}, y_{n+1}), the method is implicit — Adams–Moulton.
AB4: yn+1=yn+h24(55fn−59fn−1+37fn−2−9fn−3),AM3: yn+1=yn+h24(9fn+1+19fn−5fn−1+fn−2)\text{AB4: } y_{n+1} = y_n + \tfrac{h}{24}\big(55f_n - 59f_{n-1} + 37f_{n-2} - 9f_{n-3}\big),\qquad \text{AM3: } y_{n+1} = y_n + \tfrac{h}{24}\big(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}\big)

Fourth-order Adams–Bashforth and Adams–Moulton formulas. AB4 needs one new evaluation of f per step — against four for RK4 — but it needs four past values to start.

Derivation
Deriving Adams–Bashforth from an interpolating polynomial
Two steps. Pass a straight line through the two most recent slopes, (tn−1,fn−1)(t_{n-1}, f_{n-1}) and (tn,fn)(t_n, f_n):
p(t)=fn+fn−fn−1h (t−tn).p(t) = f_n + \frac{f_n - f_{n-1}}{h}\,(t - t_n).
Integrate it across the next step, from tnt_n to tn+1=tn+ht_{n+1} = t_n + h. The constant gives hfnhf_n; the linear term gives fn−fn−1h⋅h22\frac{f_n - f_{n-1}}{h}\cdot\frac{h^2}{2}:
yn+1=yn+hfn+h2(fn−fn−1)=yn+h(32fn−12fn−1).y_{n+1} = y_n + hf_n + \frac h2(f_n - f_{n-1}) = y_n + h\Big(\frac32f_n - \frac12f_{n-1}\Big).
That is AB2. The interpolation error of a line is O(h2)O(h^2) (lesson 7-1), integrated over a step of length hh: local error O(h3)O(h^3), global O(h2)O(h^2).
Four steps. The same procedure with the cubic through fn,fn−1,fn−2,fn−3f_n, f_{n-1}, f_{n-2}, f_{n-3} — Lagrange's formula, integrated term by term over [tn,tn+1][t_n, t_{n+1}] — produces the weights 124(55,−59,37,−9)\tfrac{1}{24}(55, -59, 37, -9). They add up to 24, as they must: for ff constant the step must be exactly hfhf. Including fn+1f_{n+1} among the interpolation points gives the implicit Adams–Moulton weights instead.
💡Predictor–corrector
The implicit Adams–Moulton formula is more accurate and more stable, but fn+1f_{n+1} depends on the unknown yn+1y_{n+1}. Avoid solving an equation: Predict yn+1y_{n+1} with Adams–Bashforth, Evaluate ff there, Correct with Adams–Moulton, Evaluate again. This PECE scheme costs two evaluations per step and is fourth order. On a circular orbit with the same budget, ABM4 is seven times more accurate than RK4 at 400 evaluations and nineteen times at 800. The difference between predictor and corrector also estimates the local error for free (Milne's device).
Starting values. A kk-step method needs kk values before it can start. They come from a one-step method of the same order — RK4 for AB4 — or, in variable-order codes, from starting at order 1 with tiny steps and raising the order as history accumulates. Changing the step size is also harder than for Runge–Kutta, since the formulas assume equal spacing; production codes (Adams in LSODE, solve_ivp's LSODA) handle both by working with interpolating polynomials directly.
💡Absolute stability: what a fixed step does to a decaying mode
Accuracy as h→0h \to 0 is one question. In practice hh is fixed, and a second question matters: does an error made now grow or shrink over many steps? To answer it once for all problems, apply each method to the test equation y′=λyy' = \lambda y with Re⁡λ<0\operatorname{Re}\lambda < 0 — a single decaying mode. (Lesson 9-4 shows why every linear system is a collection of such modes, one per eigenvalue.) A one-step method then multiplies yy by a fixed factor every step, its amplification factor R(z)R(z), which depends only on z=hλz = h\lambda. Explicit Euler gives yn+1=yn+hλyny_{n+1} = y_n + h\lambda y_n, so R(z)=1+zR(z) = 1 + z; RK4 gives the polynomial 1+z+z2/2+z3/6+z4/241 + z + z^2/2 + z^3/6 + z^4/24 found in lesson 9-2. The numerical solution decays, as the true one does, only where ∣R(z)∣≤1|R(z)| \le 1. That set of zz is the method's region of absolute stability — for Euler a disc of radius 1 centred at −1-1, reaching −2-2 on the real axis; for RK4 a larger region reaching −2.79-2.79. A multistep method has no single factor; its solutions behave like ζn\zeta^n for the roots of a polynomial, and the region is where all roots satisfy ∣ζ∣≤1|\zeta| \le 1 (see the derivation below). Here history has a price: AB2's region reaches only −1-1 on the real axis and AB4's only −0.3-0.3. High order from past values comes with a small stability region. The interactive draws all of them; keep an eye on how small the explicit ones are — the next lesson is about problems where that matters enormously.
Derivation
Stability of a multistep method on y′=λyy' = \lambda y
Apply AB2 to y′=λyy' = \lambda y: fn=λynf_n = \lambda y_n, so with z=hλz = h\lambda
yn+1=yn+z(32yn−12yn−1).y_{n+1} = y_n + z\Big(\frac32y_n - \frac12y_{n-1}\Big).
A linear recurrence like this has solutions of the form yn=ζny_n = \zeta^n. Substitute and divide by ζn−1\zeta^{n-1}:
ζ2−(1+32z)ζ+z2=0.\zeta^2 - \Big(1 + \frac32z\Big)\zeta + \frac z2 = 0.
The numerical solution is a combination of ζ1n\zeta_1^n and ζ2n\zeta_2^n, so it stays bounded exactly when both roots have ∣ζ∣≤1|\zeta| \le 1. On the real axis, a root reaches ζ=−1\zeta = -1 when 1+(1+32z)+z2=01 + (1 + \tfrac32z) + \tfrac z2 = 0, that is at z=−1z = -1: AB2 is stable for −1<z<0-1 < z < 0 — half of Euler's interval: the extra order bought from memory is paid for in stability. One of the two roots approximates eze^z; the other is a parasitic root that the method's memory creates, and it is the one that escapes first.
💡Zero-stability and the root condition
A general linear multistep method is ∑j=0kαjyn+j=h∑j=0kβjfn+j\sum_{j=0}^{k}\alpha_jy_{n+j} = h\sum_{j=0}^k\beta_jf_{n+j}, with characteristic polynomials ρ(ζ)=∑αjζj\rho(\zeta) = \sum\alpha_j\zeta^j and σ(ζ)=∑βjζj\sigma(\zeta) = \sum\beta_j\zeta^j. As h→0h \to 0 the method reduces to the recurrence ρ\rho — and its solutions behave like ζn\zeta^n for the roots ζ\zeta of ρ\rho. Zero-stability requires every root to satisfy ∣ζ∣≤1|\zeta| \le 1, with roots on the unit circle simple: the root condition. Dahlquist's equivalence theorem (1956): a consistent multistep method converges if and only if it is zero-stable. And his first barrier: a zero-stable kk-step method has order at most k+2k + 2 — you cannot buy arbitrary accuracy from history alone.
📐A consistent method that never converges
The two-step method yn+2+4yn+1−5yn=h(4fn+1+2fn)y_{n+2} + 4y_{n+1} - 5y_n = h(4f_{n+1} + 2f_n) is consistent — of order three, higher than AB2 — yet ρ(ζ)=ζ2+4ζ−5\rho(\zeta) = \zeta^2 + 4\zeta - 5 has roots 11 and −5-5. The parasitic root multiplies every error by −5-5 at every step. On y′=−yy' = -y it gives y(1)=−6.7y(1) = -6.7 with h=0.1h = 0.1, −4.7×106-4.7 \times 10^6 with h=0.05h = 0.05 and −2.9×1019-2.9 \times 10^{19} with h=0.025h = 0.025. Smaller steps make it worse, because there are more steps for the parasitic root to act. Consistency alone guarantees nothing.
Derivation
Why the roots of ρ\rho decide, as h→0h \to 0
As h→0h \to 0 the right-hand side h∑βjfn+jh\sum\beta_jf_{n+j} disappears and the method becomes the recurrence ∑jαjyn+j=0\sum_j\alpha_jy_{n+j} = 0. Try yn=ζny_n = \zeta^n:
∑jαjζn+j=ζn ρ(ζ)=0,\sum_j\alpha_j\zeta^{n+j} = \zeta^n\,\rho(\zeta) = 0,
so ζ\zeta must be a root of ρ\rho, and every solution is a combination ∑iciζin\sum_ic_i\zeta_i^n of them. (A double root contributes nζnn\zeta^n as well — which is why roots on the unit circle must be simple: n⋅1nn \cdot 1^n grows.)
One root is always 1. For the method to reproduce a constant solution — the exact solution of y′=0y' = 0 — we need ρ(1)=∑αj=0\rho(1) = \sum\alpha_j = 0; consistency guarantees it. That root carries the true solution. Every other root is created by the method's memory and carries nothing but errors. If any has ∣ζ∣>1|\zeta| > 1, those errors are multiplied by ∣ζ∣|\zeta| at every step, and since reaching a fixed time takes T/hT/h steps, the amplification ∣ζ∣T/h|\zeta|^{T/h} grows without bound as h→0h \to 0.
For the method above, ρ(ζ)=ζ2+4ζ−5=(ζ−1)(ζ+5)\rho(\zeta) = \zeta^2 + 4\zeta - 5 = (\zeta - 1)(\zeta + 5): the root 1, as required, and a parasitic −5-5. The twenty steps of h=0.05h = 0.05 multiply any error committed early on by 520≈10145^{20} \approx 10^{14} — and every step commits a new one.
You could not step twice into the same river.
– Heraclitus, as quoted in Plato's Cratylus
💡Where you will use this
Orbital Mechanics: high-order Adams and Störmer–Cowell multistep methods propagate orbits in operational orbit-determination software, because each step needs only one or two force evaluations — and force models with gravity fields of high degree are expensive. Digital control and Signals: a multistep method is a digital filter; zero-stability is the filter's stability. Ballistics and guidance: the predictor–corrector tradition continues in trajectory simulators.
📚Further reading
  • Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations (2nd ed., Cambridge University Press, 2009), chapter 2, 'Multistep methods' — order, the root condition and Dahlquist's theorems.
  • Ernst Hairer, Syvert Nørsett and Gerhard Wanner, Solving Ordinary Differential Equations I (2nd ed., Springer, 1993), chapter III, 'Multistep Methods and General Linear Methods'.
  • Germund Dahlquist, 'Convergence and stability in the numerical integration of ordinary differential equations', Mathematica Scandinavica 4, 33–53 (1956) — consistency plus zero-stability equals convergence.
↓But there is a problem
Adams methods, Runge–Kutta methods, all the explicit schemes of this unit share one more limitation. In 1952 two chemists at the University of Wisconsin, Charles Curtiss and Joseph Hirschfelder, were integrating the equations of chemical reactions in which some species react in microseconds while the overall process takes seconds. Their explicit solver crawled along with microsecond steps long after the fast reactions had finished and nothing was changing quickly any more. They called such equations stiff. Twenty years later, the circuit simulator SPICE at Berkeley met the same problem in every transistor.