Linear multistep methods and Dahlquist equivalence

Every denominator $t_{n-j}-t_{n-\mathrm{\ell }}$ in $L_j^P$ is nonzero because the nodes are distinct. Direct substitution gives $L_j^P(t_{n-r})=1$ if $j=r$ and $0$ otherwise; for $k=1$ the empty product is $1$, so this also covers the one-node predictor.

Therefore $P(s)={\sum }_{j=0}^{k-1}{L}_j^P(s)f(t_{n-j},Y_{n-j})$ is the unique polynomial of degree at most $k-1$ interpolating the known derivative samples: it has the prescribed value at each node, and the difference of two such polynomials would have $k$ distinct roots while having degree at most $k-1$.

Integrating the identity defining $P$ from $t_n$ to $t_{n+1}$ gives ${\int }_{t_n}^{t_{n+1}}P(s)ds={\sum }_{j=0}^{k-1}{b}_j^{n}f(t_{n-j},Y_{n-j})$. The Adams relation $Y_{n+1}^P=Y_n+\int P$ yields exactly the predictor coefficients stated.

The same cardinal-value argument on the $k+1$ corrector nodes gives the unique degree-at-most-$k$ interpolant; integrating it gives ${\sum }_{j=0}^k{\widehat{b}}_j^{n}f(t_{n+1-j},Y_{n+1-j})$, with the $j=0$ term evaluated at the new endpoint for the fully implicit formula.

Fix $t_n$ and expand each $y(t_n+jh)$ through degree $p$ and each $y^{\mathrm{\prime }}(t_n+jh)$ through degree $p-1$ by [F2], on the interval containing the finitely many nodes. Bounded derivatives on the compact interval make the remainders uniform over $n$ and $j$.

Substituting those expansions into [F1] and collecting the coefficient of $y^{(m)}(t_n)h^m\mathrm{/}m!$ gives ${\tau }_n(h;y)={\sum }_{m=0}^p{\mu }_m\frac{h^m}{m!}{y}^{(m)}(t_n)+O(h^{p+1})$. The derivative term contributes $m{\sum }_j{\beta }_j{j}^{m-1}$ when $m\ge 1$ and contributes nothing for $m=0$.

If ${\mu }_0,\dots ,{\mu }_p$ vanish, step 2.1 gives ${\tau }_n=O(h^{p+1})$ uniformly, so the method has order at least $p$.

Conversely, suppose the method has order at least $p$. Induct on $m=0,\dots ,p$, testing the residual on the polynomial solution $y(t)=t^m$ of $y^{\mathrm{\prime }}=m{t}^{m-1}$ (with $y=1$, $y^{\mathrm{\prime }}=0$ when $m=0$). Lower moments vanish by the induction hypothesis, so the residual is exactly ${\mu }_m{h}^m$; since $m\le p$, it can be $O(h^{p+1})$ only if ${\mu }_m=0$. Thus every required moment vanishes.

Finally, ${\mu }_0=\rho (1)$ and ${\mu }_1={\rho }^{\mathrm{\prime }}(1)-\sigma (1)$ by the polynomial definitions. Their vanishing is exactly consistency.

The companion matrix has characteristic polynomial $\rho \mathrm{/}{\alpha }_k$ by [F2]. By [F4] and [F5], write $C=SJ{S}^{-1}$ over $\mathbb{C}$, where $J$ is a block diagonal sum of Jordan blocks. Each block $J_m(\lambda )$ contributes $(z-\lambda {)}^m$ to the characteristic polynomial by [F6] and [F7]. Thus the algebraic multiplicity of each root of $\rho$ is the sum of the sizes of its Jordan blocks.

If $\rho$ has a root $\lambda$ with $\mathrm{\mid }\lambda \mathrm{\mid }>1$, the complex sequence $u_n={\lambda }^n$ solves the homogeneous recurrence. If $\lambda$ is real, this is an unbounded real solution. If $\lambda$ is nonreal, its real and imaginary parts solve the real recurrence because the coefficients are real, and $\max (\mathrm{\mid }\mathrm{Re}{u}_n\mathrm{\mid },\mathrm{\mid }\mathrm{Im}{u}_n\mathrm{\mid })\ge \mathrm{\mid }{u}_n\mathrm{\mid }\mathrm{/}\sqrt{2}$ for every $n$; hence at least one of those real solutions is unbounded along a subsequence. By [F3], the corresponding real solution has an initial state whose companion iterates are unbounded.

If a root $\lambda$ of modulus one has multiplicity $m\ge 2$, write $\rho (z)=(z-\lambda {)}^{m}q(z)$. For the forward shift $E{u}_n=u_{n+1}$, $(E-\lambda )(n^s{\lambda }^n)={\lambda }^{n+1}((n+1{)}^s-n^s)$, whose polynomial factor has degree $s-1$. Iterating $m$ times gives $(E-\lambda {)}^{m}u=0$ for $u_n=n^{m-1}{\lambda }^n$; since $q(E)$ commutes with $(E-\lambda {)}^m$, this sequence satisfies $\rho (E)u=0$, yet $\mathrm{\mid }{u}_n\mathrm{\mid }=n^{m-1}$ is unbounded. If coefficients are real and $\lambda$ is nonreal, the real and imaginary parts are real recurrence solutions and at least one is unbounded along a subsequence.

Suppose the root condition holds. If $\mathrm{\mid }\lambda \mathrm{\mid }=1$, its algebraic multiplicity is one by [F1], so its only Jordan block has size one.

For a block with $0<\mathrm{\mid }\lambda \mathrm{\mid }<1$, the binomial identity gives $J_m(\lambda {)}^n={\sum }_{j=0}^{m-1}\left(\genfrac{}{}{0px}{}{n}{j}\right){\lambda }^{n-j}{N}^j$. For each fixed $j$, the ratio of successive terms $n^j\mathrm{\mid }\lambda {\mathrm{\mid }}^{n-j}$ tends to $\mathrm{\mid }\lambda \mathrm{\mid }<1$, so this scalar sequence is bounded; if $\lambda =0$, the block powers vanish for $n\ge m$. A unit-modulus block has size one by step 2.1 and therefore has powers of modulus one. Finitely many blocks give ${\sup }_n\mathrm{\parallel }{J}^n{\mathrm{\parallel }}_{\mathrm{\infty }}<\mathrm{\infty }$.

For any matrices $A,B$, expanding entries and applying the triangle inequality gives $\mathrm{\parallel }AB{\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\parallel }B{\mathrm{\parallel }}_{\mathrm{\infty }}$. Hence $\mathrm{\parallel }{C}^n{\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }S{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\parallel }{J}^n{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\parallel }{S}^{-1}{\mathrm{\parallel }}_{\mathrm{\infty }}$, uniformly in $n$.

Under the root condition, steps 2.1, 3.1, and 4.1 bound the powers of $C$. If it fails, steps 1.2–1.3 give an unbounded homogeneous mode, and [F3] turns that mode into an initial state with unbounded companion iterates. Thus the powers are bounded exactly under the root condition, and [F1] gives the equivalent zero-stability statement.

Set $s=(t-t_{n+1})\mathrm{/}{h}_n$. The interpolation nodes become $0,-1,-(1+1\mathrm{/}{r}_n)$. Their Lagrange basis derivatives at $s=0$ are respectively $(1+2r_n)\mathrm{/}(1+r_n)$, $-(1+r_n)$, and $r_n^2\mathrm{/}(1+r_n)$, as follows by differentiating the three quadratic basis products. Since $ds\mathrm{/}dt=1\mathrm{/}{h}_n$, the derivative formula in the statement follows.

Apply [F1] separately to each component at $t_{n+1}$ through degree two. In the weighted combination $S=(1+2r_n)y(t_{n+1})-(1+r_n{)}^{2}y(t_n)+r_n^{2}y(t_{n-1})$, the constant coefficient is zero; the coefficient of $y^{\mathrm{\prime }}$ is $(1+r_n)h_n$; and the coefficient of $y^{\mathrm{\prime }\mathrm{\prime }}$ is $-\frac{1}{2}(1+r_n{)}^2{h}_n^2+\frac{1}{2}{r}_n^2(h_n+h_{n-1}{)}^2=0$.

Let $M$ bound the absolute third derivatives of every component. By [F1], the max norm of the remaining weighted remainders is at most $(M\mathrm{/}6)((1+r_n{)}^2{h}_n^3+r_n^2(h_n+h_{n-1}{)}^3)$. On the fixed compact ratio interval, $h_{n-1}=h_n\mathrm{/}{r}_n$ makes this at most $C{h}_n^3$ with a uniform finite $C$. The canceled terms in step 1.2 leave exactly these remainders after subtracting the exact derivative term, so $\mathrm{\parallel }{\tau }_n{\mathrm{\parallel }}_{\mathrm{\infty }}\le C{h}_n^3$.

Assume consistency and the root condition. Let $C$ be the companion matrix of $\rho \mathrm{/}{\alpha }_k$. By [F3], ${\sup }_{m\ge 0}\mathrm{\parallel }{C}^m{\mathrm{\parallel }}_{\mathrm{\infty }}\le M$. For a $d$-vector IVP, lift $C$ to the block companion $\mathcal{C}=C\otimes I_d$; its induced max row-sum norm is the same as that of $C$. Define $\mathrm{\parallel }U{\mathrm{\parallel }}_{\ast }:={\sup }_{m\ge 0}\mathrm{\parallel }{\mathcal{C}}^{m}U{\mathrm{\parallel }}_{\mathrm{\infty }}$. Then $\mathrm{\parallel }U{\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }U{\mathrm{\parallel }}_{\ast }\le M\mathrm{\parallel }U{\mathrm{\parallel }}_{\mathrm{\infty }}$ and $\mathrm{\parallel }\mathcal{C}U{\mathrm{\parallel }}_{\ast }\le \mathrm{\parallel }U{\mathrm{\parallel }}_{\ast }$.

Apply one step to the exact history. For an implicit method, [A2] supplies its selected root in $K$; subtracting the exact residual and using [A1] bounds its displacement from $y(t_{n+k})$ by $\mathrm{\parallel }{\tau }_n{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{/}(\mathrm{\mid }{\alpha }_k\mathrm{\mid }-h\mathrm{\mid }{\beta }_k\mathrm{\mid }L)$. The denominator is at least $\mathrm{\mid }{\alpha }_k\mathrm{\mid }\mathrm{/}2$ for small $h$. For an explicit method, ${\beta }_k=0$ and the unique output differs from $y(t_{n+k})$ by $-{\tau }_n\mathrm{/}{\alpha }_k$, without evaluating $f$ at that output. By [F1] and [F2], ${\tau }_n=O(h^{p+1})$ for order at least $p$ (take $p=1$ for consistency). In either case the output for the exact history has displacement $O(h^{p+1})$ and lies in $K$ for small $h$, since the compact exact trajectory lies in its interior.

If the root condition fails, choose a violating root $\lambda$ of $\rho$. For $\mathrm{\mid }\lambda \mathrm{\mid }>1$, the homogeneous sequence $u_n={\lambda }^n$ gives, on $h=T\mathrm{/}N$, starting perturbations $h^{1\mathrm{/}2}(u_0,\dots ,u_{k-1})$ that tend to zero, while the terminal complex magnitude is $h^{1\mathrm{/}2}\mathrm{\mid }\lambda {\mathrm{\mid }}^N\to \mathrm{\infty }$. For a repeated unit root of multiplicity $m\ge 2$, the homogeneous mode $u_n=n^{m-1}{\lambda }^n$ has bounded initial entries and terminal complex magnitude $N^{m-1}$; scaling its initial entries by $h^{1\mathrm{/}2}$ makes them tend to zero but gives terminal complex magnitude $h^{1\mathrm{/}2}{N}^{m-1}=T^{m-1}{h}^{1\mathrm{/}2-(m-1)}\to \mathrm{\infty }$. For a repeated root the claimed mode satisfies the recurrence because $(E-\lambda {)}^m(n^{m-1}{\lambda }^n)=0$, where $E{u}_n=u_{n+1}$: each application reduces the polynomial degree by one, and $\rho (E)$ has this factor. For each $N$, choose $v_n=\mathrm{Re}{u}_n$ if $\mathrm{\mid }\mathrm{Re}{u}_N\mathrm{\mid }\ge \mathrm{\mid }\mathrm{Im}{u}_N\mathrm{\mid }$, and otherwise choose $v_n=\mathrm{Im}{u}_n$. Both projections solve the real scalar recurrence, and the chosen one has $\mathrm{\mid }{v}_N\mathrm{\mid }\ge \mathrm{\mid }{u}_N\mathrm{\mid }\mathrm{/}\sqrt{2}$. Its first $k$ entries remain bounded uniformly in $N$. The real $d$-dimensional run $Y_n=h^{1\mathrm{/}2}(v_n,0,\dots ,0)$ therefore has starting errors tending to zero and terminal error diverging for the exact zero solution of $y^{\mathrm{\prime }}=0$, for every fixed $d\ge 1$. The projection may depend on $N$, as the definition of convergence permits an arbitrary sequence of starting values.

If $\rho (1)\mathrm{\ne }0$, apply the recurrence to $y^{\mathrm{\prime }}=0$, $y\equiv 1$, with exact starting values. At the first generated value $Y_k-1=-\rho (1)\mathrm{/}{\alpha }_k\mathrm{\ne }0$, so the maximum grid error does not tend to zero. Thus convergence requires $\rho (1)=0$.

Choose $\delta >0$ so that every max-norm ball of radius $\delta$ about the exact trajectory lies in $K$. For an explicit method put $B=\mathrm{\mid }{\alpha }_k{\mathrm{\mid }}^{-1}{\sum }_{j<k}\mathrm{\mid }{\alpha }_j\mathrm{\mid }$ and $D=\mathrm{\mid }{\alpha }_k{\mathrm{\mid }}^{-1}L{\sum }_{j<k}\mathrm{\mid }{\beta }_j\mathrm{\mid }$. If a history has error at most ${\eta }_0<\delta$, subtraction of its explicit formula from the formula for the exact history bounds the output difference by $(B+hD){\eta }_0$. Choose ${\eta }_0>0$ so that $(B+1){\eta }_0<\delta \mathrm{/}2$ and take $h$ small; step 1.2 then puts every such output in $K$. For an implicit method take ${\eta }_0\le \eta$ small enough that all history values lie in $K$; [A2] puts its selected outputs in $K$. Write the history as $U_n=(Y_n,\dots ,Y_{n+k-1}{)}^T$. Its update is $U_{n+1}=\mathcal{C}{U}_n+h{e}_k{R}_n(U_n,(U_{n+1}{)}_k)$, where $e_k$ inserts a vector in the last state slot and $R_n$ is the derivative sum divided by ${\alpha }_k$. For any two histories in this neighborhood, [A1], step 1.1 and the finite coefficients give $\mathrm{\parallel }{U}_{n+1}-V_{n+1}{\mathrm{\parallel }}_{\ast }\le \mathrm{\parallel }{U}_n-V_n{\mathrm{\parallel }}_{\ast }+h{K}_0(\mathrm{\parallel }{U}_n-V_n{\mathrm{\parallel }}_{\ast }+\mathrm{\parallel }{U}_{n+1}-V_{n+1}{\mathrm{\parallel }}_{\ast })$. For explicit methods the endpoint term is absent. Absorbing it when $h{K}_0\le 1\mathrm{/}2$ gives $\mathrm{\parallel }{U}_{n+1}-V_{n+1}{\mathrm{\parallel }}_{\ast }\le (1+Ch)\mathrm{\parallel }{U}_n-V_n{\mathrm{\parallel }}_{\ast }$.

Suppose $\rho (1)=0$ but ${\rho }^{\mathrm{\prime }}(1)\mathrm{\ne }\sigma (1)$. If the root condition already fails, step 1.3 applies; otherwise the root at $1$ is simple and ${\rho }^{\mathrm{\prime }}(1)\mathrm{\ne }0$. Put $\delta ={\rho }^{\mathrm{\prime }}(1)-\sigma (1)$ and test $y^{\mathrm{\prime }}=1$, $y(t)=t$, with exact starts. The numerical error satisfies ${\sum }_j{\alpha }_j{e}_{n+j}=-h\delta$. The sequence $p_n=-\delta nh\mathrm{/}{\rho }^{\mathrm{\prime }}(1)$ is a particular solution because ${\sum }_j{\alpha }_j=0$ and ${\sum }_{j}j{\alpha }_j={\rho }^{\mathrm{\prime }}(1)$. The difference $e_n-p_n$ is homogeneous and has starting size $O(h)$, so [F3] bounds it by $O(h)$ uniformly. At $n=N$, $p_N=-\delta T\mathrm{/}{\rho }^{\mathrm{\prime }}(1)\mathrm{\ne }0$, contradicting convergence. Hence ${\rho }^{\mathrm{\prime }}(1)=\sigma (1)$ is necessary.

The state vector after $n$ homogeneous steps is $C^n$ times the initial state. A uniform bound on these powers therefore gives the zero-stability estimate for every $T$ and $N$. Conversely, if the powers are unbounded, choose an exponent $n$ and unit initial state with arbitrarily large amplification, then take $T=1$ and $N\ge n+k-1$ so the entire amplified history lies on the grid; this contradicts the zero-stability bound. By [F3] and the companion construction in step 1.1, the powers are uniformly bounded exactly under the root condition.

Let $E_n=\mathrm{\parallel }{U}_n-U_n^{\ast }{\mathrm{\parallel }}_{\ast }$, where $U_n^{\ast }$ is the exact history state. As long as the preceding history remains within the radius ${\eta }_0$ of step 2.1, steps 2.1 and 1.2 give $E_{n+1}\le (1+Ch)E_n+C{h}^{p+1}$. Iterating through any first possible exit from that neighborhood yields $E_n\le e^{CT}(E_0+CT{h}^p)$, since $(1+Ch{)}^N\le e^{CT}$ and $Nh=T$. Starting errors tending to zero make the right side tend to zero; starters of order $p$ make it $O(h^p)$. Since $\mathrm{\parallel }{U}_n-U_n^{\ast }{\mathrm{\parallel }}_{\mathrm{\infty }}\le E_n$, for sufficiently small $h$ these bounds are strictly below ${\eta }_0$, so a first exit is impossible. The updates therefore exist through $t_N$, and the same estimate gives the asserted convergence and order.

Steps 1.2, 2.1, and 3.1 prove sufficiency and its order-$p$ estimate; steps 1.3, 1.4, and 2.2 prove necessity of the root condition and both consistency identities. Together with [F1] and step 2.3, these are exactly consistency and zero-stability, which proves the Dahlquist equivalence.

With $f=0$, the recurrence is $(1+2r_n)Y_{n+1}-(1+r_n{)}^2{Y}_n+r_n^2{Y}_{n-1}=0$. Substitute $Y_{n+1}=Y_n+{\mathrm{\Delta }}_{n+1}$ and $Y_{n-1}=Y_n-{\mathrm{\Delta }}_n$; the coefficient of $Y_n$ cancels, leaving $(1+2r_n){\mathrm{\Delta }}_{n+1}=r_n^2{\mathrm{\Delta }}_n$.

Since $1+2r_n>0$, division gives ${\mathrm{\Delta }}_{n+1}=q(r_n){\mathrm{\Delta }}_n$. The derivative $q^{\mathrm{\prime }}(r)=2r(1+r)\mathrm{/}(1+2r{)}^2$ is positive for $r>0$, while $q(r)\le 1$ is equivalent to $r^2-2r-1\le 0$, or $r\le 1+\sqrt{2}$.

The function $q$ is continuous and increasing, so for $r_n\le r_{\ast }<1+\sqrt{2}$ its maximum is $q_{\ast }<1$. This proves the uniform contraction claim.