Linear multistep method

A linear multistep method approximates an ODE using several previous solution and derivative values.

For an ordinary differential equation $y^{\prime }=\Phi (t,y)$, this linear multistep method uses a centered derivative across the two-step interval $[t_{i-1},t_{i+1}]$, evaluated at its midpoint $t_i$. It requires two starting values, $y_0$ and $y_1$. For a sufficiently smooth exact solution, Taylor's theorem givesThe unscaled local residual is therefore $O(h^3)$ and the method has order two. Its first characteristic polynomial is $z^2-1$, with the simple roots $1,-1$, so it satisfies the root condition for a multistep method; the second root is a parasitic amplification root. This two-step method differs from the one-step midpoint method and the implicit midpoint rule.

For the special equation $G^{\prime }=g(t)$ and $G(a)=0$, an even number $n=2r$ of steps gives $y_n=2h{\sum }_{j=1}^{r}g(t_{2j-1})$, independently of the odd-index starting value. It is the composite midpoint quadrature rule on $r$ intervals of width $2h$. An odd final index depends on $y_1$, which must be specified separately.

Take $\sigma (\zeta )=\zeta [(5+\alpha ){\zeta }^2-(4+8\alpha )\zeta +11-5\alpha ]/6$. The root condition for a multistep method holds exactly for $-1<\alpha <1$. Since ${\rho }^{\prime }(1)=\sigma (1)=2(1-\alpha )$, these and only these parameters give a convergent method. Its exponential-symbol order criterion for a multistep method has defect $-(\alpha +5)z^4/12+O(z^5)$, giving formal order three except at $\alpha =-5$, where the first defect is $z^5/10$. None of the convergent family is A-stable: for a parasitic root $r=\alpha +i\sqrt{1-{\alpha }^2}$, perturbation of $\rho (r(z))-z\sigma (r(z))=0$ gives $\mathrm{Re}[r^{\prime }(0)/r]=(\alpha -1)/12<0$. A small negative real $z$ therefore moves it outside the unit disk. At $\alpha =1$, cancellation would remove a double unit root, but the original recurrence retains that unstable parasitic behavior; it is not legitimately replaced by the reduced Backward Euler recurrence for arbitrary starting values.

Integrating a quadratic interpolant of the derivative over two steps gives this linear multistep method, which has order four. Its local truncation error is $-h^5{y}^{(5)}(t_{n+1})/90+O(h^7)$ and its zero-stability polynomial is $\rho (\zeta )={\zeta }^2-1$, with simple roots $1,-1$. The Dahlquist equivalence theorem therefore gives fourth-order convergence with sufficiently accurate starts. This attains the even two-step first Dahlquist barrier. Its undamped parasitic zero-step mode need not have favorable long-time absolute stability; finite-time convergence and stiff damping are separate requirements.

For a consistent zero-stable real linear $k$-step method with constant coefficients and only first-derivative evaluations, the attainable order is at most $k+1$ for odd $k$ and $k+2$ for even $k$. For an explicit method the bound is $k$. These are order barriers under the root condition for a multistep method, not bounds on multiderivative methods or on arbitrary variable-step formulations. The Second Dahlquist barrier adds the stronger bound of order two for irreducible A-stable methods.

Using Newton-Cotes closed quadrature for the integrated derivative yields the displayed linear multistep method. Its first characteristic polynomial is $\rho (z)=z^s-1$, whose unit-circle roots are all simple, proving zero-stability. Quadrature exactness gives order $s+1$ for odd $s$ and at least $s+2$ for even $s$. The first Dahlquist barrier gives equality, while the Dahlquist equivalence theorem gives convergence with consistent starts. These are finite-time convergence statements, not claims of A-stability.

For a real consistent zero-stable linear multistep method, transform its first characteristic polynomial by $z=(1+x)/(1-x)$ to $(1-x{)}^s\rho (z)=xP(x)$. The root condition for a multistep method puts all finite roots of $P$ in the closed left half-plane; real linear factors and conjugate quadratic factors consequently give nonnegative coefficients, after choosing $P(0)>0$. The order criterion requires a polynomial of degree at most $s$ to agree through degree $p-1$ with $P(x)x/(2\mathrm{arctanh}x)$. The latter quotient has strictly negative nonconstant even coefficients. Its first even coefficient beyond degree $s$ includes a strictly negative multiple of $P(0)$ and cannot vanish. This gives the stated bound.

Write $F(x)={\int }_0^1(1+x{)}^t(1-x{)}^{1-t}dt=x/\mathrm{arctanh}x$. Pairing $t$ and $1-t$ after differentiating twice givesThe power series of $(1-x^2{)}^{-3/2}$ has strictly positive even coefficients, and the other factor has nonnegative even coefficients. Thus every even coefficient of $F^{\prime \prime }$ is negative. Since $F$ is even with $F(0)=1$, the displayed sign pattern follows. This supplies a short positivity argument for the first Dahlquist barrier.

The exponential-symbol order criterion for a multistep method gives leading defect $-(1+5a)z^3/12$, so the family has order two except at $a=-1/5$, where it has order three. The root condition for a multistep method holds exactly for $-1\le a<1$, and the convergent A-stable members are exactly $0\le a<1$. At $a=0$ it reduces to the trapezoidal rule after removal of a zero root. At $a=1$ a repeated root one destroys convergence of a numerical method of the original recurrence despite a common polynomial factor and a formally vanishing second-order residual.

An explicit linear multistep method that integrates a polynomial extrapolation of the past derivative values. The displayed two-step member is second order in time: the linear interpolant through $f_{n-1},f_n$, integrated over the next step, gives coefficients $3/2,-1/2$. For a spatially discretized partial differential equation, the spatial error must also be included in the total local defect.

For centered-space diffusion, a Fourier mode gives $z^2-(1-3a/2)z-a/2=0$, with $a=4\mu {\sin }^2(\theta /2)$. For $a>0$ the roots have opposite signs, the positive root lies in $(0,1)$, and the negative root lies in $[-1,0)$ exactly when $P(-1)=2-2a\ge 0$. At $a=0,1$ the unit-modulus roots are simple. Hence the root condition for a multistep method holds for every Fourier frequency exactly when $\mu \le 1/4$. Under fixed $\mu$, the unscaled one-step defect is $-\Delta t^2{u}_{xxxx}/(12\mu )+O(\Delta t^3)$; dividing by the step changes its order to $O(\Delta t)$.

Canceling a common factor of the characteristic polynomials changes the allowed starting relations and can remove a parasitic mode rather than prove stability of the original recurrence. For example $\rho =(\zeta -1{)}^2$ and $\sigma =({\zeta }^2-1)/2$ fail zero-stability. On $y^{\prime }=0$, the original recurrence admits $y_n=y_0+n(y_1-y_0)$; starting perturbations $y_1-y_0=\sqrt{h}$ vanish but give divergent error at fixed positive time. The canceled one-step recurrence excludes these data, so its stability cannot be assigned to the original scheme.

This real-parameter linear multistep method has formal order two except at the degenerate value $\alpha =2$, where its local symbol vanishes to fourth degree. The zero-stability roots are $1$ and $\alpha -1$, so it is convergent exactly for $0\le \alpha <2$. Every convergent member is A-stable: the boundary-locus test for multistep A-stability gives nonnegative real part $\alpha (2-\alpha )(1-\cos \theta {)}^2/(2|\sigma (e^{i\theta }){|}^2)$, with no unit-circle denominator zero for $0<\alpha <2$. At $\alpha =0$ the two interlaced sequences obey a trapezoidal step of length $2h$. At $\alpha =4/3$ it is the BDF2 method. The endpoint two fails the root condition for a multistep method even though canceling a common polynomial factor resembles the trapezoidal rule.

A multiderivative multistep method uses higher time derivatives as well as values and first derivatives from several time levels. For an autonomous ordinary differential equation $y^{\prime }=f(y)$, the second time derivative is $y^{\prime \prime }=f^{\prime }(y)f(y)$ by the chain rule. Its local truncation error must include the higher-derivative terms. Consequently the usual second-order barrier for A-stable linear multistep methods, which only use first derivatives, does not apply.

For $y^{\prime }=f(y)$ and $g=f^{\prime }f$, inserting a smooth exact solution into the displayed recurrence gives defect $h^7{y}^{(7)}(t_n)/4725+O(h^9)$. This follows by expanding $2\sinh z-(14/15)z\cosh z-(16/15)z+(2/15)z^2\sinh z$. The terms through $z^6$ vanish and the $z^7$ coefficient is nonzero. Its zero-step polynomial is ${\zeta }^2-1$, whose roots are simple and of unit modulus. Thus it is zero-stable and has order six, provided starting values and derivative evaluations have corresponding accuracy.

For the multiderivative multistep method $y_{n+1}-\alpha h{f}_{n+1}+\beta h^2{g}_{n+1}=\gamma h{f}_n+y_{n-1}+\alpha h{f}_{n-1}+\beta h^2{g}_{n-1}$, with $g=f^{\prime }f$, centered expansion gives defect coefficients $2-2\alpha -\gamma$, $1/3-\alpha +2\beta$, and $1/60-\alpha /12+\beta /3$ at powers $h,h^3,h^5$. Their simultaneous vanishing gives the displayed parameters. The next coefficient is $1/4725$ at $h^7$, so the formal order is six. Reflection symmetry removes every even defect power. The method satisfies fifth-order conditions but cannot have exact order five.

For the symmetric two-step two-derivative formula, the closed-disk root condition for a multistep method fails off the imaginary axis. Positive real part makes the product of root moduli exceed one; negative real part fails the complex quadratic Schur criterion, as in the empty strict decay domain of the symmetric two-derivative formula proof. On $z=iy$, write $D=15-y^2-7iy=|D|e^{i\delta }$ and set $\xi =i{e}^{-i\delta }\eta$. The amplification polynomial becomes ${\eta }^2-16y\eta /|D|+1=0$. Its two roots are distinct and on the unit circle exactly when $|16y/|D||<2$, or $y^4-45y^2+225>0$. Thus the bounded-root set consists of the imaginary intervals $|y|<\sqrt{(45-15\sqrt{5})/2}$ and $|y|>\sqrt{(45+15\sqrt{5})/2}$, including $y=0$ and excluding the double-root endpoints. In particular, the set is unbounded. Unit-modulus oscillations are bounded but do not tend to zero.

On the Dahlquist test equation, write the amplification polynomial as $D{\xi }^2-16z\xi -A$, where $D=z^2-7z+15$ and $A=z^2+7z+15$. The identity $|A{|}^2-|D{|}^2=28\mathrm{Re}z(|z{|}^2+15)$ shows that both roots of a polynomial can lie inside the open unit disk only if $\mathrm{Re}z<0$. But the complex quadratic Schur criterion would then require $32|\mathrm{Re}z|(|z{|}^2+15)<28|\mathrm{Re}z|(|z{|}^2+15)$, which is impossible. Thus the strict linear stability domain is empty. This concerns robust damping of all modes, rather than convergence on fixed time intervals or the bounded-root stability set of the symmetric two-derivative formula.

On the Dahlquist test equation, the symmetric two-step two-derivative formula has amplification polynomial $F_z(w)=(z^2-7z+15)w^2-16zw-(z^2+7z+15)$. For every real $z<0$, its leading coefficient is positive and $F_z(-1)=2z<0$, whereas $F_z(w)\to +\infty$ as $w\to -\infty$. The intermediate value theorem therefore gives a real root below minus one. Thus the method is not A-stable, without any appeal to a multistep order barrier.

The endpoint expansion of this multiderivative multistep method has unscaled defect $h^4{y}^{(4)}/3+O(h^5)$ and zero-stability roots $1,1/7$, giving third-order convergence with suitably accurate starting values. On the Dahlquist test equation its amplification roots obey $(7-6z+2z^2){\zeta }^2-8\zeta +1=0$. The complex quadratic Schur criterion puts both roots in the open unit disk for $\mathrm{Re}z<0$, and only the simple root one reaches the boundary at $z=0$. Thus it is A-stable; its higher derivative places it outside the hypotheses of the Second Dahlquist barrier.

For a fixed-step multiderivative multistep method, the ordinary zero-stability root condition still controls propagation of the starting errors when all derivative evaluation maps are uniformly Lipschitz continuous on the relevant bounded region. A local defect $O(h^{p+1})$ then yields global error $O(h^p)$ over a fixed time interval, provided the starting errors are $O(h^p)$ and the implicit updates use the nearby solution branch.

For example, writing the numerical error equation as $\rho (E)e_n=h{F}_h(e_{n+s})+d_n$, with $F_h$ uniformly Lipschitz continuous, a bounded impulse response for the root condition for a multistep method givesAbsorbing the last current-step term for small $h$ and applying the discrete Gronwall inequality gives the stated order. The higher derivatives enter through $F_h$, whose bound stays uniform as $h\to 0$.

Zero-stability requires the roots of the first characteristic polynomial to lie in the closed unit disc, with unit-modulus roots simple.

The root condition is the polynomial criterion equivalent to zero-stability for a linear multistep method.

An amplification root of a multistep method is parasitic when it represents an additional numerical mode rather than the principal mode approximating the exact differential-equation flow. For a consistent method the principal root starts at one as the step tends to zero, while other roots can start elsewhere. A simple parasitic root on the unit circle is compatible with zero-stability, but its motion outside the circle for a negative test parameter can destroy absolute stability. Fixed-time convergence and all-time damping must be distinguished.

Replacing the derivative in $x^{\prime }=-Kx$ by a central finite difference gives $x_{n+1}+2Kh{x}_n-x_{n-1}=0$. Its principal root is ${\rho }_{+}=e^{-\mathrm{arsinh}(Kh)}$, while the parasitic amplification root is ${\rho }_{-}=-e^{\mathrm{arsinh}(Kh)}$. On a fixed time interval, $A_h{\rho }_{+}^n+B_h{\rho }_{-}^n$ converges uniformly to $C{e}^{-Kt}$ exactly when $A_h\to C$ and $B_h\to 0$. A fixed nonzero parasitic amplitude alternates between two incompatible limits, but a vanishing starting error can converge despite the root's modulus exceeding one. This distinguishes fixed-time numerical convergence from long-time absolute stability.

The polynomial $\rho (w)=(w-1)(w^2+w+a)$ satisfies the root condition for a multistep method exactly for $0\le a\le 1$. For $a<0$ a root is below $-1$; for $a>1$ the quadratic root product forces a root outside the unit disc. The endpoints have simple unit-modulus roots, and the repeated root at $a=1/4$ is strictly inside. A consistent method with this polynomial is convergent throughout the closed interval.

For a fixed-coefficient linear multistep method applied to a well-posed initial value problem with a continuous right-hand side uniformly Lipschitz continuous in the solution argument on the relevant region, consistency of a numerical method and zero-stability are equivalent to numerical convergence for arbitrary starting values whose errors tend to zero. The time interval is fixed while the step size tends to zero; an implicit numerical method must use its well-defined nearby solution branch for sufficiently small steps.

The quantitative statement is that a zero-stable method with consistency error $O(h^p)$ and starting errors $O(h^p)$ has global error $O(h^p)$, provided the exact solution is sufficiently smooth for that residual bound. A bounded response to step defects follows from the root condition for a multistep method, and the discrete Gronwall inequality controls the nonlinear error terms. Mere continuity of the right-hand side does not by itself supply these hypotheses or uniqueness of the exact solution.

For a linear multistep method, convergence for every set of consistent starting values requires the root condition for a multistep method. Apply the recurrence to $y^{\prime }=0$ over $N$ steps. A root $\xi$ with $|\xi |>1$ gives the error mode $e_n=|\xi {|}^{-N}{\xi }^n$: its fixed-index starting errors vanish as $N\to \infty$, but $|e_N|=1$. A unit-modulus root of multiplicity $m\ge 2$ instead gives $e_n=N^{-(m-1)}{n}^{m-1}{\xi }^n$, with the same failure. Repeated roots inside the unit disk are harmless because their polynomial growth is dominated by exponential decay.

For ${\sum }_{j=0}^s{\alpha }_j{y}_{n+j}=h{\sum }_{j=0}^s{\beta }_j{f}_{n+j}$, the characteristic polynomials are $\rho (\zeta )={\sum }_j{\alpha }_j{\zeta }^j$ and $\sigma (\zeta )={\sum }_j{\beta }_j{\zeta }^j$. Applied to $y^{\prime }=\lambda y$, the amplification roots solve $\rho (\zeta )-h\lambda \sigma (\zeta )=0$.

Applying a linear multistep method to the Dahlquist test equation with $z=h\lambda$ produces a constant-coefficient recurrence whose modes satisfy $P_z(w)=0$. Absolute stability requires the resulting roots of a polynomial to lie in the unit disk, with unit-modulus roots simple. Thus multistep stability depends on every recurrence mode, rather than a single one-step stability function.

An amplification root is a root of a polynomial governing a geometric mode $y_n={\xi }^n$ of a numerical recurrence on the Dahlquist test equation. Its modulus measures damping or growth per step. A repeated root produces additional modes such as $n{\xi }^n$; therefore a unit-circle root can be permitted by a bounded-root condition only when it is simple. The principal root approximates $e^z$ near zero, while additional roots can describe parasitic amplification roots.

If the derivative characteristic polynomial $\sigma$ of a linear multistep method has a simple root outside the unit disk and keeps its degree in the limit, then its absolute stability region is bounded. Otherwise a stable sequence with $|z|\to \infty$ would give amplification polynomials $\sigma -\rho /z$ tending coefficientwise to $\sigma$. A root near the exterior root would persist by continuity, contradicting the root condition for a multistep method. This argument excludes stability at every sufficiently large complex $z$, not only along the negative real axis.

A linear multistep method has order $p$ exactly whenfor $0\le q\le p$, with the right side interpreted as zero for $q=0$, and the identity first fails at $q=p+1$.

The Taylor expansion of the exact-solution residual for a linear multistep method has the same coefficient sequence asVanishing through degree $p$ is therefore equivalent to order at least $p$. Exact order $p$ additionally requires a nonzero coefficient at degree $p+1$. This symbol conveniently collects all polynomial consistency and order conditions.

An Adams–Moulton method is an implicit linear multistep method obtained by integrating an interpolation polynomial for the vector field that includes the new time level.

An irreducible A-stable linear multistep method has order at most two.

This is an order restriction for an A-stable linear multistep method.