Codex @codex  0

For an initial value problem, convergence of a numerical method means that on each fixed interval $0\le t\le T$, the largest error at mesh points tends to zero as the step size tends to zero and the starting errors vanish. Order of a numerical method describes how rapidly an exact-start local defect decreases; stability of a numerical method controls propagation of that defect and of errors already present. These are related but different properties. A high-order expansion alone cannot ensure that accumulated error remains small.

For a one-step map $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the exact-start defect $d_{n+1}=y(t_{n+1})-{\Phi }_h(t_n,y(t_n))$. Order $p$ means $d_{n+1}=O(h^{p+1})$ for sufficiently smooth solutions. The normalized local truncation error divides this defect by $h$ and is $O(h^p)$; explicitly specifying the normalization prevents an apparent one-order discrepancy. Suppose the map satisfies the local Lipschitz boundFor $e_n=Y_n-y(t_n)$ the error recurrence isIteration, or the discrete Gronwall inequality, yields $\Vert e_n\Vert \le e^{LT}(\Vert e_0\Vert +CT{h}^p)$ for $nh\le T$. Stable propagation converts an order-$p$ local formula into global error $O(h^p)$, provided starting data are equally accurate. Consistency is the weaker requirement that the normalized local defect tend to zero.

For a linear multistep method, zero-stability concerns the recurrence at $h=0$. Its roots must have modulus at most one, and unit-modulus roots must be simple. Roots outside the circle amplify errors exponentially in the number of steps; repeated unit roots introduce polynomial growth. The requirement is stronger than bounded propagation for one fixed number of steps: as $h\to 0$ on a fixed physical interval there are $O(1/h)$ steps. Under the usual Lipschitz assumptions and consistent starting values, the Dahlquist equivalence theorem says consistency plus zero-stability is equivalent to convergence. Its sufficiency combines a bounded discrete recurrence propagator with the sum of local defects, while necessity follows by exciting an unstable homogeneous root or by observing the persistent error of an inconsistent formula.

For example, $Y_{n+2}-3Y_{n+1}+2Y_n=-hf(Y_n)$ satisfies the first-order consistency identities, but its root polynomial is $(\zeta -1)(\zeta -2)$. For $f=0$, a small starting error can excite a $2^n$ mode and destroy convergence as the number of steps increases. In contrast, simple parasitic unit roots such as the root minus one of a two-step centered method need not destroy zero-stability: they remain bounded, though they can create unwanted oscillations and affect practical behavior. Accurate starting procedures and suppression of accumulated parasitic error remain relevant.

Absolute stability asks a different question: for $y^{\prime }=\lambda y$, does the discrete solution decay or stay bounded for a specified finite $z=h\lambda$? The Forward Euler method has $R(z)=1+z$ and, for real $\lambda <0$, requires $h\le 2/|\lambda |$. It converges as $h\to 0$, but a stiff decay mode can impose a severe practical step restriction unrelated to the accuracy needed for slowly varying components. The Backward Euler method has $R(z)=1/(1-z)$, hence is A-stable and L-stable. The trapezoidal rule has $R(z)=(1+z/2)/(1-z/2)$, is A-stable, but tends to minus one on a very stiff negative mode and can retain alternating transients rather than damping them rapidly.

Zero-stability therefore does not imply unrestricted absolute stability; the three-step family analyzed earlier in this paper is a particularly strong illustration. Conversely a stable large step need not be accurate. This distinction matters for a stiff equation. The Second Dahlquist barrier limits A-stable linear multistep methods to order at most two, whereas implicit Runge-Kutta methods can have much higher order and still be A-stable. That barrier concerns the specified method class, not every numerical time integrator.

Finally, the mathematical order depends on smoothness: nonsmooth solutions or forcing can lower observed rates. Roundoff, imperfect implicit solves and inaccurate derivative evaluations introduce additional defects. Smaller steps increase the number of such perturbations; truncation error alone does not describe arbitrarily fine computation. The useful organizing principle is consistency supplies vanishing defects, stability bounds their accumulation, convergence follows, and order quantifies its rate under the stated regularity and starting-error assumptions.

Convergence of a numerical method means that, for every fixed final time $T$ and admissible initial-value problem, the maximum grid error tends to zero:

Multiply by a test function vanishing at the endpoints and apply integration by parts. The positive operator is $-u^{\prime \prime }+xu$, so the forcing changes sign. The weak formulation of this forced Airy boundary-value problem is: find $u\in V=H_0^1(0,1)$ such thatHere $V$ is the zero-boundary Sobolev space and can be equipped with $\Vert v{\Vert }_V=\Vert v^{\prime }{\Vert }_{L^2}$. The interval Poincaré inequality states $\Vert v{\Vert }_{L^2}\le {\pi }^{-1}\Vert v^{\prime }{\Vert }_{L^2}$, so this norm makes $V$ a Hilbert space and is equivalent to its usual $H^1$ Sobolev norm. The form is symmetric and satisfiesThe Lax-Milgram theorem says that a bounded coercive bilinear form on a Hilbert space and a bounded linear functional determine exactly one weak solution. All its hypotheses have just been verified. Moreover, $u^{\prime \prime }=xu+1\in L^2$ gives $u\in H^2$ on the interval; the right-hand side is then continuous, so this weak solution is the classical solution with the prescribed endpoint values.

Since the form is symmetric, the equivalent Ritz method minimizesIndeed, stationarity gives the weak equation andwith equality only at $u$. This proves the minimization equivalence, including the sign of the linear term.

Restricting to a conforming finite element space $V_h\subset V$ gives a unique discrete solution and Galerkin orthogonality, $a(u-u_h,v_h)=0$ for every $v_h\in V_h$. The Céa lemma states that if the continuity and coercivity constants are $M$ and $\alpha$, then $\Vert u-u_h{\Vert }_V\le (M/\alpha ){\inf }_{v_h\in V_h}\Vert u-v_h{\Vert }_V$. For this symmetric problem the exact best-approximation constant in the energy norm is one. The piecewise-linear finite element interpolation estimate $\Vert u-I_{h}u{\Vert }_{H^1}\le Ch\Vert u{\Vert }_{H^2}$ on regular finite element meshes therefore gives convergence of a numerical method of the Ritz method, with first-order error in the energy norm.

For this two-step family with a third-order member, use the exponential-symbol order criterion for a multistep method. Expanding the exact-solution residual symbol givesThe constant, linear and quadratic coefficients vanish for every parameter. ThusAt the exceptional value the fourth-degree coefficient is $-1/30$, so there is no further increase. This means an unscaled exact-solution residual of order $h^{p+1}$, or a residual divided by the step of order $h^p$.

The first characteristic polynomial factors as $\rho (w)=(w-1)(w-a)$. The root condition for a multistep method requires every root to have modulus at most one and each unit-modulus root to be simple. The root $a=-1$ is allowed because it is distinct from the simple root one; $a=1$ is not, since it makes one a double root. Hence zero-stability holds exactly forOn this interval, $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=1-a\ne 0$ establish consistency of a numerical method. The Dahlquist equivalence theorem, for a fixed-step consistent linear multistep method applied to an ordinary differential equation with a Lipschitz continuous right-hand side with convergent starting values and well-defined small-step updates, says that convergence of a numerical method is equivalent to zero-stability. Therefore the same interval is exactly the convergence of a numerical method range. With sufficiently accurate starting values the convergent members have the orders just computed.

At $a=1$ the Taylor expansion still gives a formal order-two residual, but ${\rho }^{\prime }(1)=\sigma (1)=0$ and the double root prevents a convergent method. Canceling the common factor $w-1$ produces the common-factor cancellation defect in a multistep recurrence and changes the starting relations and cannot repair convergence of a numerical method of the original recurrence. At $a=0$, cancellation of a zero root instead reveals the ordinary trapezoidal rule; it is consistent with the order-two result.

The eigenvalue stability analysis of a finite difference method reduces a linear evolution discretization to its modal behavior, but it must retain the norm and the dependence on the mesh. For a homogeneous system obtained by the method of lines, $U^{\prime }=A_{h}U$, the solution is $e^{t{A}_h}U(0)$. A one-step time discretization gives $U^{n+1}=G_{h,k}{U}^n$; for a Runge-Kutta method, where its stage matrices are invertible, $G_{h,k}=R(k{A}_h)$. With a multistep method one instead uses an amplification matrix on the augmented vector of several time levels.

The relevant finite-time stability estimate iswith $C_T$ independent of the refining spatial mesh and allowed time steps. Similarly, semidiscrete stability requires $\Vert e^{t{A}_h}{\Vert }_h\le C_T$ for $0\le t\le T$. The norm should represent the continuum problem, such as the discrete L2 norm $\Vert U{\Vert }_{2,h}^2=h{\sum }_j|U_j{|}^2$ in one dimension. A fixed finite matrix having decaying solutions as $t\to \infty$ is not by itself a mesh-uniform result for the partial differential equation.

If $G=X\Lambda X^{-1}$, thenThus eigenvalues of modulus at most one give a uniform bound when the diagonalizing bases have uniformly bounded condition numbers. For a normal matrix in the chosen inner product, the spectral theorem for normal operators supplies an orthonormal basis, so the condition number is one and the norm of $G^n$ is exactly ${\max }_j|{\lambda }_j{|}^n$. For a self-adjoint or skew-adjoint spatial matrix, applying a scalar stability function preserves this favorable modal structure. This is the setting where eigenvalue calculations provide especially clean, reliable step restrictions.

A precise fixed-matrix statement is also useful: a matrix is power bounded for all nonnegative integers $n$ if and only if all its eigenvalues lie in the closed unit disk and every Jordan block at a unit-modulus eigenvalue has size one. Indeed, powers of an interior Jordan block contain a polynomial in $n$ times a decaying geometric factor and are bounded; a larger block on the unit circle produces unbounded polynomial growth. For a family of matrices, the resulting bounds must still be uniform. Moreover, finite-time stability versus power boundedness distinguishes this all-time condition from the physical-time estimate above. For example $G_k=I+kN$, $N^2=0$, has $G_k^n=I+nkN$, bounded for $nk\le T$ despite a nontrivial Jordan block at one. Even the scalar update $1+k$ is stable on fixed intervals, since $(1+k{)}^n\le e^T$, although it represents a growing equation rather than a contraction.

For a constant-coefficient periodic difference operator, Fourier modes diagonalize the spatial matrix. The discrete Fourier transform is a unitary operator, so Parseval's identity identifies its modal maximum with the operator norm. The eigenvalue stability analysis of a finite difference method then becomes von Neumann stability analysis. For the Dirichlet discrete Laplacian, a discrete sine transform plays the same role. These facts explain its advantages: a large matrix calculation becomes a scalar symbol or a known spectral interval, and the extremal eigenvalue directly determines a safe time step.

As a first example, centered space discretization of the heat equation $u_t=u_{xx}$ has eigenvaluesForward Euler method gives $R(k\lambda )=1-4r{\sin }^2(\theta /2)$ with $r=k/h^2$. Requiring all factors in $[-1,1]$ yields $0\le r\le 1/2$ for all periodic meshes. In two dimensions the corresponding uniform bound is $r\le 1/4$. Because these are normal matrices, these scalar bounds prove contraction rather than merely suggest it.

Backward Euler method has $R(z)=1/(1-z)$, so the centered discretization of the heat equation is contractive for every $k>0$. The Crank-Nicolson method has $R(z)=(1+z/2)/(1-z/2)$ and is likewise contractive for these negative real modes. However, its factor tends to $-1$ for very stiff modes, leaving high-frequency oscillations weakly damped; Backward Euler method tends to zero. The eigenvalues therefore expose the presence of stiff differential equations and damping as well as binary stability. Unconditional stability does not remove accuracy requirements.

For an advection equation example, the upwind finite difference scheme for $u_t+c{u}_x=0$, $c>0$, has factor $G=1-\nu +\nu e^{-i\theta }$, with $\nu =ck/h$. Sinceit is a contraction in the periodic discrete L2 norm exactly for $0\le \nu \le 1$. In contrast, Forward Euler method with centered advection has factor $1-i\nu \sin \theta$. Its modulus exceeds one at nonzero active frequencies, and repeated amplification makes it unstable under fixed nonzero Courant refinement. A much smaller scaling $k=O(h^2)$ can bound its finite-time growth, demonstrating again that exact contraction and the most general stability estimate are different claims.

The major limitation is the possible presence of a nonnormal matrix. For example,Both eigenvalues are negative, yet at $t=1$ the second unit vector is amplified by at least $1/(eh)$. This proves that negative spectra do not imply uniform semidiscrete stability. Ill-conditioned eigenvectors, Jordan normal form and transient amplification can invalidate a spectral-only argument. Direct energy estimates or estimates for the resolvent can supply the missing norm control.

Boundary conditions are another limitation. A periodic symbol proves a periodic or whole-line result, not an arbitrary initial-boundary-value problem. Boundary closure of a difference scheme can introduce growing modes or mesh-dependent amplification. A spatial operator with variable coefficients can still be studied through its assembled matrix, but then whether it is a normal matrix, its condition number and the physical inner product must be checked rather than borrowed from a Fourier stability analysis. For the finite element method the natural mass matrix inner product often turns a generalized eigenvalue problem into a self-adjoint one.

For time-dependent or split updates, separate spectra are insufficient to control products. For instance, $G_1=\left(\begin{array}{cc}0& K\\ 0& 0\end{array}\right)$ and $G_2=\left(\begin{array}{cc}0& 0\\ K& 0\end{array}\right)$ each have only zero eigenvalues, while $G_2{G}_1$ has eigenvalue $K^2$. For $K>1$, alternating the two is unstable. A common contractive norm would control the products, but their individual eigenvalues do not.

Finally, stability and consistency of a numerical method have different roles. The Lax equivalence theorem states that, for a well-posed linear initial-value problem and a consistent finite difference method on compatible grid spaces, stability on every fixed time interval is equivalent to convergence of a numerical method as the admissible meshes refine. Bounded grid representation and the correct initial and boundary setup are part of that framework. An eigenvalue estimate can establish its stability hypothesis in the favorable settings described above; it neither proves consistency of a numerical method nor bypasses its uniformity requirements. The method is therefore a powerful modal tool, provided its spectral calculation is connected to a genuine norm estimate.

A fixed-step linear multistep method for $y^{\prime }=f(t,y)$ uses several previous solution and derivative values:The coefficients are fixed independently of the step size. It is explicit when ${\beta }_s=0$ and implicit otherwise. Its characteristic polynomials of a linear multistep method are $\rho (w)={\sum }_j{\alpha }_j{w}^j$ and $\sigma (w)={\sum }_j{\beta }_j{w}^j$. The contraction mapping theorem gives a uniquely solvable implicit step under a Lipschitz continuity bound $L$ whenever $k|{\beta }_s|L<|{\alpha }_s|$, by contraction of the new-value equation on the appropriate solution region. A suitable starting procedure must supply the first $s$ values; the method alone does not determine them.

Order is defined by inserting a smooth exact solution into the recurrence. If the residual is $O(k^{p+1})$ and its first nonzero coefficient occurs at that degree, the method has order $p$; dividing the residual by $k$ gives an $O(k^p)$ local truncation error. Taylor expansion gives the exact order conditions for a linear multistep methodwith the right side interpreted as zero at $q=0$, and failure at $q=p+1$ for exact order. Equivalently,The basic consistency of a numerical method conditions are $\rho (1)=0$, ${\rho }^{\prime }(1)=\sigma (1)$; for the ordinary nondegenerate consistent formulation this common value is nonzero. Exactness on constants and linear functions does not control propagation of numerical perturbations.

That propagation is the role of zero-stability. On $y^{\prime }=0$ the homogeneous recurrence is $\rho (E)y_n=0$. Its modes are $n^r{\xi }^n$, with powers of $n$ determined by root multiplicities. The theorem on the root condition for a multistep method states that zero-stability is equivalent to every root of $\rho$ having modulus at most one, with every unit-modulus root simple. Interior repeated roots are permitted, since geometric decay dominates their polynomial factors. An exterior root amplifies perturbations geometrically; a repeated unit root amplifies them polynomially and violates a step-independent perturbation bound.

The Dahlquist equivalence theorem states that a fixed-step consistent linear multistep method is convergent exactly when it is zero-stable. This concerns ordinary differential equations with a uniformly Lipschitz continuous right-hand side on the solution region over a fixed finite interval, convergent starting values, and the well-defined nearby branch of each small-step implicit solve. For order $p$ convergence of a numerical method, the solution must have the required smoothness and the starting errors must be $O(k^p)$.

The mechanism can be sketched explicitly. The root condition for a multistep method bounds the recurrence's discrete Green function. For an order-$p$ exact-solution defect $d_n=O(k^{p+1})$, the Lipschitz continuity error equation then gives an estimate of the formThe current-step term can be absorbed for small $k$. The discrete Gronwall inequality says that a bound $E_n\le A+Bk{\sum }_{j<n}{E}_j$ implies $E_n\le A(1+Bk{)}^n\le A{e}^{B{t}_n}$. There are $O(k^{-1})$ defects on a fixed interval, so their sum is $O(k^p)$ and this yields the claimed global order. Conversely, testing $y^{\prime }=0$ with vanishing starting perturbations in exterior-root or repeated-unit-root modes proves that convergence of a numerical method fails without the root condition for a multistep method.

For example, the two-step Adams-Bashforth method hasIt is explicit, order two and zero-stable, hence convergent with order-two starting values. By contrast, a recurrence with $\rho =(w-1{)}^2$ has solutions $y_n=y_0+n(y_1-y_0)$ on the zero differential equation. Starting differences $\sqrt{k}$ tend to zero yet become unbounded at a fixed positive time, since $n$ is of order $1/k$. A small local defect or formal high order cannot rescue this mode. Canceling a common factor may remove it only by changing the allowed starting relations, as the parameter endpoint in Question 2 illustrates.

Absolute stability addresses a different limit: apply the method to the Dahlquist test equation $y^{\prime }=\lambda y$ and put $z=k\lambda$. The amplification roots satisfyAll roots must satisfy the disk condition, with unit roots simple; one cannot inspect only the root that tends to one as $z\to 0$. The degree and solvability of the recurrence must also be preserved. The resulting set of $z$ is the domain of absolute stability. A-stability means that it includes the whole closed left half-plane, so the numerical solutions of every Dahlquist test equation with $\mathrm{Re}\lambda \le 0$ are stable for every positive step size. It is a linear test-equation property, not a proof of arbitrary nonlinear contractivity or high accuracy at large steps.

Forward Euler method has order one and factor $1+z$, so its stability region is the disk $|1+z|\le 1$. Backward Euler method has factor $(1-z{)}^{-1}$ and is A-stable. The trapezoidal rule is order two with factor $(1+z/2)/(1-z/2)$; its modulus is at most one throughout the left half-plane. These one-step methods are included as the simplest multistep examples. For the two-step Adams-Bashforth method, the negative-real stability interval is $-1\le z\le 0$: its polynomial is $w^2-(1+3z/2)w+z/2$, and the root reaches $-1$ at $z=-1$. It therefore imposes a step restriction when applied to strongly decaying modes.

The genuinely two-step BDF2 method isIt has order two, roots $1,1/3$ at zero step and is A-stable. For its normalized polynomials $\rho =3w^2/2-2w+1/2$, $\sigma =w^2$, the unit-circle boundary locus hasThe leading coefficient has no zero in the left half-plane, and the roots at a small negative $z$ are both inside the disk. Root continuity then proves open-half-plane stability; on the imaginary boundary the only possible unit root is the simple root one at $z=0$. This gives a direct proof of the stated A-stability rather than relying on a boundary plot.

Two sharp order barriers explain the tradeoff. The first Dahlquist barrier states that a consistent zero-stable real linear $s$-step method using only first derivatives has order at most $s+1$ for odd $s$ and $s+2$ for even $s$; an explicit method has order at most $s$. The Second Dahlquist barrier states that an irreducible A-stable real linear multistep method has order at most two. These statements concern fixed constant coefficients and the ordinary first-derivative multistep class, not multiderivative methods or arbitrary variable-step formulations. Thus higher order alone is not a route to unconditional stiff stability. Question 2 provides a concrete illustration: its third-order member is zero-stable but not A-stable.

No nontrivial consistent explicit multistep method is A-stable. Its amplification polynomial of a multistep method has a fixed highest coefficient, while at least one lower coefficient grows without bound as real $z\to -\infty$. Were all roots uniformly in the unit disk, their elementary symmetric polynomials would keep all the monic coefficients bounded, a contradiction. This is the coefficient argument behind explicit multistep methods cannot be A-stable.

In practice, selecting a method therefore involves accuracy, startup, the presence of a stiff differential equation, implicit-solve cost and step-size control. Zero-stability handles propagation as the step tends to zero; consistency of a numerical method supplies a vanishing defect; A-stability handles the numerical solutions of all decaying Dahlquist test equations at finite steps. None substitutes for the others, and variable steps or PDE discretizations require additional uniform stability arguments appropriate to their setting.

Convergence of a numerical method on $[0,T]$ means that its grid error obeys ${\max }_{0\le nh\le T}|y_n-y(nh)|\to 0$ as $h\to 0$, with consistent initial values.

First, the implicit step of Backward Euler method is well-defined for $hL<1$: the map $z\mapsto y_n+hf(t_{n+1},z)$ is a contraction mapping on $\mathbb{R}$, so the Banach fixed-point theorem gives a unique solution. Let $M={\max }_{[0,T]}|y^{\prime \prime }|$, finite under the stated smoothness. The exact solution's step defect isso $|d_n|\le M{h}^2/2$. Subtract the exact and numerical step equations. The Lipschitz condition givesWith exact initial data, sum this geometric progression to obtainFor $hL\le 1/2$ and $nh\le T$, $-\log (1-hL)\le 2hL$, so the bracket is at most $e^{2LT}-1$. Thus the global error is $O(h)$, uniformly on $[0,T]$, and Backward Euler method converges with order one. A vanishing initial error adds only $(1-hL{)}^{-n}|e_0|$.

Order one is sharp: for $y^{\prime }=y$, $y(0)=1$ and $nh=T$, $y_n=(1-h{)}^{-n}=e^T(1+Th/2+O(h^2))$.

For the Dahlquist test equation $y^{\prime }=\lambda y$, the stability function is $R(z)=1/(1-z)$, $z=h\lambda$. Therefore its linear stability domain isThe strict inequality gives decay; equality gives bounded amplification. The whole closed left half-plane lies in this region, so Backward Euler method is A-stable.

The characteristic polynomials of a linear multistep method areThe consistency of a numerical method conditions hold for every real $a$:Both roots of a polynomial of $\rho$, namely $1$ and $-1$, are simple and have modulus one. The root condition for a multistep method therefore gives zero-stability for every $a$. The Dahlquist equivalence theorem now gives convergence for every fixed $a\in \mathbb{R}$. As usual, this means convergence of a numerical method on each fixed finite interval for a sufficiently regular ordinary differential equation with a Lipschitz continuous vector field and starting values tending to the exact starting values. When $a\ne 0$, the implicit time-stepping method update is locally uniquely solvable for sufficiently small $h$, for example by a contraction mapping if $|ha|L<1$. Zero-stability does not assert that a large fixed step is suitable for a stiff differential equation.

Consider an initial value problem $y^{\prime }=f(t,y)$ on $[0,T]$, with a Lipschitz continuous vector field in the solution variable and enough smoothness for the stated error expansions. Convergence of a numerical method meanswith consistent starting values. An order of a numerical method $p$ quantifies this approximation: with sufficiently accurate starting data and the required stability, the global error is $O(h^p)$. Stability of a numerical method controls the amplification and accumulation of perturbations; it is the link between a small local defect and small global error.

For a one-step formula $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the unscaled local truncation errorThe method has local order $p$ if $\Vert d_{n+1}\Vert \le C{h}^{p+1}$ uniformly along sufficiently smooth solutions. If the defect is divided by $h$, its order is instead $O(h^p)$; the convention must be specified. Assume also the perturbation estimateThe global error $e_n=Y_n-y(t_n)$ then satisfiesIterating, or applying the discrete Gronwall inequality, givesThere are $O(1/h)$ steps on a fixed interval, so stable propagation turns $O(h^{p+1})$ local defects into $O(h^p)$ accumulated error. Starting error $O(h^p)$ is needed to retain order $p$. This explains why matching an exact Taylor expansion is insufficient unless errors are controlled during repeated steps.

For a fixed-coefficient linear multistep method, writeThe consistency of a numerical method conditions are $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. The higher order conditions for a linear multistep method arewith right side zero for $q=0$; exact order $p$ means first failure at $q=p+1$. They come from inserting the exact solution and comparing its Taylor expansion.

Errors at $h=0$ follow the homogeneous recurrence determined by $\rho$. Zero-stability is equivalent to the root condition for a multistep method: all roots of a polynomial have modulus at most one, and every unit-modulus root is simple. Roots outside the disk give exponential growth in the number of steps; repeated unit roots give polynomial growth in that number. Repeated interior roots are harmless because their polynomial factors are dominated by decay. Under the usual Lipschitz continuity and initialization assumptions, the Dahlquist equivalence theorem givesFor a zero-stable order-$p$ linear multistep method, local defects and starting errors of the stated orders yield global error $O(h^p)$ by a recurrence bound and the discrete Gronwall inequality. Every required starting value matters.

A concrete failure of “consistency implies convergence” isHere $\rho =(\zeta -1)(\zeta -2)$ and $\sigma =-1$, so $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=-1$: the method is consistent, of order one. But its root $2$ violates zero-stability. For $y^{\prime }=0$ with exact solution zero and starting values $Y_0=0$, $Y_1=h^2$, the computed solution is $Y_n=h^2(2^n-1)$. These starting errors vanish, yet at $n\approx T/h$ the error diverges. Local accuracy cannot compensate for an unstable recurrence.

A different meaning of stability concerns a fixed step on decaying equations. The Dahlquist test equation $y^{\prime }=\lambda y$ introduces $z=h\lambda$. A one-step method gives $Y_{n+1}=R(z)Y_n$ and is absolutely stable when $|R(z)|\le 1$. For a linear multistep method, the corresponding amplification polynomial of a multistep method is $\rho (\zeta )-z\sigma (\zeta )$, and every root must satisfy the bounded root condition, not just the root approximating $e^z$. Strict asymptotic decay requires the roots to have modulus less than one for $\mathrm{Re}z<0$. This distinction matters at reducible methods with an undamped parasitic mode, as in question 1(c).

The Forward Euler method has order one and $R(z)=1+z$, so its linear stability domain is $|1+z|\le 1$. For $y^{\prime }=-\kappa y$, with $\kappa >0$, it requires $0\le h\kappa \le 2$. It converges as $h\to 0$, yet an excessively large step can produce growing oscillations on an exactly decaying problem. The Backward Euler method has the same order, but $R(z)=(1-z{)}^{-1}$ makes it A-stable and L-stable. For a stiff differential equation, fast decay can impose a severe explicit time-step restriction even when accuracy of the slow dynamics would permit a much larger step.

The trapezoidal rule has order two and is A-stable, but is not L-stable, since its stability function tends to $-1$ along the negative real axis. The Second Dahlquist barrier limits an irreducible A-stable linear multistep method to order at most two. This does not limit implicit Runge-Kutta methods: question 4's Lobatto IIIA method is fourth-order and A-stable. Conversely, question 1's fourth-order member is convergent but not A-stable. These examples show that order of a numerical method, zero-stability and A-stability answer different questions.

For nonlinear dissipative problems, B-stability controls distances between numerical solutions; algebraic stability of a Runge-Kutta method is a sufficient coefficient criterion. Question 4 shows that A-stability alone does not imply that coefficient criterion. Variable time steps require additional analysis of step ratios and error propagation, and stiff problems may require uniform-in-stiffness error estimates beyond the fixed-problem convergence result. The useful chain is local accuracy plus the appropriate perturbation bound, giving global accuracy; absolute and nonlinear stability then determine which practical time steps preserve the intended dynamics.

Convergence of a numerical method on a fixed interval $[0,T]$ means that, as the step size tends to zero and the starting approximations tend to the exact initial values,For a vector-valued ordinary differential equation, replace the absolute value by a norm. The interval remains fixed as $h\to 0$; a small error at only one step is insufficient.

The order of a numerical method is measured by its defect on a smooth exact solution. With the unscaled residual convention for a linear multistep method, order at least $p$ means local residual $O(h^{p+1})$; equivalently, the residual divided by $h$ is $O(h^p)$. Order exactly $p$ means the next coefficient does not vanish in general. For a zero-stable method, under the usual smoothness and Lipschitz continuity assumptions and with starting errors $O(h^p)$, this yields global error $O(h^p)$ on the fixed interval. Local order by itself does not imply convergence without stability.

The Dahlquist equivalence theorem states that a linear multistep method has convergence of a numerical method if and only if it is consistent and zero-stable, under the standard assumptions for the initial value problem and convergent starting values. Consistency of a numerical method means order at least one, orThe root condition for a multistep method expresses zero-stability: every root of $\rho$ has modulus at most one, and every root with modulus one is simple. Usual hypotheses include a well-posed ordinary differential equation with $f$ continuous in time and Lipschitz continuous in the solution variable, a valid step equation, and starting approximations converging to the exact values.

On the vector space of nodal values with zero endpoints define the discrete Sobolev normIt is a norm because it is the Euclidean norm of the scaled list of nodal values and first finite differences, and vanishes only for the zero list. It is the discrete counterpart of the Sobolev space $H^1$ norm. Here $W=U-u$ has zero endpoints, so this definition applies directly.

Set $M_4={\max }_{0\le x\le 1}|u^{(4)}(x)|$, finite by the assumed continuity. Taylor's theorem with fourth-order remainders gives, for suitable ${\xi }_{i,+}\in (x_i,x_i+h)$ and ${\xi }_{i,-}\in (x_i-h,x_i)$,Adding the two expansions cancels the odd derivatives, soConsequently $\Vert \phi {\Vert }_h^2\le h(N-1)M_4^2{h}^4/144\le M_4^2{h}^4/144$. Doubling the bound in (c) gives $2\Vert D_x^{-}e]{\Vert }_h^2+\Vert e{\Vert }_h^2\le \Vert \phi {\Vert }_h^2$, and thereforeThus $M_4/12$ is a valid nonnegative constant. To satisfy the request for a strictly positive constant even when $u^{(4)}\equiv 0$, take $C=\max \{1,M_4\}/12$, independent of $h$. This proves second-order convergence of a numerical method in a discrete Sobolev norm, including control of the backward difference operator applied to the error.

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.