Multistep Methods, Systems, and Stability

Euler, improved Euler, and Runge–Kutta are one-step methods: the value $y_{n + 1}$ depends only on the data at the single preceding point $(t_{n}, y_{n})$ . Once a few values $y_{1}, y_{2}, \dots, y_{n}$ have been computed, though, they are all available, and a multistep method uses several of them to compute $y_{n + 1}$ .¹ The saving is arithmetic: Runge–Kutta requires four evaluations of $f$ per step, whereas a multistep method of the same order reuses slopes it has already computed and needs only one or two new evaluations.

Throughout, the problem is $y^{'} = f (t, y)$ , $y (t_{0}) = y_{0}$ , the step size $h$ is constant, and $f_{j} = f (t_{j}, y_{j})$ .

Adams methods

The exact solution satisfies

ϕ (t_{n + 1}) - ϕ (t_{n}) = t_{n} \int t_{n + 1} ϕ^{'} (t) d t = t_{n} \int t_{n + 1} f (t, ϕ (t)) d t .

An Adams method approximates the integrand $ϕ^{'} (t)$ by an interpolating polynomial $P_{k} (t)$ of degree $k$ through $k + 1$ known data points, then integrates that polynomial exactly.² Which points are used splits the family in two.

An Adams method fits a polynomial through past slope values

f_{n - 3}, \dots, f_{n}

and integrates it across the next step

[t_{n}, t_{n + 1}]

(shaded); Adams–Moulton also uses the not-yet-known

f_{n + 1}

, making it implicit.

Adams–Bashforth: explicit

Interpolating through past points only, ending at $(t_{n}, f_{n})$ , gives an explicit formula. The degree-one case uses $(t_{n}, f_{n})$ and $(t_{n - 1}, f_{n - 1})$ ; integrating the line across $[t_{n}, t_{n + 1}]$ yields the second-order formula. Higher degrees use more past points.

Adams–Moulton: implicit

Including the new point $(t_{n + 1}, f_{n + 1})$ among the interpolation nodes gives an implicit formula, since $f_{n + 1} = f (t_{n + 1}, y_{n + 1})$ contains the unknown $y_{n + 1}$ .

At equal order the two families have errors of the same power of $h$ , but the Adams–Moulton constant is much smaller. For the fourth-order pair the Adams–Moulton proportionality constant is under one-tenth of the Adams–Bashforth constant.³ Adams–Moulton is more accurate, but each step must solve an equation for $y_{n + 1}$ .

Predictor–corrector methods

Combining the two families keeps the accuracy of the implicit formula while avoiding the implicit solve. Predict $y_{n + 1}$ with an explicit Adams–Bashforth step, evaluate $f_{n + 1}$ at the predicted value, then correct with the Adams–Moulton formula, which is now an explicit evaluation because $f_{n + 1}$ is known.

One predictor–corrector step: Adams–Bashforth predicts

y_{n + 1}^{P}

, a single evaluation gives

f_{n + 1}

, and Adams–Moulton corrects to

y_{n + 1}

before the window slides forward.

The corrector may be reapplied, but needing it more than once or twice signals that $h$ is too large and should be cut.

On the same problem the four fourth-order approaches rank cleanly by accuracy and cost.⁴

Method	Value at $t = 0.4$	Error	Cost per step
Adams–Bashforth (explicit)	$5.783631$	$- 0.010596$	$1$ evaluation
Predictor–corrector	$5.792672$	$- 0.001554$	$2$ evaluations
Backward differentiation (BDF4)	$5.796763$	$+ 0.002537$	implicit solve
Adams–Moulton (implicit)	$5.794268$	$+ 0.000042$	implicit solve

The exact value is $5.794226$ . Adams–Bashforth is fastest and least accurate; Adams–Moulton alone is most accurate but implicit; the predictor–corrector sits between, explicit yet far better than the predictor alone.

Backward differentiation formulas

A second implicit family interpolates the solution $ϕ (t)$ rather than its derivative, differentiates the polynomial, and sets $P_{k}^{'} (t_{n + 1}) = f (t_{n + 1}, y_{n + 1})$ .

The first-order BDF is again the backward Euler method. These formulas are the standard tool for stiff problems, discussed below, because of their stability, not their raw accuracy.⁵

Starting values

Every multistep method needs several values before it can begin: the fourth-order Adams–Bashforth formula needs $y_{1}, y_{2}, y_{3}$ in addition to $y_{0}$ . The usual remedy runs a one-step method of comparable order — typically fourth-order Runge–Kutta — for the first few steps, then hands off to the cheaper multistep formula. A one-step method is self-starting; a $k$ -step method is not.

The trade-off between the two styles has several axes.

Property	One-step (Runge–Kutta)	Multistep (Adams)
Evaluations of $f$ per step	four (RK4)	one or two (past a warm start)
Self-starting	yes	no; needs a one-step starter
Changing the step size	easy	awkward (fixed mesh assumed)
Error estimate	needs an embedded pair	easy from successive formulas
Off-mesh interpolation	not built in	polynomial gives it directly

Systems and higher-order equations

Every method extends to systems by promoting scalars to vectors. A system of two first-order equations

x^{'} = f (t, x, y), y^{'} = g (t, x, y), x (t_{0}) = x_{0}, y (t_{0}) = y_{0}

is written $x^{'} = f (t, x)$ , $x (t_{0}) = x_{0}$ , with $x = (x, y)^{⊤}$ and $f = (f, g)^{⊤}$ . The formulas are unchanged apart from replacing $y$ by $x$ and $f$ by $f$ .⁶ The Euler and Runge–Kutta steps become

x_{n + 1} = x_{n} + h f_{n}, x_{n + 1} = x_{n} + \frac{h}{6} (k_{n 1} + 2 k_{n 2} + 2 k_{n 3} + k_{n 4}),

with the four vector slopes defined exactly as before. This matters beyond systems proper: any $n$ th-order equation reduces to a first-order system by naming the derivatives as new variables, so a single vector integrator handles every scalar higher-order problem too.

Errors: truncation versus round-off

Order $p$ describes truncation error only. For a $p$ th-order method the local error is proportional to $h^{p + 1}$ and the global error on a finite interval to $h^{p}$ , so shrinking $h$ shrinks truncation error.⁷ But every arithmetic operation carries a round-off error from finite precision, and smaller $h$ means more steps, so accumulated round-off grows as $h$ shrinks.

The two errors move oppositely in $h$ , so the total error has a minimum at an optimal step $h_{opt}$ : below it, round-off dominates and further shrinking $h$ makes the answer worse.

Truncation error falls and round-off error rises as the step size shrinks; their sum is least at

h_{opt}

, below which more steps only accumulate more round-off.

For the test problem computed with only four significant digits, Euler's method is most accurate at $t = 1$ around $N \approx 1600$ steps; using more steps makes the answer worse as round-off takes over. Retaining more digits pushes $h_{opt}$ smaller, but the qualitative U-shape is unavoidable.

Stability and stiff equations

A numerical method can be stable or unstable independently of the problem. The cleanest illustration is the linear test equation $y^{'} = r y$ with constant $r$ , whose exact solution decays when $r < 0$ . Applying each method and asking whether a perturbation $δ$ at step $n$ grows or shrinks by the next step gives an amplification factor.

Euler. $y_{n + 1} = y_{n} (1 + r h)$ , amplification $1 + r h$ . Stable only when $∣1 + r h ∣ < 1$ , which for $r < 0$ requires $h < 2/∣ r ∣$ .
Backward Euler. $y_{n + 1} = y_{n} / (1 - r h)$ , amplification $1/ (1 - r h)$ . For $r < 0$ this is below $1$ for every $h$ : unconditionally stable.

So the explicit method inherits a step-size restriction that the problem itself does not impose; the implicit method does not.⁸ Plotting the set of complex $z = h r$ for which the amplification has magnitude below one gives the method's region of absolute stability.

The absolute-stability region of forward Euler is the disk

∣1 + h r ∣ < 1

in the complex plane, a unit disk centred at

- 1

; the step must keep

h r

inside it, whereas backward Euler is stable across the whole left half-plane.

Worked example. The problem $y^{'} = - 100 y + 100 t + 1$ , $y (0) = 1$ has exact solution $y = e^{- 100 t} + t$ . The transient $e^{- 100 t}$ is negligible past a thin boundary layer near $t = 0$ , after which the solution is essentially the line $y = t$ . Accuracy would permit a large step once past the layer, but stability of the Euler method needs $h < 2/∣ - 100 ∣ = 0.02$ . Comparing the Euler method just past that limit ( $h = 0.025$ ) with backward Euler ( $h = 0.1$ ) against the exact solution:

$t$	Exact	Euler, $h = 0.025$	Backward Euler, $h = 0.1$
$0.0$	$1.000000$	$1.000000$	$1.000000$
$0.1$	$0.100045$	$5.162500$	$0.190909$
$0.2$	$0.200000$	$25.8289$	$0.208264$
$0.4$	$0.400000$	$657.241$	$0.400068$

The Euler iterates oscillate with exploding amplitude and are worthless past the stability limit; backward Euler, unconditionally stable, tracks the solution cleanly at four times the step.⁹

On the stiff problem

y^{'} = - 100 y + 100 t + 1

, the exact solution is essentially

y = t

, yet the Euler method with

h

past the stability limit oscillates with growing amplitude.

The point is general: no numerical method makes an unstable problem stable, but a poorly chosen method or step can introduce instabilities absent from the original problem. Stiff systems appear wherever a model couples fast and slow modes, such as chemical kinetics and circuits with disparate time constants, and their efficient solution rests on implicit formulas rather than explicit ones.

Boyce, §8.4 — one-step versus multistep methods; a multistep formula uses data at more than the last mesh point. ↩
Boyce, §8.4 — Adams methods approximate $ϕ^{'} (t)$ by an interpolating polynomial through $k + 1$ data points and integrate it across $[t_{n}, t_{n + 1}]$ . ↩
Boyce, §8.4 — the fourth-order Adams–Moulton proportionality constant is under one-tenth that of Adams–Bashforth, so the implicit family is markedly more accurate at equal order. ↩
Boyce, §8.4, Examples 1–2 — comparison of the fourth-order Adams–Bashforth, Adams–Moulton, predictor–corrector, and backward differentiation results at $t = 0.4$ for $y^{'} = 1 - t + 4 y$ , $y (0) = 1$ . ↩
Boyce, §8.4 — backward differentiation formulas interpolate $ϕ (t)$ , and their popularity for stiff equations traces to Gear's work in the 1970s. ↩
Boyce, §8.5 — extension of the Euler and Runge–Kutta formulas to systems in vector form; also Simmons, §77. ↩
Boyce, §8.6 — for a method of order $p$ the local error is proportional to $h^{p + 1}$ and the global error to $h^{p}$ ; round-off error grows as $h$ shrinks, giving an optimal step size. Also Simmons, §74. ↩
Boyce, §8.6 — stability analysis of $y^{'} = r y$ : the Euler amplification $1 + r h$ requires $h < 2/∣ r ∣$ for $r < 0$ , while backward Euler's $1/ (1 - r h)$ is unconditionally stable. ↩
Boyce, §8.6, Example 2 — the stiff problem $y^{'} = - 100 y + 100 t + 1$ , $y (0) = 1$ , on which Euler with $h = 0.025$ is unstable while backward Euler with $h = 0.1$ is accurate. ↩