Linear multistep method

• Dr. Ernst Hairer, Department of Mathematics, University of Geneva, Switzerland

• Dr. Gerhard Wanner, Department of Mathematics, University of Geneva, Switzerland

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 for ordinary differential equations, and particular methods are well suited for solving non-stiff and stiff equations as well as Hamiltonian systems over long time intervals. In contrast to one-step methods, multistep methods increase efficiency by using the information from several previous solution approximations.

Methods for nonstiff problems

Consider an initial value problem

<math eq:diffeq>

\dot y=f(t,y), \quad y(t_0)=y_0 , </math> where the solution \(y(t)\) can be scalar- or vector-valued. Let \(t_0,t_1,t_2,\ldots\) be a sequence of grid points, which for simplicity is supposed to be equidistant with step size \(h=t_{n+1}-t_n\). A numerical method is an algorithm that yields approximations \(y_n\) to the solution \(y(t_n)\) at the grid points.

Figure 1: Explicit Adams methods

Explicit Adams methods

These methods are based on the following idea: suppose that \(k\) values \( y_{n-k+1},\,y_{n-k+2},\ldots, y_{n}\) approximating \(y(t_{n-k+1}),\, y(t_{n-k+2}),\ldots,y(t_{n}) \) are known. Compute the derivatives

<math leq:slopes>

f_{n+j}=f(t_{n+j},y_{n+j}), \quad j=-k+1,\ldots ,0 </math> and replace in the integrated form of (<ref>eq:diffeq</ref>)

<math leq:diffeqint>

y(t_{n+1})=y(t_{n})+\int _{t_{n}}^{t_{n+1}}f(t,y(t))\,dt </math> the integrand \(f(t,y(t))\) by the polynomial \(p(t)\) interpolating the values (<ref>leq:slopes</ref>). Then evaluate the integral analytically and obtain the next approximation to the solution, \(y_{n+1}\). After advancing the scheme by one step, this procedure can be repeated to obtain \(y_{n+2},\, y_{n+3}\), and so on. See Figure 1 for an illustration using the logistic growth equation \(\dot y=a y (1-y)\).

Newton's formula for polynomial interpolation leads to the following equations (for more details see Hairer, Nørsett & Wanner (1993), p. 357-358):

<math leq:adamsex>

\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)\UNIQ3dceff6c77138a71-MathJax-1-QINU The accuracy of methods having the same order can be distinguished by their ''error constant'' UNIQ3dceff6c77138a71-MathJax-2-QINU ====Stability and convergence==== The basic theory for linear multistep methods has been developped by Dahlquist (1959) and is exposed in the classical book of Henrici (1962). The main difference to one-step methods is that consistency alone does not guarantee convergence. A linear multistep method is called ''stable'', if for the differential equation UNIQ3dceff6c77138a71-MathJax-19-QINU all solutions of the difference equation (<ref>eq:lmm</ref>) remain bounded. This is equivalent to the property that all zeros of UNIQ3dceff6c77138a71-MathJax-20-QINU satisfy UNIQ3dceff6c77138a71-MathJax-21-QINU with those lying on the unit circle being simple. The highest attainable order of linear UNIQ3dceff6c77138a71-MathJax-22-QINU-step methods is UNIQ3dceff6c77138a71-MathJax-23-QINU. However, these methods are of no practical use because of the following result. '''First Dahlquist barrier.''' The order UNIQ3dceff6c77138a71-MathJax-24-QINU of a stable linear UNIQ3dceff6c77138a71-MathJax-25-QINU-step method satisfies *UNIQ3dceff6c77138a71-MathJax-26-QINUUNIQ3dceff6c77138a71-MathJax-27-QINU if UNIQ3dceff6c77138a71-MathJax-28-QINU is even, *UNIQ3dceff6c77138a71-MathJax-29-QINUUNIQ3dceff6c77138a71-MathJax-30-QINU if UNIQ3dceff6c77138a71-MathJax-31-QINU is odd, *UNIQ3dceff6c77138a71-MathJax-32-QINUUNIQ3dceff6c77138a71-MathJax-33-QINUUNIQ3dceff6c77138a71-MathJax-34-QINUUNIQ3dceff6c77138a71-MathJax-35-QINU if UNIQ3dceff6c77138a71-MathJax-36-QINU, in particular if the method is explicit. '''Convergence theorem.''' If a linear multistep method (<ref>eq:lmm</ref>) is stable and of order UNIQ3dceff6c77138a71-MathJax-37-QINU, then it is convergent of order UNIQ3dceff6c77138a71-MathJax-38-QINU. Convergence of order UNIQ3dceff6c77138a71-MathJax-39-QINU means that for sufficiently accurate starting approximations UNIQ3dceff6c77138a71-MathJax-40-QINU the global error satisfies UNIQ3dceff6c77138a71-MathJax-41-QINU on compact intervals UNIQ3dceff6c77138a71-MathJax-42-QINU. The constant symbolized by the UNIQ3dceff6c77138a71-MathJax-43-QINU notation is bounded by UNIQ3dceff6c77138a71-MathJax-44-QINU, where UNIQ3dceff6c77138a71-MathJax-45-QINU is a Lipschitz constant of the vector field UNIQ3dceff6c77138a71-MathJax-46-QINU. ==Methods for stiff problems== [[Image:sta1.gif|right|thumb|280px|Figure |Explicit Adams method for stiff problem]] Continuing the computation of the above logistic growth equation UNIQ3dceff6c77138a71-MathJax-47-QINU, using UNIQ3dceff6c77138a71-MathJax-48-QINU, with the explicit Adams process of order 3 leads to a bad surprise. Instead of approaching the steady-state solution UNIQ3dceff6c77138a71-MathJax-49-QINU, the numerical solution exhibits a violent instability phenomenon (see Figure 3). [[Stiff systems|Stiff differential equations]] are characterized by the fact that the use of explicit methods requires very small step sizes, and certain implicit methods (like the implicit Euler method) perform much better. Stiff differential equations arise frequently as singularly perturbed problems in chemical reaction systems and in electrical circuitry, and as space discretizations of parabolic partial differential equations. ====Stiff stability analysis==== The foregoing stability analysis for nonstiff problems, which is uniquely based on the polynomial UNIQ3dceff6c77138a71-MathJax-50-QINU, is not sufficient for stiff differential equations. The stiff stability analysis, initiated by Dahlquist (1963), is done by linearization in the neighborhood of a stationary solution, translation and, in the case of systems, diagonalization, which leads to :<math eq:testeq> \dot y=\lambda y , </math> where UNIQ3dceff6c77138a71-MathJax-51-QINU stands for an eigenvalue of the Jacobian of UNIQ3dceff6c77138a71-MathJax-52-QINU and can be a complex number. In the above case, we have UNIQ3dceff6c77138a71-MathJax-53-QINU. Method (<ref>eq:lmm</ref>) applied to this equation gives UNIQ3dceff6c77138a71-MathJax-54-QINU, whose stability depends on the roots of [[Image:ps-adams-rootlocus.jpeg|thumb|500px|right|Figure ??|Root-locus curve for explicit Adams method ]] :<math eq:lintesteq> \rho(\zeta) -\mu\sigma(\zeta)=0 ,\qquad \mu=h\lambda . </math> The ''stability domain'' UNIQ3dceff6c77138a71-MathJax-55-QINU is the set of complex numbers UNIQ3dceff6c77138a71-MathJax-56-QINU, for which all roots of (<ref>eq:lintesteq</ref>) satisfy UNIQ3dceff6c77138a71-MathJax-57-QINU with those lying on te unit circle being simple. For the explicit Adams method with UNIQ3dceff6c77138a71-MathJax-58-QINU, equation (<ref>eq:lintesteq</ref>) becomes :<math eq-left> { \zeta ^3-\zeta ^2-\mu\Bigl( \frac{23}{12}\zeta ^2-\frac{16}{12}\zeta+\frac{5}{12}\Bigr) =0}, </math> or :<math eq-right> \mu=\frac{\zeta ^3-\zeta ^2}{\frac{23}{12}\zeta ^2-\frac{16}{12}\zeta +\frac{5}{12}} . </math> [[Image:ps-adams-pc-rootlocus.jpeg|thumb|160px|right|Figure ??|Root-locus curve for predictor-corrector ]] An elegant way of avoiding the high degree equation in (<ref>eq-left</ref>) is to insert the limiting values UNIQ3dceff6c77138a71-MathJax-59-QINU, for UNIQ3dceff6c77138a71-MathJax-60-QINU real, into (<ref>eq-right</ref>). This leads to the ''root locus curve'' for UNIQ3dceff6c77138a71-MathJax-61-QINU, parts of which must form the boundary of the stability domain UNIQ3dceff6c77138a71-MathJax-62-QINU. These domains are drawn for four explicit Adams methods in Figure 4 and shrink rapidly with increasing UNIQ3dceff6c77138a71-MathJax-63-QINU. In the example from the beginning of this section, the value of UNIQ3dceff6c77138a71-MathJax-64-QINU was UNIQ3dceff6c77138a71-MathJax-65-QINU, far beyond the stability bound UNIQ3dceff6c77138a71-MathJax-66-QINU. Implicit Adams methods for UNIQ3dceff6c77138a71-MathJax-67-QINU have finite stability domains as well, but they are considerably larger than those for explicit Adams methods. However, if they are solved by a predictor-corrector approach, a good deal of the stability is lost, as seen in Figure 5 for the case UNIQ3dceff6c77138a71-MathJax-68-QINU. ====Backward differentiation formulas (BDF)==== This class of methods has been discovered, together with the phenomenon of stiffness, by Curtiss & Hirschfelder (1952) in calculations for chemistry, and their extreme importance for stiff problems is recognized since the work of Gear (1971). In contrast to the Adams methods, where polynomial interpolations were used for the slopes UNIQ3dceff6c77138a71-MathJax-69-QINU, we here interpolate the UNIQ3dceff6c77138a71-MathJax-70-QINU-values UNIQ3dceff6c77138a71-MathJax-71-QINU and require for the interpolating polynomial UNIQ3dceff6c77138a71-MathJax-72-QINU the ''collocation condition'' UNIQ3dceff6c77138a71-MathJax-73-QINU. The differentiation of Newton's interpolation formula then leads to the methods :<math eq:bdfgen> \sum_{j=1} ^k {1 \over j}\nabla^j y_{n+1} = hf_{n+1} , </math> (UNIQ3dceff6c77138a71-MathJax-74-QINU denotes the backward difference operator) or, when developed, :<math eq:bdfdev> \begin{array}{ll} k =1:\quad & y_{n+1} - y_n = hf_{n+1} \UNIQ3dceff6c77138a71-MathJax-3-QINU ====Convergence==== The interest of applying BDF methods (or other A-stable or AUNIQ3dceff6c77138a71-MathJax-75-QINU-stable multistep methods) relies in the fact that they can be applied to stiff differential equations with step sizes that are much larger than the inverse of a Lipschitz constant of UNIQ3dceff6c77138a71-MathJax-76-QINU. In this situation, classical convergence theory (for non-stiff problems) does not provide any information. It is nevertheless possible to prove estimates UNIQ3dceff6c77138a71-MathJax-77-QINU of the global error (for sufficiently accurate starting approximations and for AUNIQ3dceff6c77138a71-MathJax-78-QINU-stable methods of order UNIQ3dceff6c77138a71-MathJax-79-QINU) in the following situations: *''singularly perturbed problems'' UNIQ3dceff6c77138a71-MathJax-80-QINU: if the eigenvalues of UNIQ3dceff6c77138a71-MathJax-81-QINU are in the negative half plane, then the constant symbolized by UNIQ3dceff6c77138a71-MathJax-82-QINU is independent of UNIQ3dceff6c77138a71-MathJax-83-QINU. *''discretized parabolic equations'' and stiff differential equations satisfying a one-sided Lipschitz constant. Precise statements and proofs can be found, for example, in Sections V.7, V.8, and VI.2 of the monograph Hairer & Wanner (1996). ==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. ====Second order Hamiltonian systems==== An important class of conservative systems are second order differential equations :<math ode2> \ddot q = f(q), \qquad f(q)=-M^{-1}\nabla U(q) </math> with a constant, symmetric, positive definite mass matrix UNIQ3dceff6c77138a71-MathJax-84-QINU, and a smooth potential UNIQ3dceff6c77138a71-MathJax-85-QINU. This equation is Hamiltonian with energy UNIQ3dceff6c77138a71-MathJax-86-QINU, where UNIQ3dceff6c77138a71-MathJax-87-QINU is the momentum. The exact flow is known to be symplectic, and the energy is preserved along solutions. Typical examples are UNIQ3dceff6c77138a71-MathJax-88-QINU-body systems, e.g., planetary motion in astronomy or simulations in molecular dynamics. ====Symmetric linear multistep methods==== It is possible to write (<ref>ode2</ref>) as a first order system and to apply any multistep method considered above. Eliminating the velocity leads to a formula :<math lmm2> \sum_{j=0}^k \alpha_{j}\,q_{n+j} = h^2 \sum_{j=0}^k \beta_{j}\, f(q_{n+j}) , </math> with coefficients UNIQ3dceff6c77138a71-MathJax-89-QINU that are different from those above. We again introduce generating polynomials UNIQ3dceff6c77138a71-MathJax-4-QINU A linear multistep method (<ref>lmm2</ref>) has ''order'' UNIQ3dceff6c77138a71-MathJax-90-QINU for differential equations (<ref>ode2</ref>) if UNIQ3dceff6c77138a71-MathJax-5-QINU and it is ''stable'' if all zeros of UNIQ3dceff6c77138a71-MathJax-91-QINU satisfy UNIQ3dceff6c77138a71-MathJax-92-QINU with those lying on the unit circle being at most double. Order UNIQ3dceff6c77138a71-MathJax-93-QINU and stability guarantee ''convergence'', and the global error satisfies UNIQ3dceff6c77138a71-MathJax-94-QINU on compact intervals UNIQ3dceff6c77138a71-MathJax-95-QINU. For the integration of Hamiltonian systems, we are mainly interested in long-time integration. For example, the solution of the harmonic oscillator UNIQ3dceff6c77138a71-MathJax-96-QINU satisfies UNIQ3dceff6c77138a71-MathJax-97-QINU and remains bounded for all times. Can we have a similar behavior for the numerical solution, when UNIQ3dceff6c77138a71-MathJax-98-QINU is fixed and UNIQ3dceff6c77138a71-MathJax-99-QINU? When applied to the harmonic oscillator, method (<ref>lmm2</ref>) becomes a linear difference relation with characteristic equation UNIQ3dceff6c77138a71-MathJax-100-QINU. 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., :<math symm> \alpha_{j} = \alpha_{k-j}, \qquad \beta_{j}=\beta_{k-j}. </math> Lambert & Watson (1976) initiated the study of the long-time behavior of symmetric multistep methods. Quinlan & Tremaine (1990) constructed high order methods (<ref>lmm2</ref>) and applied them successfully to long-time integrations of the outer solar system. One of their methods of order UNIQ3dceff6c77138a71-MathJax-101-QINU with UNIQ3dceff6c77138a71-MathJax-102-QINU is given by :<math quinlan> \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) ,\\[2mm] 120960\,\sigma (\zeta ) &\!\! = \!\! & 192481 (\zeta^7 + \zeta ) + 6582 (\zeta^6 +\zeta^2) +816783 (\zeta^5 + \zeta^3) - 156812 \zeta^4 . \end{array} </math> 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 computations with high accuracy.

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 (<ref>lmm2</ref>) 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 (<ref>ode2</ref>) 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})\).

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 (most of it can be downloaded from http://www.netlib.org/ode/) 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 (see Brenan, Campbell & Petzold 1989). The most recent version and various modifications can be downloaded from the web page of Linda Petzold at http://www.cs.ucsb.edu/~cse/software.html
• MEBDF is an implementation of modified extended BDF methods (Cash & Considine 1992). Various modifications can be obtained from http://www2.imperial.ac.uk/~jcash/

References

• Brenan K.E., Campbell S.L. and Petzold L.R. (1989 and 1996) Numerical solution of initial-value problems in differential-algebraic equations. North-Holland, New York, and SIAM, Philadelphia.

• Cash J.R. and Considine S. (1992) An MEBDF code for stiff initial value problems. ACM Trans. Math. Software 18:142-158.

• Curtiss C.F. and Hirschfelder J.O. (1952) Integration of stiff equations. Proc. Nat. Acad. Sci. 38:235-243.

• Dahlquist G. (1956) Convergence and stability in the numerical integration of ordinary differential equations. Math. Scand. 4:33-53.

• Dahlquist G. (1959) Stability and error bounds in the numerical integration of ordinary differential equations. Trans. of the Royal Inst. of Techn., Nr. 130, Stockholm, Sweden.

• Dahlquist G. (1963) A special stability problem for linear multistep methods. BIT 3:27-43.

• Gear C.W. (1971) Numerical initial value problems in ordinary differential equations. Prentice Hall.

• 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. [1]

• Hairer E, Nørsett S.P. and Wanner, G. (1993) Solving Ordinary Differential Equations I: Nonstiff Problems. Second edition. Springer-Verlag, Berlin. [2]

• Hairer E. and Wanner, G. (1996) Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Second edition. Springer-Verlag, Berlin. [3]

• Henrici P. (1962) Discrete Variable Methods in Ordinary Differential Equations. Volume 1, John Wiley & Sons, New York.

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

Sponsored by: Eugene M. Izhikevich, Editor-in-Chief of Scholarpedia, the peer-reviewed open-access encyclopedia

Sponsored by: Prof. Barbara Zubik-Kowal, Department of Mathematics, Boise State University, Idaho, USA

Reviewed by: Anonymous

Accepted on: 2010-03-28 10:20:28 GMT