Linear multistep method

From Scholarpedia

This article has not been peer-reviewed or accepted for publication yet; It may be unfinished, contain inaccuracies, or unapproved changes.

Author: Dr. Ernst Hairer, Department of Mathematics, University of Geneva, Switzerland
Author: Dr. Gerhard Wanner, University of Geneva

Dr. Ernst Hairer accepted the invitation on 3 October 2008 (self-imposed deadline: 3 April 2009).

Linear multistep methods constitute an important class of numerical integrators, and particular methods are well suited for solving non-stiff and stiff differential equations as well as Hamiltonian systems over long time intervals.


Contents

[edit]

Methods for nonstiff problems

Let k\ge 1 be a fixed integer, t_0,t_1,t_2,\ldots be a sequence of grid points, which we suppose for the moment to be equidistant with step-size h=t_{n+1}-t_n, and

(1)
\dot y(t)=f(t,y(t))

a differential equation (or a system of differential equations) which we want to solve.

Figure 1: Explicit Adams methods
[edit]

Explicit Adams methods

These methods are based on the following idea: suppose that k values y_n,\,y_{n+1},\ldots, y_{n+k-1} approximating y(t_n),\, y(t_{n+1}),\ldots,y(t_{n+k-1}) are known, where n=0 for the first step. We compute the derivatives

(2)
f_{n+j}=f(t_{n+j},y_{n+j})\quad j=n,\ldots ,n+k-1

and replace in the integrated form of (1)

(3)
y(t_{n+k})=y(t_{n+k-1})+\int _{t_{n+k-1}}^{t_{n+k}}f(t,y(t))\,dt

the integrand f(t,y(t)) by the polynomial p(t) interpolating these values. Then we can evaluate this integral analytically and obtain the next solution value y_{n+k}. After advancing the scheme by one step, we obtain similarly y_{n+k+1}, y_{n+k+2} and so on. This procedure is illustrated here for the differential equation \dot y=a y (1-y), which represents the phenomenon of logistic growth.

Explicit formulas are obtained using Newton's formula for polynomial interpolation as follows

(4)
\begin{array}{ll} k=1:\quad & y_{n+1}=y_n+hf_n\\ k=2: \quad & y_{n+1}=y_n+h\bigl( {3\over 2}f_n- {1\over 2}f_{n-1}\bigr)\\[1mm] k=3: \quad & y_{n+1}=y_n+h\bigl( {23\over 12}f_n-{16\over 12}f_{n-1}+{5\over 12}f_{n-2} \bigr)\\[1mm] k=4: \quad & y_{n+1}=y_n+h\bigl( {55\over 24}f_n-{59\over 24}f_{n-1}+{37\over 24}f_{n-2} -{9\over 24}f_{n-3}\bigr) . \end{array}

For more details of this calculation see Hairer, Norsett & Wanner (1993), p.\,357-358.

[edit]

Implicit Adams methods

If f_{n+k}=f(t_{n+k},y_{n+k}) were known, we could use the interpolation polynomial stretched over the entire interval t_n,\ldots , t_{n+k} and thereby very much increase the precision of the result. The formulas so obtained have the form

(5)
\begin{array}{ll} k=0:\quad & y_{n+1}=y_n+hf_{n+1} \\[1mm] k=1: \quad & y_{n+1}=y_n+h\bigl( {1\over 2}f_{n+1}+ {1\over 2}f_n\bigr)\\[1mm] k=2: \quad & y_{n+1}=y_n+h\bigl( {5\over 12}f_{n+1}+{8\over 12}f_n-{1\over 12}f_{n-1} \bigr)\\[1mm] k=3: \quad & y_{n+1}=y_n+h\bigl( {9\over 24}f_{n+1}+{19\over 24}f_n -{5\over 24}f_{n-1} +{1\over 24}f_{n-2}\bigr) . \end{array}

and represent an implicit equation for computing y_{n+k}. For more details on their calculation see Hairer, Norsett & Wanner (1993), p.\,359-360. The special cases k=0 and k=1 are the implicit Euler method and the trapezoidal rule, respectively. They are actually one-step methods.

Figure 2: Predictor-corrector scheme
[edit]

Predictor-Corrector schemes

An elegant way of solving the implicit equations in (5) is to predict a value \widehat y_{n+k} by the explicit Adams method (with the same value of k) and replace in the corrector" (5) the missing value f_{n+k} by \widehat f_{n+k}=f(t_{n+k},\widehat y_{n+k}). This is the commonly used implementation of the implicit Adams methods for non-stiff problems. The complete algorithm created by this idea is illustrated as follows:

[edit]

General formulation and consistency

A linear multistep method for the numerical solution of ordinary differential equations \dot y = f(y) is given by

(6)
\sum_{i=0}^k \alpha_i y_{n+i} =h \sum_{i=0}^k \beta _i f(y_{n+i}) ,

where ....

[edit]

Stability and convergence

Figure 3: Explicit Adams method for stiff problem
[edit]

Methods for stiff problems

We obtain a bad surprise if we compute three more values of the above logistic growth problem with the explicit Adams method of order 3:

[edit]

Examples, singularly perturbed problems, discretized parabolic equations

stiff systems

[edit]

Backwad differentiation formulas (BDF)

[edit]

A-stability

[edit]

Convergence

[edit]

Methods for Hamiltonian systems

For long-time integration of conservative differential equations special integrators are usually more appropriate. The study of such integrators is the main topic of geometric numerical integration (Hairer, Lubich & Wanner 2006), and there are interesting geometric integrators among multistep methods.

[edit]

Second order Hamiltonian systems

An important class of conservative systems are second order differential equations

(7)
\ddot q = f(q), \qquad f(q)=-M^{-1}\nabla U(q)

with a constant, symmetric, positive definite mass matrix M, and a smooth potential U(q). This equation is Hamiltonian with energy \textstyle H(p,q)={1\over 2}\, p^{\mathsf T} M^{-1} p +U(q), where p=M\dot q is the momentum. The exact flow is known to be symplectic, and the energy is preserved along solutions. Typical examples are N-body systems, e.g., planetary motion in astronomy or simulations in molecular dynamics.

[edit]

Symmetric linear multistep methods

It is possible to write (7) as a first order system and to apply any multistep method considered above. Eliminating the velocity leads to a formula

(8)
\sum_{j=0}^k \alpha_{j}\,q_{n+j} = h^2 \sum_{j=0}^k \beta_{j}\, f(q_{n+j}) ,

with coefficients \alpha_{j},\beta_{j} that are different from those above. We again introduce generating polynomials

\rho (\zeta ) = \sum_{j=0}^k \alpha_{j}\, \zeta^j , \qquad \sigma (\zeta ) = \sum_{j=0}^k \beta_{j}\, \zeta^j .

A linear multistep method (8) has order r for differential equations (7) if

\rho (e^h ) - h^2 \sigma (e^h ) = O (h^{r+2}) \qquad \hbox{for}\quad h\to 0 ,

and it is stable if all zeros of \rho (\zeta ) satisfy |\zeta |\le 1 with those lying on the unit circle being at most double. Order r and stability guarantee convergence, and the global error satisfies q_{n}-q(t_{0}+nh) = O (h^r) on compact intervals nh \le T.

For the integration of Hamiltonian systems, we are mainly interested in long-time integration. For example, the solution of the harmonic oscillator \ddot q = -\omega^2 q satisfies \textstyle\frac 12 (\dot q(t) ^2 + \omega^2 q(t)^2)= const and remains bounded for all times. Can we have a similar behavior for the numerical solution, when h is fixed and n\to\infty? When applied to the harmonic oscillator, method (8) becomes a linear difference relation with characteristic equation \rho (\zeta ) - (h\omega)^2 \sigma (\zeta ) =0. It turns out that near energy conservation for the harmonic oscillator can be achieved only if all roots have modulus one. This in turn implies that the method has to be symmetric, i.e.,

(9)
\alpha_{j} = \alpha_{k-j}, \qquad \beta_{j}=\beta_{k-j}.

Lambert & Watson (1976) initiated the study of the long-time behavior of symmetric multistep methods. Quinlan & Tremaine (1990) constructed high order methods (8) and applied them successfully to long-time integrations of the outer solar system. One of their methods of order 8 with k=8 is given by

(10)
\begin{array}{rcl} \rho (\zeta ) &\!\! = \!\! & (\zeta -1)^2 (\zeta^6 + 2\zeta^5 +3\zeta^4 + 3.5 \zeta^3 + 3\zeta^2 +2\zeta +1)\\[1mm] 120960\,\sigma (\zeta ) &\!\! = \!\! & 192481 (\zeta^7 + \zeta ) + 6582 (\zeta^6 +\zeta^2) +816783 (\zeta^5 + \zeta^3) - 156812 \zeta^4 \end{array}

Notice that the method is explicit (i.e., \beta_{k}=0) and therefore its implementation with constant steps is straight-forward. Due to possible numerical resonances, symmetric multistep methods are recommended only for high accuracy computations.

[edit]

Long-time behavior

For the study of qualitative and quantitative features of the numerical solution on intervals that are far beyond the applicability of classical convergence theory, backward error analysis is the appropriate tool. It turns out that a symmetric method has to be s-stable (i.e., apart from the double root at 1, all zeros of \rho (\zeta ) are simple and of modulus one). It has been shown in Hairer & Lubich (2004) that symmetric, s-stable multistep methods (8) of order r have a long-time behavior that it very similar to that of symplectic one-step methods. In particular, their numerical solutions satisfy:

  • for a Hamiltonian system (7) the total energy \textstyle H(p,q)=\frac 12 p^{\mathsf T} M^{-1} p +U(q) is conserved up to O (h^r) on intervals of length O (h^{-r-2});
  • quadratic first integrals of the form \dot q^{\mathsf T} A q (such as the angular momentum in N-body problems) are conserved up to O (h^r) on intervals of length O (h^{-r-2}) (here, numerical approximations \dot q_{n} for the derivative are computed by symmetric finite differences);
  • for completely integrable Hamiltonian systems, action variables are in general conserved up to O (h^r) and angle variables are bounded by O (t_{n} h^r) on intervals of length O (h^{-r}).
[edit]

Implementation and software

The efficient implementation of linear multistep methods is a nontrivial task. A code has to be able to automatically select the step size (which can be done either by formulas with grid-dependent coefficients or by using a Nordsieck vector formulation) and to adjust an optimal order. Among the available software let us mention the following:

  • ODE/STEP is a widely used code for the solution of non-stiff differential equations. It is based on Adams methods and documented in the monograph of Shampine & Gordon.
  • DIVA is an Adams integrator by Fred Krogh. It can directly integrate second order differential equations.
  • LSODE is the ``Livermore Solver of Alan Hindmarsh. It solves stiff or nonstiff differential equations of first order.
  • DASSL is a differential-algebraic equation solver, based on BDF schemes. The most recent version and various modifications can be downloaded from the web page of Linda Petzold.

Most of these codes are publicly available under http://www.netlib.org/ode/

[edit]

References

  • Hairer E and Lubich C (2004) Symmetric multistep methods over long times. Numer. Math. 97:699-723.
  • Hairer E, Lubich C. and Wanner, G. (2006) Geometric numerical integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Second edition. Springer-Verlag, Berlin.
  • Hairer E, Nørsett S.P. and Wanner, G. (1993) Solving Ordinary Differential Equations I: Nonstiff Problems. Second edition. Springer-Verlag, Berlin.
  • Hairer E. and Wanner, G. (1996) Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Second edition. Springer-Verlag, Berlin. [1]
  • Lambert J.D. and Watson I.A. (1976) Symmetric multistep methods for periodic initial value problems. J. Inst. Maths. Applics. 18:189-202.
  • Quinlan G.D. and Tremaine S. (1990) Symmetric multistep methods for the numerical integration of planetary orbits. Astron. J. 100:1694-1700.
Invited by: Dr. Eugene M. Izhikevich, Editor-in-Chief of Scholarpedia, the peer-reviewed open-access encyclopedia
Retrieved from "http://www.scholarpedia.org/article/Linear_multistep_method"
For authors