Multistep Methods - Executable Engineering: Interactive Numerical Methods in Python

Multistep methods exploit the history of previously computed solution points $(t_{\le i}, y_{\le i})$ when advancing to the next step $y_{i+1}$ . Because they rely on a multi-point history, multistep schemes are not self-starting and require bootstrapping (e.g., via single-step Runge-Kutta or variable-order ramp-up) for the first few steps.

Starting from the exact integral formulation of an ODE:

\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t} &= f(t, y) \\ y_{i+1} &= y_i + \int_{t_i}^{t_{i+1}} f(t, y) \, \mathrm{d}t \end{aligned}

(1)

recall that numerical quadrature rules approximate an integral as a weighted sum of function evaluations at discrete nodes. Generalizing this idea across $s$ historical steps leads to the General Linear Multistep Formula:

\sum_{j=0}^{s} \alpha_j y_{i+j} = h \sum_{j=0}^{s} \beta_j f(t_{i+j}, y_{i+j})

(2)

where the coefficients $\alpha_j$ and $\beta_j$ are selected to maximize order of accuracy and numerical stability. By convention, we solve for the newest state $y_{i+s}$ and normalize $\alpha_s = 1$ .

Explicit Linear Multistep Methods¶

When $\beta_s = 0$ , the new solution $y_{i+s}$ depends exclusively on previously known states and derivatives, making the scheme explicit (such as the Adams-Bashforth family).

Explicit multistep formulae are frequently paired with an implicit multistep formula (such as Adams-Moulton) in predictor-corrector architectures: the explicit step generates an initial prediction based on past trajectory history, and the implicit corrector refines it using the slope over the new interval.

On their own, explicit multistep methods struggle with variable step sizes, cannot self-start, and have smaller stability regions than single-step Runge-Kutta schemes of the same order.

Backward Differentiation Formulae (BDF)¶

Setting $\beta_{j \neq s} = 0$ and $\beta_s \neq 0$ yields the most important class of implicit linear multistep methods: the Backward Differentiation Formulae (BDF), which serve as the default stiff ODE solvers in modern scientific software (method='BDF' in scipy.integrate.solve_ivp).

Starting from $\frac{\mathrm{d}y}{\mathrm{d}t} = f(t, y)$ , the BDF coefficients $\alpha_j$ are derived by differentiating a Lagrange interpolating polynomial fitted to the history $(t_n, y_n), \dots, (t_{n+1}, y_{n+1})$ and matching the derivative at $t_{n+1}$ .

For a uniform step size $h$ , the first five BDF orders are:

\begin{aligned} \text{BDF1:} \quad y_{n+1} &= y_n + h f(t_{n+1}, y_{n+1}) \\ \text{BDF2:} \quad y_{n+1} &= \frac{4}{3} y_n - \frac{1}{3} y_{n-1} + \frac{2}{3} h f(t_{n+1}, y_{n+1}) \\ \text{BDF3:} \quad y_{n+1} &= \frac{18}{11} y_n - \frac{9}{11} y_{n-1} + \frac{2}{11} y_{n-2} + \frac{6}{11} h f(t_{n+1}, y_{n+1}) \\ \text{BDF4:} \quad y_{n+1} &= \frac{48}{25} y_n - \frac{36}{25} y_{n-1} + \frac{16}{25} y_{n-2} - \frac{3}{25} y_{n-3} + \frac{12}{25} h f(t_{n+1}, y_{n+1}) \\ \text{BDF5:} \quad y_{n+1} &= \frac{300}{137} y_n - \frac{300}{137} y_{n-1} + \frac{120}{137} y_{n-2} - \frac{25}{137} y_{n-3} + \frac{12}{137} y_{n-4} + \frac{60}{137} h f(t_{n+1}, y_{n+1}) \end{aligned}

(3)

Notice that BDF1 is identical to the Backward Euler method! Key advantages of BDF solvers include:

Lagrange polynomial differentiation generalizes naturally to variable step sizes.
They are implicit with strong stability properties ideally suited for stiff equations (up to order 5; BDF6 is rarely used, and BDF orders $\ge 7$ are unconditionally unstable).
Solvers can dynamically vary the order $s \in \{1, \dots, 5\}$ to self-start from BDF1 and adapt whenever the physical regime changes.