A-Stability and the Dahlquist Barrier

The step of an integrator is bounded twice: once by accuracy, once by stability. For non-stiff problems the accuracy bound is usually tighter, so a stability analysis looks academic. For stiff problems the reverse holds, and the shape of the stability region becomes the property that decides whether a method is usable at all. The strongest requirement one can impose is that a decaying mode be damped for every step size, and the theory that follows from that requirement contains a hard limit on the order of multistep methods.

The stability region as a set in the complex plane

Apply a one-step method to the scalar test equation $y’ = \lambda y$. Any such method produces $y_{n+1} = R(h\lambda) y_n$ with $R$ determined by its coefficients. The region of absolute stability is the set of $z = h\lambda \in \mathbb{C}$ for which $|R(z)| \le 1$. A method is A-stable when this set contains the entire open left half plane, ${\operatorname{Re} z < 0}$. For a linear multistep method the recurrence is not a single rational function but a polynomial relation between successive values, and the stability region is defined instead by the requirement that the roots of the characteristic polynomial $\sum_{j=0}^{k} (\alpha_j - z\beta_j)\zeta^{,j} = 0$ lie in the closed unit disc, with $\alpha_j, \beta_j$ the coefficients of the method.

For explicit Euler, $R(z) = 1 + z$, and the region is the disc $|1 + z| \le 1$, of radius 1 centered at $-1$. Implicit Euler, $R(z) = 1/(1-z)$, has the region $|1 - z| \ge 1$: everything outside a disc of radius 1 centered at $+1$, which includes the whole left half plane. The trapezoid rule, $R(z) = (1 + z/2)/(1 - z/2)$, has $|R| \le 1$ precisely in the left half plane and nowhere else in the interior of the right. Each of the three has a different fate for the mode at $z \to -\infty$: explicit Euler diverges, implicit Euler sends it to zero, the trapezoid rule sends it to $-1$, which is bounded but not damped.

What A-stability does not cover: the mode at infinity

The behavior at $z \to -\infty$ matters because a discretized parabolic operator produces eigenvalues with magnitude scaling like $1/\Delta x^2$, and the stiffest of them should be annihilated, not merely kept bounded. A method is L-stable when it is A-stable and $R(z) \to 0$ as $z \to -\infty$. Backward Euler is L-stable; the trapezoid rule is not. The practical consequence appears in the first steps after an abrupt change in initial condition: an L-stable method kills the incompatible fast transient in one step, while the trapezoid rule rings, oscillating with amplitude that decays slowly because no eigenvalue of the discretization damps the mode it maps near $-1$.

A refinement, stiff decay, quantifies the same property by requiring that the amplification of the stiffest modes be small rather than merely bounded. For a linear multistep method the condition is stated in terms of the root with the largest magnitude in the limit $z \to -\infty$; for Runge-Kutta methods it is stated as a bound on $|R(z)|$ for large negative real $z$.

The point of the distinction is that A-stability, taken literally over the whole left half plane, constrains the method far more than stiff problems strictly require, which is why methods that are A-stable but not L-stable are useful for some problems and unsuitable for others.

The second Dahlquist barrier

The central result for multistep methods is a pair of theorems. The first barrier states that an $A$-stable linear multistep method has order at most two, if it also satisfies the requirement that the error constant be small; the precise statement is that an A-stable method of order greater than two has an error constant exceeding a negative bound, so it is “A-stable but of no practical use” for accuracy. The second and the decisive one states without qualification: a linear multistep method that is A-stable has order at most two.

The second-order methods that achieve it are the trapezoid rule and the two-step BDF, and both come with caveats. The trapezoid rule is not L-stable. The two-step BDF has a stability region that includes the left half plane but not in an entirely satisfactory way. Its behavior near the imaginary axis is poor, which matters for problems with oscillatory eigenvalues, and its stability is limited in a way that makes it unattractive relative to backward differentiation formulas of higher step count.

The step beyond order two requires giving something up. The backward differentiation formulas form the standard family: apply the multistep rule with $\beta_j = 0$ for $j < k$ so that the right-hand side is evaluated only at the newest point, making the method implicit in one unknown per step. BDF of order $k$ has order $k$, and its stability region covers the left half plane for $k = 1, \ldots, 6$ except for angles near the imaginary axis, which grow with $k$. Order 6 is the last stable one; BDF7 has no stability region containing the negative real axis at all, and it is useless for stiff problems. This is the price the barrier forces: orders above two are available for multistep methods only by retreating from full A-stability, and the retreat is quantified as an angular loss near the imaginary axis.

The angular loss is not a small correction. The included angle around the negative real axis shrinks as the order rises, which is the same statement as saying that the family of stiff problems for which the method remains stable narrows. A problem with eigenvalues near the imaginary axis needs the angle to cover them, and past the limit the method becomes unstable rather than inaccurate. That is the practical form of the barrier: not a warning about order, but a geometric condition that can be checked against a scatter plot of the problem’s eigenvalues before any integration runs.

The check requires knowing the engine that produces those eigenvalues. For a scheme derived from a spatial discretization, the Jacobian inherits the stencil’s spectrum, and for a problem with both diffusion and advection the spectrum is a mixture of a real negative part from the diffusion and an imaginary part from the advection. Adding diffusion pushes eigenvalues leftward, which widens the angle available for a given BDF order; removing it leaves the spectrum on the imaginary axis, where the angular loss bites immediately. The parabolic-only case in explicit diffusion and the parabolic step limit is the one where the spectrum is real and the BDF family is at its most favorable.

Why Runge-Kutta methods escape the barrier

The barrier applies to linear multistep methods, not to one-step methods with more stages. An implicit Runge-Kutta method’s stability function $R(z)$ is a rational approximation to $e^z$ whose numerator and denominator degrees are bounded by the stage count, and by choosing coefficients it is possible to make a rational function of, say, degree 2 over degree 2 that is A-stable and of order 4, exhibited by the two-stage Gauss-Legendre collocation method, which is A-stable, symplectic, and of order $2s$ for $s$ stages applied to non-stiff problems. Its stability function is the diagonal Padé approximant to $e^z$, and those approximants are A-stable for all degrees.

Higher-stage Radau methods are L-stable at odd stage counts, which makes them standard choices for stiff problems; the price is that all stages couple and the implicit system is full rather than triangular, so each step requires a Newton solve on a system $s$ times larger than the state. The Newton structure and its convergence conditions are treated in Newton iterations for the implicit spatial operator. The two-stage Gauss method has $R(z)$ of modulus exactly 1 on the imaginary axis, making it nondissipative for oscillatory problems but not damping, a distinction that matters when the problem has modes that should decay.

What the barrier forbids for one class it permits for another because of a structural difference in the stability analysis. A multistep method’s characteristic polynomial is constrained by the order conditions and by the requirement that all $k$ roots lie inside the unit disc for every $z$ in the left half plane, and the root locus available under those constraints is too poor. A one-step method with stages has the freedom of a rational function with independent numerator and denominator, and the Padé family supplies all the necessary rational approximations.

Reading this as a design constraint

The practical translation is a decision table. If the problem is non-stiff, order and accuracy matter first and A-stability is irrelevant. If the problem is stiff and accuracy demands more than second order in a multistep scheme, BDF orders up to six are the normal choice, and the angular loss near the imaginary axis should be checked against the eigenvalues of the problem’s Jacobian. If the problem is stiff with significant oscillatory modes, the BDF angular loss is a real limitation, and implicit Runge-Kutta with L-stability, or a specialized scheme, becomes appropriate. If the problem’s stiff eigenvalues are all real and negative, the BDF angular loss is harmless.

The barrier is a statement about the worst case, and its practical weight depends on how much of the left half plane the application actually occupies. The standard way to see where a given discretization sits is the eigenvalue plot: eigenvalues of the Jacobian in a plane, with the method’s stability region superimposed. The analysis of the eigenvalues that justify calling a problem stiff in the first place is in stiffness, and the step-size consequences of staying inside the region are in step-size control.