---
title: Multistep Methods, Systems, and Stability
module: Numerical Methods
moduleNumber: 7
lessonNumber: 2
order: 702
summary: >
  One-step methods discard everything but the last point. Multistep methods fit
  a polynomial to several past values and integrate it forward: the explicit
  Adams–Bashforth formulas, the implicit and more accurate Adams–Moulton
  formulas, and predictor–corrector pairs that combine them. The same rules
  extend verbatim to systems in vector form. A separate concern is stability:
  round-off can dominate truncation, and stiff equations force a tiny step for
  stability even when accuracy would allow a large one.
topics: [Numerical Methods]
sources:
  - book: Boyce
    ref: "§8.4 Multistep Methods; §8.5 Systems of First-Order Equations"
  - book: Boyce
    ref: "§8.6 More on Errors; Stability"
  - book: Simmons
    ref: "Ch. 14 §74 Errors; §77 Systems"
draft: false
---

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}$.[^boyce-ms] 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

$$
\phi(t_{n+1}) - \phi(t_n) = \int_{t_n}^{t_{n+1}} \phi'(t)\, \d t
= \int_{t_n}^{t_{n+1}} f\bigl(t, \phi(t)\bigr)\, \d t.
$$

An **Adams method** approximates the integrand $\phi'(t)$ by an interpolating
polynomial $P_k(t)$ of degree $k$ through $k+1$ known data points, then
integrates that polynomial exactly.[^boyce-adams] Which points are used splits
the family in two.

$$
% caption: 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.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (-0.1,0) -- (6.4,0) node[right, black!70] {$t$};
\draw[->, black] (0,-0.1) -- (0,3.0) node[above, black!70] {$f$};
% shade the integration step
\fill[acc!10] (4.4,0) rectangle (5.6,3.0);
% interpolating polynomial through past slopes
\draw[very thick, domain=0.6:5.7, samples=60] plot (\x, {2.2 - 0.62*\x + 0.12*\x*\x});
% data points
\foreach \x/\lab in {1.0/{$f_{n-3}$}, 2.2/{$f_{n-2}$}, 3.4/{$f_{n-1}$}, 4.4/{$f_n$}} {
  \fill[black!75] (\x, {2.2 - 0.62*\x + 0.12*\x*\x}) circle (1.8pt);
  \draw[black] (\x,0.05) -- (\x,-0.05);
  \node[black!70, anchor=north, font=\scriptsize] at (\x,-0.06) {\lab};
}
% open point for f_{n+1}
\draw[acc, thick, fill=white] (5.6, {2.2 - 0.62*5.6 + 0.12*5.6*5.6}) circle (1.8pt);
\draw[black] (5.6,0.05) -- (5.6,-0.05);
\node[black!70, anchor=north, font=\scriptsize] at (5.6,-0.06) {$f_{n+1}$};
\node[black!70, anchor=south] at (5.0,2.55) {integrate here};
\end{tikzpicture}
$$

### 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.

> **Definition (Adams–Bashforth formulas).** Explicit multistep formulas from
> integrating a polynomial through past slopes. The second- and fourth-order
> members are
>
> $$
> y_{n+1} = y_n + h\left(\tfrac{3}{2} f_n - \tfrac{1}{2} f_{n-1}\right),
> $$
> $$
> y_{n+1} = y_n + \tfrac{h}{24}\bigl(55 f_n - 59 f_{n-1} + 37 f_{n-2} - 9 f_{n-3}\bigr),
> $$
>
> with local truncation errors proportional to $h^3$ and $h^5$. The first-order
> Adams–Bashforth formula reduces to Euler's method.

### 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}$.

> **Definition (Adams–Moulton formulas).** Implicit multistep formulas whose
> interpolation includes the new point. The second- and fourth-order members are
>
> $$
> y_{n+1} = y_n + \tfrac{h}{2}\bigl(f_{n+1} + f_n\bigr),
> $$
> $$
> y_{n+1} = y_n + \tfrac{h}{24}\bigl(9 f_{n+1} + 19 f_n - 5 f_{n-1} + f_{n-2}\bigr),
> $$
>
> with local truncation errors proportional to $h^3$ and $h^5$. The first-order
> Adams–Moulton formula is the backward Euler method.

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.[^boyce-am] 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.

$$
% caption: 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.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize,
  box/.style={draw, minimum width=27mm, minimum height=11mm, align=center}]
\definecolor{acc}{HTML}{4A6FA5}
\node[box] (known) at (0,0) {known slopes\\$f_{n-3}, \dots, f_n$};
\node[box, draw=acc, text=acc] (pred) at (3.9,0) {predict $y_{n+1}^{P}$\\(Adams--Bashforth)};
\node[box] (eval) at (7.8,0) {evaluate\\$f_{n+1}$};
\node[box, draw=acc, text=acc] (corr) at (7.8,-2.0) {correct $y_{n+1}$\\(Adams--Moulton)};
\node[box] (adv) at (3.9,-2.0) {advance\\$n \gets n+1$};
\draw[->, thick] (known) -- (pred);
\draw[->, thick] (pred) -- (eval);
\draw[->, thick] (eval) -- (corr);
\draw[->, thick] (corr) -- (adv);
\draw[->, thick] (adv) -- (known);
\end{tikzpicture}
$$

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

> **Worked example.** Advance $y' = 1 - t + 4y$, $y(0) = 1$ from $t = 0.3$ to
> $t = 0.4$ with $h = 0.1$, using fourth-order Runge–Kutta starting values
> $y_0 = 1$, $y_1 = 1.6089333$, $y_2 = 2.5050062$, $y_3 = 3.8294145$ and the
> corresponding slopes $f_0 = 5$, $f_1 = 7.3357332$, $f_2 = 10.820025$,
> $f_3 = 16.017658$.
>
> Predict with the fourth-order Adams–Bashforth formula:
>
> $$
> y_4^{P} = y_3 + \tfrac{h}{24}\bigl(55 f_3 - 59 f_2 + 37 f_1 - 9 f_0\bigr) = 5.7836305.
> $$
>
> Evaluate the slope there, $f_4 = 1 - 0.4 + 4(5.7836305) = 23.734522$, then
> correct with the fourth-order Adams–Moulton formula, now explicit because
> $f_4$ is known:
>
> $$
> y_4 = y_3 + \tfrac{h}{24}\bigl(9 f_4 + 19 f_3 - 5 f_2 + f_1\bigr) = 5.7926721.
> $$
>
> The exact value is $\phi(0.4) = 5.7942260$. The predictor alone is off by
> $-0.0105955$; one correction cuts that to $-0.0015539$, a factor of about $7$
> from a single extra evaluation of $f$.

On the same problem the four fourth-order approaches rank cleanly by accuracy
and cost.[^boyce-pcex]

| 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** $\phi(t)$ rather than its
derivative, differentiates the polynomial, and sets $P_k'(t_{n+1}) = f(t_{n+1},
y_{n+1})$.

> **Definition (Backward differentiation formulas).** Implicit multistep formulas
> obtained by differentiating an interpolant of $y$ and matching the slope at the
> new point. The first-, second-, and fourth-order members are
>
> $$
> y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}),
> $$
> $$
> y_{n+1} = \tfrac{1}{3}\bigl(4 y_n - y_{n-1} + 2h f(t_{n+1}, y_{n+1})\bigr),
> $$
> $$
> y_{n+1} = \tfrac{1}{25}\bigl(48 y_n - 36 y_{n-1} + 16 y_{n-2} - 3 y_{n-3} + 12 h f(t_{n+1}, y_{n+1})\bigr).
> $$

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.[^boyce-bdf]

## 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), \qquad y' = g(t, x, y), \qquad x(t_0) = x_0, \; y(t_0) = y_0
$$

is written $\mathbf{x}' = \mathbf{f}(t, \mathbf{x})$, $\mathbf{x}(t_0) =
\mathbf{x}_0$, with $\mathbf{x} = (x, y)^\top$ and $\mathbf{f} = (f, g)^\top$.
The formulas are unchanged apart from replacing $y$ by $\mathbf{x}$ and $f$ by
$\mathbf{f}$.[^boyce-sys] The Euler and Runge–Kutta steps become

$$
\mathbf{x}_{n+1} = \mathbf{x}_n + h\,\mathbf{f}_n, \qquad
\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{h}{6}\bigl(\mathbf{k}_{n1} + 2\mathbf{k}_{n2} + 2\mathbf{k}_{n3} + \mathbf{k}_{n4}\bigr),
$$

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.

> **Worked example.** Integrate the system
>
> $$
> x' = x - 4y, \qquad y' = -x + y, \qquad x(0) = 1, \; y(0) = 0,
> $$
>
> whose exact solution is $x = \tfrac{1}{2}(e^{-t} + e^{3t})$,
> $y = \tfrac{1}{4}(e^{-t} - e^{3t})$, out to $t = 0.2$. Two Euler steps with
> $h = 0.1$, writing $\mathbf{f} = (x - 4y,\; -x + y)$:
>
> $$
> \mathbf{x}_1 = (1, 0) + 0.1\,(1, -1) = (1.1, -0.1),
> $$
> $$
> \mathbf{x}_2 = (1.1, -0.1) + 0.1\,(1.5, -1.2) = (1.25, -0.22).
> $$
>
> A single Runge–Kutta step with $h = 0.2$ uses the four vector slopes
> $\mathbf{k}_{01} = (1, -1)$, $\mathbf{k}_{02} = (1.5, -1.2)$,
> $\mathbf{k}_{03} = (1.63, -1.27)$, $\mathbf{k}_{04} = (2.342, -1.58)$:
>
> $$
> \mathbf{x}_1 = (1, 0) + \tfrac{0.2}{6}\bigl(\mathbf{k}_{01} + 2\mathbf{k}_{02} + 2\mathbf{k}_{03} + \mathbf{k}_{04}\bigr) = (1.320067, -0.250667).
> $$
>
> The exact values are $(1.320425, -0.250847)$. Euler is off by about $5.3\%$
> and $12.3\%$; the single RK4 step is within one-tenth of a percent. RK4 uses
> twice the evaluations of the Euler run yet is roughly $200$ times more
> accurate.

## 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.[^boyce-err] 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.

> **Definition (Round-off error).** The error $R_n = y_n - Y_n$ between the exact
> value $y_n$ a formula would give in perfect arithmetic and the value $Y_n$
> actually computed at finite precision. The total error is bounded by $|E_n| +
> |R_n|$, truncation plus round-off.

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

$$
% caption: Truncation error falls and round-off error rises as the step size
% shrinks; their sum is least at $h_{\mathrm{opt}}$, below which more steps only
% accumulate more round-off.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (-0.1,0) -- (5.2,0) node[right, black!70] {$h$};
\draw[->, black] (0,-0.1) -- (0,3.7) node[above, black!70] {error};
% truncation: rises with h
\draw[black!70, very thick, domain=0.5:4.7, samples=40] plot (\x, {0.55*\x});
\node[black!70, anchor=west] at (4.0,2.55) {$E_n$ (truncation)};
% round-off: falls with h
\draw[very thick, domain=0.5:4.7, samples=40] plot (\x, {0.9/\x});
\node[anchor=west] at (0.55,2.35) {$R_n$ (round-off)};
% total = sum
\draw[acc, very thick, dashed, domain=0.5:4.7, samples=60] plot (\x, {0.55*\x + 0.9/\x});
\node[acc, anchor=south west] at (2.9,3.05) {total};
% h_opt marker
\draw[black, dashed] (1.28,0) -- (1.28,1.4);
\node[black!70, anchor=north] at (1.28,-0.06) {$h_{\mathrm{opt}}$};
\end{tikzpicture}
$$

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_{\text{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' = ry$ with constant $r$,
whose exact solution decays when $r < 0$. Applying each method and asking whether
a perturbation $\delta$ at step $n$ grows or shrinks by the next step gives an
**amplification factor**.

- **Euler.** $y_{n+1} = y_n(1 + rh)$, amplification $1 + rh$. Stable only when
  $|1 + rh| < 1$, which for $r < 0$ requires $h < 2/|r|$.
- **Backward Euler.** $y_{n+1} = y_n/(1 - rh)$, amplification $1/(1 - rh)$. 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.[^boyce-stab] Plotting the set of
complex $z = hr$ for which the amplification has magnitude below one gives the
method's **region of absolute stability**.

$$
% caption: The absolute-stability region of forward Euler is the disk
% $|1 + hr| < 1$ in the complex plane, a unit disk centred at $-1$; the step
% must keep $hr$ inside it, whereas backward Euler is stable across the whole
% left half-plane.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (-3.3,0) -- (1.4,0) node[right, black!70] {$\mathrm{Re}(hr)$};
\draw[->, black] (0,-1.9) -- (0,1.9) node[above, black!70] {$\mathrm{Im}(hr)$};
% stability disk centred at (-1,0), radius 1
\fill[acc!12] (-1,0) circle (1);
\draw[acc, thick] (-1,0) circle (1);
% center and endpoints
\fill[acc] (-1,0) circle (1.3pt);
\draw[black] (-2,0.06) -- (-2,-0.06) node[anchor=north, black!70] {$-2$};
\node[acc, anchor=south] at (-1,0.06) {stable};
\node[black, anchor=west] at (0.15,1.1) {unstable};
\end{tikzpicture}
$$

> **Definition (Stiff equation).** A problem for which stability forces a much
> smaller step size than accuracy alone would need. Explicit methods become
> impractical; implicit methods (backward Euler, the backward differentiation
> formulas) are the standard tools.

> **Worked example.** The problem $y' = -100y + 100t + 1$, $y(0) = 1$ has exact
> solution $y = e^{-100t} + t$. The transient $e^{-100t}$ 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.[^boyce-stiff]

$$
% caption: On the stiff problem $y' = -100y + 100t + 1$, the exact solution is
% essentially $y = t$, yet the Euler method with $h$ past the stability limit
% oscillates with growing amplitude.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\definecolor{warn}{HTML}{C0392B}
\draw[->, black] (-0.1,0) -- (5.2,0) node[right, black!70] {$t$};
\draw[->, black] (0,-1.9) -- (0,2.2) node[above, black!70] {$y$};
% exact: y = t (rescaled)
\draw[black!75, very thick] (0,0.2) -- (5.0,1.9);
\node[black!75, anchor=south] at (4.0,1.5) {exact $y \approx t$};
% unstable euler zigzag about the line, growing
\draw[warn, thick]
  (0.3,0.32) -- (0.7,0.05) -- (1.1,0.75) -- (1.5,-0.15) -- (1.9,1.1) -- (2.3,-0.55) -- (2.7,1.7) -- (3.1,-1.15) -- (3.5,2.05);
\foreach \p in {(0.3,0.32),(0.7,0.05),(1.1,0.75),(1.5,-0.15),(1.9,1.1),(2.3,-0.55),(2.7,1.7),(3.1,-1.15),(3.5,2.05)}
  \fill[warn] \p circle (1.3pt);
\node[warn, anchor=west] at (3.4,-1.35) {Euler, $h > 2/|r|$};
\end{tikzpicture}
$$

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-ms]: **Boyce**, §8.4 — one-step versus multistep methods; a multistep formula uses data at more than the last mesh point.
[^boyce-adams]: **Boyce**, §8.4 — Adams methods approximate $\phi'(t)$ by an interpolating polynomial through $k+1$ data points and integrate it across $[t_n, t_{n+1}]$.
[^boyce-am]: **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-pcex]: **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 + 4y$, $y(0) = 1$.
[^boyce-bdf]: **Boyce**, §8.4 — backward differentiation formulas interpolate $\phi(t)$, and their popularity for stiff equations traces to Gear's work in the 1970s.
[^boyce-sys]: **Boyce**, §8.5 — extension of the Euler and Runge–Kutta formulas to systems in vector form; also **Simmons**, §77.
[^boyce-err]: **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-stab]: **Boyce**, §8.6 — stability analysis of $y' = ry$: the Euler amplification $1 + rh$ requires $h < 2/|r|$ for $r<0$, while backward Euler's $1/(1-rh)$ is unconditionally stable.
[^boyce-stiff]: **Boyce**, §8.6, Example 2 — the stiff problem $y' = -100y + 100t + 1$, $y(0)=1$, on which Euler with $h = 0.025$ is unstable while backward Euler with $h = 0.1$ is accurate.