.

Definitions

Numerical methods for ordinary differential equations approximate solutions to initial value problems of the form

\( y' = f(t,y), \quad y(t_0) = y_0. \)

The result is approximations for the value of y(t) at discrete times \( t_i \) :

\( y_i \approx y(t_i) \quad\text{where}\quad t_i = t_0 + i h, \)

where h is the time step (sometimes referred to as \Delta t ) and i is an integer.

Multistep methods use information from the previous s steps to calculate the next value. In particular, a linear multistep method uses a linear combination of \( y_i \) and \( f(t_i,y_i) \) to calculate the value of y for the desired current step. Thus, a linear multistep method is a method of the form

\( \begin{align} & y_{n+s} + a_{s-1} \cdot y_{n+s-1} + a_{s-2} \cdot y_{n+s-2} + \cdots + a_0 \cdot y_n \\ & \qquad {} = h\cdot\left( b_s \cdot f(t_{n+s},y_{n+s}) + b_{s-1} \cdot f(t_{n+s-1},y_{n+s-1}) + \cdots + b_0 \cdot f(t_n,y_n) \right), \end{align} \)

The coefficients \( a_0, \dotsc, a_{s-1} \) and \(b_0, \dotsc, b_s \) determine the method. The designer of the method chooses the coefficients, balancing the need to get a good approximation to the true solution against the desire to get a method that is easy to apply. Often, many coefficients are zero to simplify the method.

One can distinguish between explicit and implicit methods. If \( b_s = 0 \) , then the method is called "explicit", since the formula can directly compute \( y_{n+s} \) . If \( b_s \ne 0 \) then the method is called "implicit", since the value of \( y_{n+s} \) depends on the value of \( f(t_{n+s}, y_{n+s}) \) , and the equation must be solved for \( y_{n+s} \) . Iterative methods such as Newton's method are often used to solve the implicit formula.

Sometimes an explicit multistep method is used to "predict" the value of \( y_{n+s} \) . That value is then used in an implicit formula to "correct" the value. The result is a predictor–corrector method.
Examples

Consider for an example the problem

\( y' = f(t,y)=y, \quad y(0) = 1. \)

The exact solution is \( y(t) = \mathrm{e}^t . \)
One-step Euler

A simple numerical method is Euler's method:

\( y_{n+1} = y_n + hf(t_n, y_n). \, \)

Euler's method can be viewed as an explicit multistep method for the degenerate case of one step.

This method, applied with step size \( h = \tfrac{1}{2} \) on the problem y' = y , gives the following results:

\( \begin{align} y_1 &= y_0 + hf(t_0, y_0) = 1 + \tfrac{1}{2}\cdot1 = 1.5, \\ y_2 &= y_1 + hf(t_1, y_1) = 1.5 + \tfrac{1}{2}\cdot1.5 = 2.25, \\ y_3 &= y_2 + hf(t_2, y_2) = 2.25 + \tfrac{1}{2}\cdot2.25 = 3.375, \\ y_4 &= y_3 + hf(t_3, y_3) = 3.375 + \tfrac{1}{2}\cdot3.375 = 5.0625. \end{align} \)

Two-step Adams–Bashforth

Euler's method is a one-step method. A simple multistep method is the two-step Adams–Bashforth method

\( y_{n+2} = y_{n+1} + \tfrac{3}{2} hf(t_{n+1},y_{n+1}) - \tfrac{1}{2} hf(t_n,y_n). \)

This method needs two values, \( y_{n+1} \) and \( y_n \) , to compute the next value, \( y_{n+2} \) . However, the initial value problem provides only one value, \( y_0 = 1 \) . One possibility to resolve this issue is to use the y_1 computed by Euler's method as the second value. With this choice, the Adams–Bashforth method yields (rounded to four digits):

\( \begin{align} y_2 &= y_1 + \tfrac32 hf(t_1, y_1) - \tfrac12 hf(t_0, y_0) = 1.5 + \tfrac32\cdot\tfrac12\cdot1.5 - \tfrac12\cdot\tfrac12\cdot1 = 2.375, \\ y_3 &= y_2 + \tfrac32 hf(t_2, y_2) - \tfrac12 hf(t_1, y_1) = 2.375 + \tfrac32\cdot\tfrac12\cdot2.375 - \tfrac12\cdot\tfrac12\cdot1.5 = 3.7812, \\ y_4 &= y_3 + \tfrac32 hf(t_3, y_3) - \tfrac12 hf(t_2, y_2) = 3.7812 + \tfrac32\cdot\tfrac12\cdot3.7812 - \tfrac12\cdot\tfrac12\cdot2.375 = 6.0234. \end{align} \)

The exact solution at t = t_4 = 2 is \( \mathrm{e}^2 = 7.3891\ldots \) , so the two-step Adams–Bashforth method is more accurate than Euler's method. This is always the case if the step size is small enough.

Families of multistep methods

Three families of linear multistep methods are commonly used: Adams–Bashforth methods, Adams–Moulton methods, and the backward differentiation formulas (BDFs).
Adams–Bashforth methods

The Adams–Bashforth methods are explicit methods. The coefficients are \( a_{s-1}=-1 and a_{s-2} = \cdots = a_0 = 0 \) , while the \( b_j \) are chosen such that the methods has order s (this determines the methods uniquely).

The Adams–Bashforth methods with s = 1, 2, 3, 4, 5 are (Hairer, Nørsett & Wanner 1993, §III.1; Butcher 2003, p. 103):

\( \begin{align} y_{n+1} &= y_n + hf(t_n, y_n) , \qquad\text{(This is the Euler method)} \\ y_{n+2} &= y_{n+1} + h\left( \frac{3}{2}f(t_{n+1}, y_{n+1}) - \frac{1}{2}f(t_n, y_n) \right) , \\ y_{n+3} &= y_{n+2} + h\left( \frac{23}{12} f(t_{n+2}, y_{n+2}) - \frac{4}{3} f(t_{n+1}, y_{n+1}) + \frac{5}{12}f(t_n, y_n)\right) , \\ y_{n+4} &= y_{n+3} + h\left( \frac{55}{24} f(t_{n+3}, y_{n+3}) - \frac{59}{24} f(t_{n+2}, y_{n+2}) + \frac{37}{24} f(t_{n+1}, y_{n+1}) - \frac{3}{8} f(t_n, y_n) \right) , \\ y_{n+5} &= y_{n+4} + h\left( \frac{1901}{720} f(t_{n+4}, y_{n+4}) - \frac{1387}{360} f(t_{n+3}, y_{n+3}) + \frac{109}{30} f(t_{n+2}, y_{n+2}) - \frac{637}{360} f(t_{n+1}, y_{n+1}) + \frac{251}{720} f(t_n, y_n) \right) . \end{align} \)

The coefficients \( b_j \) can be determined as follows. Use polynomial interpolation to find the polynomial p of degree s-1 such that

\( p(t_{n+i}) = f(t_{n+i}, y_{n+i}), \qquad \text{for } i=0,\ldots,s-1. \)

The Lagrange formula for polynomial interpolation yields

\( p(t) = \sum_{j=0}^{s-1} \frac{(-1)^{s-j-1}f(t_{n+j}, y_{n+j})}{j!(s-j-1)!h^{s-1}} \prod_{i=0 \atop i\ne j}^{s-1} (t-t_{n+i}). \)

The polynomial p is locally a good approximation of the right-hand side of the differential equation y' = f(t,y) that is to be solved, so consider the equation y' = p(t) instead. This equation can be solved exactly; the solution is simply the integral of p. This suggests taking

\( y_{n+s} = y_{n+s-1} + \int_{t_{n+s-1}}^{t_{n+s}} p(t)\,dt. \)

The Adams–Bashforth method arises when the formula for p is substituted. The coefficients \( b_j \) turn out to be given by

\( b_{s-j-1} = \frac{(-1)^j}{j!(s-j-1)!} \int_0^1 \prod_{i=0 \atop i\ne j}^{s-1} (u+i) \,du, \qquad \text{for } j=0,\ldots,s-1. \)

Replacing f(t, y) by its interpolant p incurs an error of order hs, and it follows that the s-step Adams–Bashforth method has indeed order s (Iserles 1996, §2.1)

The Adams–Bashforth methods were designed by John Couch Adams to solve a differential equation modelling capillary action due to Francis Bashforth. Bashforth (1883) published his theory and Adams' numerical method (Goldstine 1977).

Adams–Moulton methods

The Adams–Moulton methods are similar to the Adams–Bashforth methods in that they also have \( a_{s-1} = -1\) and \( a_{s-2} = \cdots = a_0 = 0 \) . Again the b coefficients are chosen to obtain the highest order possible. However, the Adams–Moulton methods are implicit methods. By removing the restriction that \( b_s = 0 \) , an s-step Adams–Moulton method can reach order s+1 , while an s-step Adams–Bashforth methods has only order s.

The Adams–Moulton methods with s = 0, 1, 2, 3, 4 are (Hairer, Nørsett & Wanner 1993, §III.1; Quarteroni, Sacco & Saleri 2000):

\( \begin{align} y_n &= y_{n-1} + h f(t_n,y_n) , \qquad\text{(This is the backward Euler method)}\\ y_{n+1} &= y_n + \frac{1}{2} h \left( f(t_{n+1},y_{n+1}) + f(t_n,y_n) \right) , \qquad\text{(This is the trapezoidal rule)}\\ y_{n+2} &= y_{n+1} + h \left( \frac{5}{12} f(t_{n+2},y_{n+2}) + \frac{2}{3} f(t_{n+1},y_{n+1}) - \frac{1}{12} f(t_n,y_n) \right) , \\ y_{n+3} &= y_{n+2} + h \left( \frac{3}{8} f(t_{n+3},y_{n+3}) + \frac{19}{24} f(t_{n+2},y_{n+2}) - \frac{5}{24} f(t_{n+1},y_{n+1}) + \frac{1}{24} f(t_n,y_n) \right) , \\ y_{n+4} &= y_{n+3} + h \left( \frac{251}{720} f(t_{n+4},y_{n+4}) + \frac{646}{720} f(t_{n+3},y_{n+3}) - \frac{264}{720} f(t_{n+2},y_{n+2}) + \frac{106}{720} f(t_{n+1},y_{n+1}) - \frac{19}{720} f(t_n,y_n) \right) . \end{align} \)

The derivation of the Adams–Moulton methods is similar to that of the Adams–Bashforth method; however, the interpolating polynomial uses not only the points \( t_{n-1},\dots, t_{n-s} \) , as above, but also \( t_n \) . The coefficients are given by

\( b_{s-j} = \frac{(-1)^j}{j!(s-j)!} \int_0^1 \prod_{i=0 \atop i\ne j}^{s} (u+i-1) \,du, \qquad \text{for } j=0,\ldots,s. \)

The Adams–Moulton methods are solely due to John Couch Adams, like the Adams–Bashforth methods. The name of Forest Ray Moulton became associated with these methods because he realized that they could be used in tandem with the Adams–Bashforth methods as a predictor-corrector pair (Moulton 1926); Milne (1926) had the same idea. Adams used Newton's method to solve the implicit equation (Hairer, Nørsett & Wanner 1993, §III.1).

Backward differentiation formulas (BDF)

Main article: Backward differentiation formula

The BDF methods are implicit methods with \( b_{s-1} = \cdots = b_0 = 0 \) and the other coefficients chosen such that the method attains order s (the maximum possible). These methods are especially used for the solution of stiff differential equations.

Analysis

The central concepts in the analysis of linear multistep methods, and indeed any numerical method for differential equations, are convergence, order, and stability.

Consistency and order

The first question is whether the method is consistent: is the difference equation

\( \begin{align} & y_{n+s} + a_{s-1} y_{n+s-1} + a_{s-2} y_{n+s-2} + \cdots + a_0 y_n \\ & \qquad {} = h \bigl( b_s f(t_{n+s},y_{n+s}) + b_{s-1} f(t_{n+s-1},y_{n+s-1}) + \cdots + b_0 f(t_n,y_n) \bigr), \end{align} \)

a good approximation of the differential equation y' = f(t,y) ? More precisely, a multistep method is consistent if the local truncation error goes to zero faster than the step size h as h goes to zero, where the local truncation error is defined to be the difference between the result \( y_{n+s} \) of the method, assuming that all the previous values \( y_{n+s-1}, \ldots, y_n \) are exact, and the exact solution of the equation at time \( t_{n+s} \) . A computation using Taylor series shows out that a linear multistep method is consistent if and only if

\( \sum_{k=0}^{s-1} a_k = -1 \quad\text{and}\quad \sum_{k=0}^s b_k = s + \sum_{k=0}^{s-1} ka_k. \)

All the methods mentioned above are consistent (Hairer, Nørsett & Wanner 1993, §III.2).

If the method is consistent, then the next question is how well the difference equation defining the numerical method approximates the differential equation. A multistep method is said to have order p if the local error is of order \( O(h^{p+1}) \) as h goes to zero. This is equivalent to the following condition on the coefficients of the methods:

\( \sum_{k=0}^{s-1} a_k = -1 \quad\text{and}\quad q \sum_{k=0}^s k^{q-1} b_k = s^q + \sum_{k=0}^{s-1} k^q a_k \text{ for } q=1,\ldots,p. \)

The s-step Adams–Bashforth method has order s, while the s-step Adams–Moulton method has order s+1 (Hairer, Nørsett & Wanner 1993, §III.2).

These conditions are often formulated using the characteristic polynomials

\( \rho(z) = z^s + \sum_{k=0}^{s-1} a_k z^k \quad\text{and}\quad \sigma(z) = \sum_{k=0}^s b_k z^k. \)

In terms of these polynomials, the above condition for the method to have order p becomes

\( \rho(\mathrm{e}^h) - h\sigma(\mathrm{e}^h) = O(h^{p+1}) \quad \text{as } h\to 0. \)

In particular, the method is consistent if it has order one, which is the case if \( \rho(1)=0 \) and \( \rho'(1)=\sigma(1). \)

Stability and convergence

The numerical solution of a one-step method depends on the initial condition y_0 , but the numerical solution of an s-step method depend on all the s starting values, \( y_0, y_1, \ldots, y_{s-1} \) . It is thus of interest whether the numerical solution is stable with respect to perturbations in the starting values. A linear multistep method is zero-stable for a certain differential equation on a given time interval, if a perturbation in the starting values of size ε causes the numerical solution over that time interval to change by no more than Kε for some value of K which does not depend on the step size h. This is called "zero-stability" because it is enough to check the condition for the differential equation y' = 0 (Süli & Mayers 2003, p. 332).

If the roots of the characteristic polynomial ρ all have modulus less than or equal to 1 and the roots of modulus 1 are of multiplicity 1, we say that the root condition is satisfied. A linear multistep method is zero-stable if and only if the root condition is satisfied (Süli & Mayers 2003, p. 335).

Now suppose that a consistent linear multistep method is applied to a sufficiently smooth differential equation and that the starting values \( y_1, \ldots, y_{s-1} \) all converge to the initial value \( y_0 \) as \( h \to 0 \) . Then, the numerical solution converges to the exact solution as \( h \to 0 \) if and only if the method is zero-stable. This result is known as the Dahlquist equivalence theorem, named after Germund Dahlquist; this theorem is similar in spirit to the Lax equivalence theorem for finite difference methods. Furthermore, if the method has order p, then the global error (the difference between the numerical solution and the exact solution at a fixed time) is \( O(h^p) \) (Süli & Mayers 2003, p. 340).

Furthermore, if the method is convergent, the method is said to be strongly stable if z=1 is the only root of modulus 1. If it is convergent and all roots of modulus 1 are not repeated, but there is more than one such root, it is said to be relatively stable. Note that 1 must be a root for the method to be convergent; thus convergent methods are always one of these two.

To assess the performance of linear multistep methods on stiff equations, consider the linear test equation y' = λy. A multistep method applied to this differential equation with step size h yields a linear recurrence relation with characteristic polynomial

\( \pi(z; h\lambda) = (1 - h\lambda\beta_s) z^s + \sum_{k=0}^{s-1} (\alpha_k - h\lambda\beta_k) z^k = \rho(z) - h\lambda\sigma(z). \)

This polynomial is called the stability polynomial of the multistep method. If all of its roots have modulus less than one then the numerical solution of the multistep method will converge to zero and the multistep method is said to be absolutely stable for that value of hλ. The method is said to A-stable if it is absolutely stable for all hλ with negative real part. The region of absolute stability is the set of all hλ for which the multistep method is absolutely stable (Süli & Mayers 2003, pp. 347 & 348). For more details, see the section on stiff equations and multistep methods.

Example

Consider the Adams–Bashforth three-step method

\( y_{n+1} = y_n + h\left( {23\over 12} f(t_n, y_n) - {16 \over 12} f(t_{n-1}, y_{n-1}) + {5\over 12}f(t_{n-2}, y_{n-2})\right). \)

The characteristic equation is thus

\( \rho(z)=z^3-z^2=z^2(z-1)\, \)

which has roots z=0, 1, and the conditions above are satisfied. As z=1 is the only root of modulus 1, the method is strongly stable.

First and second Dahlquist barriers

These two results were proved by Germund Dahlquist and represent an important bound for the order of convergence and for the A-stability of a linear multistep method. The first Dahlquist barrier was proved in Dahlquist (1956) and the second in Dahlquist (1963).

First Dahlquist barrier

A zero-stable and linear q-step multistep method cannot attain an order of convergence greater than q + 1 if q is odd and greater than q + 2 if q is even. If the method is also explicit, then it cannot attain an order greater than q (Hairer, Nørsett & Wanner 1993, Thm III.3.5).

Second Dahlquist barrier

There are no explicit A-stable and linear multistep methods. The implicit ones have order of convergence at most 2. The trapezoidal rule has the smallest error constant amongst the A-stable linear multistep methods of order 2.

.

Linear multistep method