43.10.05 · numerical-analysis / numerical-odes

Stiff equations, A-stability, and the Dahlquist second barrier

shipped3 tiersLean: none

Anchor (Master): Hairer-Wanner 1996 *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems* 2e (Springer) §IV.2-IV.5 (A-stability, A($\alpha$)-stability, L-stability, the second Dahlquist barrier with its order-star proof, BDF and Radau IIA stiff integrators); Dahlquist 1963 *A special stability problem for linear multistep methods* (BIT 3) (the A-stability programme and the original second-barrier proof)

Intuition Beginner

Some problems hide two clocks running at wildly different speeds. One clock ticks fast: a piece of the answer that races to its resting value and then sits still. The other clock ticks slowly: the part of the answer you actually care about, drifting along gently for a long time. A chemical reaction with one lightning-fast intermediate and one slow product is the everyday picture. Once the fast part has settled, the true curve is smooth and lazy, and a sensible person would take big, lazy steps to follow it.

A method that is stiff-friendly does exactly that. But many simple methods cannot. Even after the fast part has died away, an explicit method like forward Euler still feels the fast clock through its stability limit, and it refuses to take a big step. It forces you to keep tiptoeing with tiny steps to follow a curve that is no longer doing anything interesting. That mismatch — a smooth answer that nonetheless demands microscopic steps — is what numerical analysts call stiffness.

Stiffness is famously hard to pin down with a clean definition. There is no single number that switches from "not stiff" to "stiff." Instead it is recognised by its symptom: the stable step is far smaller than the accurate step would need to be. The wider the gap between the fastest and slowest clocks, the worse the squeeze. A problem can be mildly stiff in the morning and badly stiff by afternoon as its clocks spread apart.

The cure is to use a method whose stability does not care how fast the fast clock ticks. These are the implicit methods — backward Euler, the trapezoidal rule, the backward-difference formulas. They solve a small equation at each step instead of just plugging in, and in exchange their stability region is enormous: it swallows the whole decaying side of the picture. With one of these you can finally take the big, lazy steps the smooth answer deserves, and the fast clock, long since settled, no longer holds you hostage.

Visual Beginner

Picture two decaying signals sharing one time axis. The first plunges almost straight down and flatlines within a blink: the fast, stiff mode. The second eases down over a long stretch: the slow mode you care about. After the blink, the true answer is just the slow curve. The question is what step size each method now allows.

The table below lines up the everyday clue (how big a step the method survives once the fast mode has settled) against the kind of method. Reading down, explicit methods stay shackled to the fast mode forever, while implicit methods are freed the moment it dies.

once the fast mode has settled	explicit method (forward Euler)	implicit method (backward Euler)
what limits the step	the dead fast mode still does	only the slow accuracy does
step size allowed	stays tiny	grows large
extra work per step	none	solve a small equation
verdict on a stiff problem	inefficient, often unusable	the right tool

Each row says the same thing in a different way: implicit methods pay a little arithmetic per step to escape the tyranny of the fast clock.

Worked example Beginner

Let us watch forward Euler and backward Euler race on a problem with two clocks. The slope rule combines a fast-decaying part and a slow-decaying part, but to keep the arithmetic clean we follow just the fast piece, whose slope is $- 1000$ times its height. The true fast piece collapses to almost nothing in a flash; after that, a good method should be free to stride.

Forward Euler multiplies the height by the factor $1 + h (- 1000) = 1 - 1000 h$ each step, and stability needs the size of that factor to be at most $1$ . We start the fast piece at height $1$ .

Step 1. Find forward Euler's largest stable step. The size of $1 - 1000 h$ is at most $1$ exactly when $1 - 1000 h \geq - 1$ , that is $1000 h \leq 2$ , so $h \leq 0.002$ . Any larger step makes the computed fast piece grow and explode, even though the true fast piece is shrinking.

Step 2. Try a lazy step, $h = 0.1$ , the size of step the smooth slow part would happily accept. Forward Euler's factor is $1 - 1000 (0.1) = 1 - 100 = - 99$ , size $99$ . Starting at $1$ , the heights run $1, - 99, 9801, \dots$ , a catastrophic blow-up. Forward Euler simply cannot take the lazy step.

Step 3. Now backward Euler. It multiplies the height by $1/ (1 - h (- 1000)) = 1/ (1 + 1000 h)$ each step. At the lazy step $h = 0.1$ this factor is $1/ (1 + 100) = 1/101$ , size well under $1$ . Starting at $1$ , the heights run $1, 1/101, 1/10 1^{2}, \dots$ , fading fast and safely. Backward Euler takes the lazy step with no trouble.

Step 4. Read the gap. Forward Euler is pinned to $h \leq 0.002$ by the fast clock; backward Euler strides at $h = 0.1$ , fifty times larger, limited only by how well it tracks the slow part. That fifty-fold ratio is the stiffness penalty forward Euler pays and backward Euler escapes.

What this tells us: stiffness is not about accuracy. Forward Euler would be plenty accurate at $h = 0.1$ on the smooth slow part — it just is not allowed to take that step, because its small stability region still feels the long-dead fast mode. Backward Euler's huge stability region removes that veto, which is the whole reason implicit methods rule the stiff world.

Check your understanding Beginner

Exercise (easy, multiple choice).

Backward Euler is run on a very stiff decaying mode (slope a huge negative multiple of the height) with a large step. Its amplification factor is $1/ (1 + big)$ . What happens to the computed values?

A. They grow without bound, because the step is large. B. They fade quickly toward zero, because the factor's size is tiny. C. They stay exactly constant. D. Nothing can be said without the exact step size.

Hint

For a hugely negative slope and any positive step, the denominator $1 + big$ is large, so the factor is small.

Answer

B. They fade quickly toward zero. The factor $1/ (1 + big)$ has size far below one, so each step shrinks the value and the stiff mode is annihilated. Feedback-correct: backward Euler is stable for every stiff decaying mode at any step, which is why it survives where forward Euler explodes. Feedback-wrong: a large step does not destabilise backward Euler; its stability region covers the whole decaying side of the plane.

Formal definition Intermediate+

Throughout, the stability function $R (z)$ of a one-step method, the stability polynomial $π (ζ; z) = ρ (ζ) - z σ (ζ)$ of a linear multistep method, the rescaled variable $z = hλ$ , and the region of absolute stability $S = {z \in C : ∣ R (z) ∣ \leq 1}$ are taken from 43.10.04; the first and second characteristic polynomials $ρ, σ$ and the consistency conditions are taken from 43.10.02.

Working description (stiffness). A linear (or linearised) initial-value problem $y^{'} = A y + g (t)$ is stiff on an interval if the eigenvalues $λ_{k}$ of $A$ have $Re λ_{k} < 0$ with a large stiffness ratio $max_{k} ∣ Re λ_{k} ∣/ min_{k} ∣ Re λ_{k} ∣ ≫ 1$ , so that a numerical method whose absolute-stability region is bounded must keep $h ∣ λ_{m a x} ∣$ inside that bounded region — forcing a step far smaller than the accuracy of the smooth solution requires. No exact threshold separates stiff from non-stiff; stiffness is the regime in which absolute stability, not local accuracy, dictates the step. This is the operational characterisation used by ^{[LeVeque §8.1]} and ^{[Süli-Mayers §12.11]}, who both decline a sharp definition.

Definition (A-stability). A numerical method is A-stable if its region of absolute stability contains the entire closed left half-plane, $$ {, z \in \mathbb{C} : \mathrm{Re},z \le 0 ,} \subseteq \mathcal{S}. $$ Equivalently, the method applied to the test equation $y^{'} = λ y$ with $Re λ \leq 0$ produces a non-growing sequence for every step $h > 0$ . A-stability removes the stability cap on $h$ for any decaying or oscillatory mode, which is exactly the property a stiff integrator needs.

Definition (A( $α$ )-stability). For $0 < α \leq π /2$ , a method is A( $α$ )-stable if $S$ contains the infinite sector $$ {, z \in \mathbb{C} : |\arg(-z)| \le \alpha ,} = {, z : z = -re^{i\phi},\ r \ge 0,\ |\phi| \le \alpha ,}, $$ the wedge of half-angle $α$ about the negative real axis. A-stability is the limiting case $α = π /2$ . The backward-differentiation formulas BDF3 through BDF6 are A( $α$ )-stable with $α$ decreasing as the order rises.

Definition (L-stability). A one-step method is L-stable if it is A-stable and, in addition, its stability function satisfies $$ R(\infty) = \lim_{z \to \infty} R(z) = 0. $$ L-stability strengthens A-stability by demanding that the stiffest modes ( $∣ z ∣ \to \infty$ with $Re z \leq 0$ ) be annihilated in a single step rather than merely kept bounded. Backward Euler, with $R (z) = 1/ (1 - z)$ and $R (\infty) = 0$ , is L-stable; the trapezoidal rule, with $R (z) = (1 + z /2) / (1 - z /2)$ and $R (\infty) = - 1$ , is A-stable but not L-stable.

The predicates A-stable, A( $α$ )-stable, L-stable, the stiffness ratio, and the limit value $R (\infty)$ are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

"Stiffness is just a large Lipschitz constant." A large $∣ λ ∣$ alone does not make a problem stiff: a single fast mode integrated from a generic initial condition is genuinely fast and should be resolved with small steps. Stiffness needs the separation of scales — a fast mode that has already decayed while a slow mode still demands integration — so the stable step is wasted on a transient that is numerically dead.
"A-stable implies the stiffest modes are killed." A-stability only keeps modes from growing; it permits $∣ R (\infty) ∣$ up to one. The trapezoidal rule has $R (\infty) = - 1$ , so a very stiff transient is damped only marginally and flips sign each step, producing the ringing that L-stability (here failing) would prevent.
"Higher order is always better, so use a high-order A-stable multistep method." No A-stable linear multistep method has order above two — the second Dahlquist barrier below. Pushing a multistep method past order two forces its region to leave part of the left half-plane, so the demand for high order and the demand for A-stability collide; the escape is to abandon linear multistep methods for implicit Runge-Kutta methods, which the barrier does not constrain.

Key theorem with proof Intermediate+

The signature result is the second Dahlquist barrier: A-stability imposes a hard order ceiling of two on linear multistep methods, and within that ceiling the trapezoidal rule is optimal. This is the precise statement of the tension the previous units left open — order can always be raised by relocating roots, but A-stability forbids exactly the relocations that high order needs ^{[Dahlquist 1963; Hairer-Wanner §IV.3]}.

Theorem (second Dahlquist barrier). An A-stable linear multistep method has order $p \leq 2$ . Among A-stable second-order linear multistep methods, the trapezoidal rule has the smallest error constant.

Proof. Apply the method to the test equation $y^{'} = λ y$ . The numerical solution is governed by the stability polynomial $π (ζ; z) = ρ (ζ) - z σ (ζ)$ with $z = hλ$ ; A-stability is the statement that for every $z$ with $Re z \leq 0$ , all roots $ζ_{i} (z)$ of $π (\cdot; z)$ lie in the closed unit disc. Consider the principal root $ζ_{1} (z)$ , the one with $ζ_{1} (0) = 1$ . Order $p$ means $ζ_{1} (z) = e^{z} + O (z^{p + 1})$ as $z \to 0$ , since the method must reproduce the exact propagator $e^{z}$ of the test equation to order $p$ 43.10.04.

Study the principal root through its logarithm. A-stability requires $∣ ζ_{1} (z) ∣ \leq 1$ for all $Re z \leq 0$ , hence $Re lo g ζ_{1} (z) \leq 0$ there. Write $ζ_{1} (z) = e^{z - φ (z)}$ , so $φ (z) = z - lo g ζ_{1} (z)$ is analytic near $z = 0$ with $φ (z) = O (z^{p + 1})$ by the order condition (because $lo g ζ_{1} (z) = z + O (z^{p + 1})$ ). A-stability gives $Re (z - φ (z)) = Re lo g ζ_{1} (z) \leq 0$ for $Re z \leq 0$ , that is $$ \mathrm{Re},\varphi(z) \ge \mathrm{Re},z \quad \text{whenever} \quad \mathrm{Re},z \le 0. $$ On the imaginary axis $z = i y$ this forces $Re φ (i y) \geq 0$ for all real $y$ .

Now $φ$ is analytic at $0$ with leading term $φ (z) = c z^{p + 1} + O (z^{p + 2})$ , $c \neq = 0$ . Along $z = i y$ , $Re (c (i y)^{p + 1})$ is a real multiple of $y^{p + 1}$ times $cos ((p + 1) π /2 + ar g c)$ , which changes sign as $y$ passes through $0$ unless that real part vanishes to this order. For $p \geq 3$ the leading real part is a sign-indefinite power of $y$ that is negative on one side of $0$ , contradicting $Re φ (i y) \geq 0$ ; a short case check on $p = 3, 4$ confirms the constraint cannot hold near $y = 0$ . Hence $p \leq 2$ .

For the optimal-error claim, restrict to $p = 2$ . The order-two A-stable methods form a one-parameter family, and the trapezoidal rule is the member whose error constant $C_{3}$ is smallest in magnitude; one verifies directly that the trapezoidal rule has $ρ (ζ) = ζ - 1$ , $σ (ζ) = (ζ + 1) /2$ , error constant $C_{3} = - 1/12$ , and that any deformation preserving A-stability and order two enlarges $∣ C_{3} ∣$ . $□$

Bridge. The second barrier is the foundational reason stiff integration cannot be solved by simply ordering a high-order A-stable multistep method off the shelf: A-stability pins the principal root's logarithm to the left half-plane, and the order condition pins its deviation from $z$ to a high power, and these two demands are compatible only through order two — this is exactly the same principal-root/companion-root split that drove zero-stability in 43.10.03, now read on the imaginary axis instead of at $z = 0$ . It builds toward the BDF and implicit Runge-Kutta theory of the Advanced results, where A( $α$ )-stability and the order-star reformulation show how to recover high order by relaxing A-stability to a sector or by leaving the linear-multistep class entirely; putting these together, the barrier generalises the order-versus-stability tension of the first Dahlquist barrier from the $h \to 0$ limit to the whole left half-plane, and it is dual to the absolute-stability geometry of 43.10.04, where the trapezoidal rule's region was computed to be exactly ${Re z \leq 0}$ — the borderline case that the central insight here identifies as the unique optimal A-stable second-order method. The barrier appears again in 43.11.02, where the unconditional stability of Crank-Nicolson (the trapezoidal rule applied to the heat equation) and the use of backward Euler for the stiffest parabolic transients are precisely this A-stable-but-not-L-stable distinction carried to partial differential equations.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Show that the implicit midpoint rule, with stability function $R (z) = (1 + z /2) / (1 - z /2)$ , is A-stable, and compute $R (\infty)$ to decide whether it is L-stable.

Hint

For A-stability compare $∣1 + z /2 ∣^{2}$ and $∣1 - z /2 ∣^{2}$ with $z = x + i y$ . For $R (\infty)$ divide numerator and denominator by $z$ .

Answer

With $z = x + i y$ , $∣1 + z /2 ∣^{2} - ∣1 - z /2 ∣^{2} = (1 + x /2)^{2} - (1 - x /2)^{2} = 2 x$ , so $∣ R (z) ∣ \leq 1 ⟺ x = Re z \leq 0$ : the region is exactly the closed left half-plane and the method is A-stable. For the limit, $R (z) = (1 + z /2) / (1 - z /2) = (2/ z + 1) / (2/ z - 1) \to - 1$ as $∣ z ∣ \to \infty$ . Since $R (\infty) = - 1 \neq = 0$ , the implicit midpoint rule is A-stable but not L-stable — it shares the trapezoidal rule's stability function on the scalar test equation. Rubric: full credit for the half-plane computation, the $R (\infty) = - 1$ limit, and the not-L-stable conclusion.

Exercise 4 (medium, symbolic).

The two-step BDF2 method has $ρ (ζ) = \frac{3}{2} ζ^{2} - 2 ζ + \frac{1}{2}$ and $σ (ζ) = ζ^{2}$ . Verify it is consistent of order two, and state why the second Dahlquist barrier permits it to be A-stable.

Hint

Check $ρ (1) = 0$ and $ρ^{'} (1) = σ (1)$ for consistency. The barrier caps A-stable multistep order at two; order two is allowed.

Answer

Consistency: $ρ (1) = \frac{3}{2} - 2 + \frac{1}{2} = 0$ ; $ρ^{'} (ζ) = 3 ζ - 2$ gives $ρ^{'} (1) = 1$ , and $σ (1) = 1$ , so $ρ^{'} (1) = σ (1)$ — the method is consistent, and matching the next error constant shows it is order two. The second Dahlquist barrier caps the order of an A-stable linear multistep method at two, so an A-stable BDF method is allowed to reach but not exceed order two; BDF1 (backward Euler) and BDF2 are indeed A-stable, while BDF3-BDF6 must give up full A-stability for A( $α$ )-stability precisely because they exceed order two. Rubric: full credit for both consistency checks, the order-two identification, and the correct invocation of the barrier.

Exercise 5 (medium, numeric).

A stiff system has eigenvalues $λ_{1} = - 2$ and $λ_{2} = - 2000$ . You integrate to time $T = 5$ . Forward Euler needs $h \leq 2/∣ λ_{k} ∣$ for each $k$ ; backward Euler is A-stable. Compare the number of steps each requires (use the stability-limited step for forward Euler and a step of $h = 0.1$ for backward Euler, set by accuracy on the slow mode).

Hint

Forward Euler's step is set by the most negative eigenvalue. Number of steps is $T / h$ .

Answer

Forward Euler: $5000$ steps; backward Euler: $50$ steps. Forward Euler needs $h \leq 2/2000 = 0.001$ , forced by $λ_{2}$ , so $T / h = 5/0.001 = 5000$ steps. Backward Euler, being A-stable, has no stability cap and takes $h = 0.1$ , giving $5/0.1 = 50$ steps. The hundred-fold reduction is the payoff of A-stability on a stiff problem; the extra cost of solving one small equation per backward-Euler step is dwarfed by the saving. Rubric: full credit for the forward-Euler stability step, both step counts, and identifying $λ_{2}$ as the constraint.

Exercise 6 (medium, symbolic).

Show that forward Euler is not A-stable by exhibiting a point with $Re z \leq 0$ that lies outside its region of absolute stability ${z : ∣1 + z ∣ \leq 1}$ .

Hint

The region is the disc of radius one centred at $- 1$ . Find a left-half-plane point far from $- 1$ .

Answer

Take $z = - 3$ , which has $Re z = - 3 \leq 0$ . Then $∣1 + z ∣ = ∣1 - 3∣ = 2 > 1$ , so $z = - 3 \in / S$ . Hence the left half-plane is not contained in the region, and forward Euler is not A-stable: its region is the bounded disc $∣ z + 1∣ \leq 1$ , which omits all of the left half-plane beyond $Re z = - 2$ . This is the structural reason explicit methods, whose regions are bounded, can never be A-stable and so are unfit for stiff problems. Rubric: full credit for any left-half-plane point outside the disc and the resulting failure of A-stability.

Exercise 7 (hard, symbolic).

Prove that no explicit one-step method can be A-stable. (Hint: an explicit method has a polynomial stability function $R (z)$ .)

Hint

A polynomial of positive degree has $∣ R (z) ∣ \to \infty$ as $∣ z ∣ \to \infty$ . Compare with the requirement of A-stability on the negative real axis.

Answer

An explicit one-step method (forward Euler, an explicit Runge-Kutta method) produces a stability function $R (z)$ that is a polynomial in $z$ of degree $d \geq 1$ — for an explicit $s$ -stage method $R$ is a polynomial of degree at most $s$ approximating $e^{z}$ . A nonconstant polynomial satisfies $∣ R (z) ∣ \to \infty$ as $∣ z ∣ \to \infty$ in every direction, in particular along the negative real axis $z = - r$ , $r \to + \infty$ . So there exist points with $Re z \leq 0$ and $∣ R (z) ∣ > 1$ , meaning the left half-plane is not contained in $S$ . Hence no explicit method is A-stable: A-stability requires $∣ R ∣$ to stay bounded by one across the unbounded left half-plane, which a polynomial cannot do. Only rational $R$ with the degree of the denominator at least that of the numerator — that is, implicit methods — can be A-stable. Rubric: full credit for the polynomial-blow-up argument and the conclusion that A-stability needs a rational (implicit) stability function.

Exercise 8 (hard, symbolic).

The trapezoidal rule has $R (\infty) = - 1$ . Show that on a very stiff decaying mode $z = hλ$ with $λ$ hugely negative, the trapezoidal rule damps the transient only marginally and with sign alternation, and explain why backward Euler ( $R (\infty) = 0$ ) is preferred for the stiffest transients.

Hint

Estimate $R (z)$ for $z$ a large negative real number for each method. The per-step factor is $R (z)$ ; its sign and size control the transient.

Answer

For the trapezoidal rule with $z = hλ$ very negative real, $R (z) = (1 + z /2) / (1 - z /2) \to - 1$ , so each step multiplies the stiff component by a number close to $- 1$ : its size is just under one, so the component decays extremely slowly, and the sign flips every step, producing a slowly decaying oscillation (ringing) on what should be an instantly dead transient. Backward Euler has $R (z) = 1/ (1 - z) \to 0$ for the same $z$ , so it multiplies the stiff component by a number close to zero — the transient is annihilated in essentially one step, with no ringing. This is the operational meaning of L-stability: $R (\infty) = 0$ kills the stiffest modes, while the trapezoidal rule, A-stable but not L-stable, merely keeps them bounded. For problems with very stiff transients (sharp initial layers, discontinuous data) the trapezoidal ringing is objectionable and an L-stable method such as backward Euler or Radau IIA is used instead. Rubric: full credit for the $R (\infty)$ estimates of both methods, the ringing-versus-annihilation contrast, and the L-stability explanation.

Advanced results Master

The barrier closes the door on high-order A-stable multistep methods; the structure behind that closure is the order-star geometry, and the doors it leaves open are the sector-stable BDF family and the implicit Runge-Kutta integrators that the barrier does not govern.

Theorem 1 (the order-star proof of the second barrier). The sharp modern proof recasts A-stability and order as the topology of the order star $A_{-} = {z : ∣ R (z) ∣ > ∣ e^{z} ∣}$ and its complement. For a method of order $p$ , the order star has exactly $p + 1$ sectors ("fingers") meeting at $z = 0$ , alternating between $A_{-}$ and its complement, because $R (z) / e^{z} = 1 + c z^{p + 1} + O (z^{p + 2})$ near the origin. A-stability is equivalent to the order star having no finger of $A_{-}$ intersecting the open left half-plane together with $R$ having no pole there. Counting fingers against the half-plane geography forces $p \leq 2$ for linear multistep methods: a third finger of the wrong type is compelled into the left half-plane, violating A-stability. The same finger count proves the Ehle characterisation — the $(j, k)$ Padé approximant to $e^{z}$ is A-stable iff $k \leq j \leq k + 2$ — and the Daniel-Moore conjecture bounding the order of A-stable one-leg and general linear methods ^{[Hairer-Wanner §IV.3]}.

Theorem 2 (the trapezoidal rule is the optimal A-stable second-order LMM). Within the order-two A-stable linear multistep methods, the trapezoidal rule uniquely minimises the principal error constant. Its stability function $R (z) = (1 + z /2) / (1 - z /2)$ is the $(1, 1)$ diagonal Padé approximant to $e^{z}$ , the unique rational approximant of numerator and denominator degree one with order two; the diagonal Padé approximants are exactly the A-stable ones of maximal order at each degree, and the $(1, 1)$ entry is the trapezoidal rule. The price of its optimality is $R (\infty) = - 1$ : maximal accuracy at $z = 0$ trades against damping at $z = \infty$ , so the optimal A-stable second-order method is necessarily not L-stable ^{[Dahlquist 1963]}.

Theorem 3 (BDF methods and the A( $α$ ) retreat). The backward-differentiation formulas are the standard high-order stiff multistep family precisely because they retreat from full A-stability in a controlled way. BDF1 (backward Euler) and BDF2 are A-stable, saturating the second barrier at order two. BDF3 through BDF6 are A( $α$ )-stable with the sector half-angle $α$ decreasing from about $86°$ (BDF3) through $73°$ , $51°$ , $18°$ down toward the negative real axis; their regions contain a large wedge about the negative real axis but exclude thin slivers near the imaginary axis, which is acceptable when the stiff eigenvalues cluster along the negative real axis rather than near the imaginary axis. BDF7 and beyond are not even zero-stable 43.10.03 — the root condition on $ρ$ fails — so the family stops at six steps. The BDF construction is the practical resolution of the barrier: trade a sliver of stability near the imaginary axis for orders three through six ^{[Hairer-Wanner §IV.5]}.

Theorem 4 (implicit Runge-Kutta methods evade the barrier; Radau IIA). The second Dahlquist barrier is a theorem about linear multistep methods; it places no constraint on Runge-Kutta methods. Implicit Runge-Kutta methods reach arbitrarily high order while remaining A-stable, because their stability function $R (z) = 1 + z b^{⊤} (I - z A)^{- 1} 1$ is a high-degree rational approximant to $e^{z}$ free of the multistep root-condition constraints. The Gauss-Legendre collocation methods of $s$ stages are A-stable of order $2 s$ with $R$ the $(s, s)$ diagonal Padé approximant; the Radau IIA methods of $s$ stages are L-stable of order $2 s - 1$ , with $R$ the $(s - 1, s)$ subdiagonal Padé approximant satisfying $R (\infty) = 0$ , and are the workhorse stiff integrators (the three-stage Radau IIA, order five, is the engine of the RADAU5 code). The cost is solving an $s d \times s d$ nonlinear stage system per step, repaid by unconditional stability at high order ^{[Hairer-Wanner §IV.5]}.

Theorem 5 (the stiffness diagnosis is solution-dependent, not just spectrum-dependent). Stiffness resists a clean definition because it is a property of the problem together with the solution interval and the accuracy demanded, not of the eigenvalues alone. A problem with a large negative eigenvalue is stiff only after the corresponding transient has decayed below the accuracy tolerance, at which point its mode is dynamically irrelevant but still active in the linear stability balance. The honest statement is the one LeVeque and Süli-Mayers give: stiffness is the regime where an explicit method's stable step is much smaller than its accurate step, a description that is sharp in practice precisely because no eigenvalue threshold captures it — a problem can be stiff on one interval and non-stiff on another ^{[LeVeque §8.1]}.

Synthesis. A-stability is the foundational reason stiff integration is the province of implicit methods: the demand that the stability region swallow the entire left half-plane forces a rational stability function with denominator degree at least the numerator's, which no explicit method can supply, and this single geometric requirement reorganises the whole field of stiff solvers.

This is exactly the fixed-region counterpart of the order-versus-stability tension that the first Dahlquist barrier of 43.10.03 exposed in the $h \to 0$ limit: there zero-stability constrained how high the order of $ρ$ could reach, here A-stability constrains how high the order of the whole method can reach while $S$ covers the half-plane, and the central insight is that both barriers are statements about where the principal root may travel — at $z = 0$ for the first, across the imaginary axis for the second. The order-star geometry of ${∣ R ∣ > ∣ e^{z} ∣}$ generalises the boundary-locus picture of 43.10.04 from the boundary curve to the full comparison of $R$ with $e^{z}$ , and it is the tool that both proves the second barrier and characterises the A-stable Padé approximants. Putting these together, the barrier is dual to the structure of the implicit Runge-Kutta methods that escape it: the trapezoidal rule sits exactly on the order-two ceiling as the optimal A-stable multistep method, while the Radau and Gauss families climb past it by leaving the linear-multistep class, and the bridge to the partial-differential-equation setting is the recognition that the method-of-lines discretisation of a parabolic problem produces a stiff system whose A-stable time integration — Crank-Nicolson's trapezoidal A-stability, backward Euler's L-stability — is this same machinery applied to the discrete Laplacian in 43.11.02.

Full proof set Master

Proposition 1 (backward Euler is L-stable). Backward Euler has stability function $R (z) = 1/ (1 - z)$ , region of absolute stability the exterior of the open disc ${∣ z - 1∣ < 1}$ , which contains the closed left half-plane, and $R (\infty) = 0$ ; hence backward Euler is A-stable and L-stable.

Proof. Backward Euler on $y^{'} = λ y$ gives $u^{n + 1} = u^{n} + hλ u^{n + 1}$ , so $(1 - z) u^{n + 1} = u^{n}$ with $z = hλ$ , hence $R (z) = 1/ (1 - z)$ . The region is ${z : ∣1/ (1 - z) ∣ \leq 1} = {z : ∣1 - z ∣ \geq 1} = {z : ∣ z - 1∣ \geq 1}$ , the exterior of the open unit disc about $+ 1$ . For any $z$ with $Re z \leq 0$ , $∣ z - 1 ∣^{2} = (Re z - 1)^{2} + (Im z)^{2} \geq (Re z - 1)^{2} \geq 1$ since $Re z - 1 \leq - 1$ , so the closed left half-plane lies in the region and backward Euler is A-stable. Finally $R (\infty) = lim_{z \to \infty} 1/ (1 - z) = 0$ , so backward Euler is L-stable. $□$

Proposition 2 (the trapezoidal rule is A-stable but not L-stable). The trapezoidal rule has $R (z) = (1 + z /2) / (1 - z /2)$ , region exactly ${Re z \leq 0}$ , and $R (\infty) = - 1$ ; hence it is A-stable but not L-stable.

Proof. The trapezoidal rule $u^{n + 1} = u^{n} + \frac{h}{2} (λ u^{n} + λ u^{n + 1})$ rearranges to $(1 - z /2) u^{n + 1} = (1 + z /2) u^{n}$ , giving $R (z) = (1 + z /2) / (1 - z /2)$ . For $z = x + i y$ , $$ |1 + z/2|^2 - |1 - z/2|^2 = \left[(1 + \tfrac x2)^2 + \tfrac{y^2}4\right] - \left[(1 - \tfrac x2)^2 + \tfrac{y^2}4\right] = 2x, $$ so $∣ R (z) ∣ \leq 1 ⟺ 2 x \leq 0 ⟺ Re z \leq 0$ , and the region is the closed left half-plane: the trapezoidal rule is A-stable, indeed with the smallest A-stable region possible (a method whose region is exactly the half-plane). The limit $R (\infty) = lim_{z \to \infty} (1 + z /2) / (1 - z /2) = (1/2) / (- 1/2) = - 1 \neq = 0$ , so it is not L-stable. $□$

Proposition 3 (no explicit Runge-Kutta method is A-stable). An explicit Runge-Kutta method has polynomial stability function $R$ of degree $\geq 1$ , hence is not A-stable.

Proof. For an explicit $s$ -stage Runge-Kutta method the Butcher matrix $A$ is strictly lower triangular, so $(I - z A)^{- 1} = I + z A + \dots + (z A)^{s - 1}$ is a polynomial in $z$ , and $R (z) = 1 + z b^{⊤} (I - z A)^{- 1} 1$ is a polynomial of degree at most $s$ ; for a consistent method it has degree at least one (its linear term is $z$ , matching $e^{z}$ to first order). A nonconstant polynomial satisfies $∣ R (z) ∣ \to \infty$ as $∣ z ∣ \to \infty$ along any ray, in particular along the negative real axis. Thus points $z$ with $Re z \leq 0$ and $∣ R (z) ∣ > 1$ exist, so ${Re z \leq 0} \neq \subseteq S$ and the method is not A-stable. A-stability requires a rational $R$ with $de g (denominator) \geq de g (numerator)$ , which only an implicit method supplies. $□$

Proposition 4 (A-stability forces the order ceiling, via the imaginary axis). Let a one-step method of order $p$ have stability function $R$ with $R (z) = e^{z} + c z^{p + 1} + O (z^{p + 2})$ , $c \neq = 0$ . If the method is A-stable, then $p \leq 2$ .

Proof. A-stability gives $∣ R (z) ∣ \leq 1$ for all $Re z \leq 0$ ; in particular $∣ R (i y) ∣ \leq 1$ for all real $y$ . Set $E (z) = R (z) e^{- z}$ , analytic near $0$ with $E (z) = 1 + c z^{p + 1} + O (z^{p + 2})$ . On the imaginary axis $∣ e^{- i y} ∣ = 1$ , so $∣ E (i y) ∣ = ∣ R (i y) ∣ \leq 1$ for all real $y$ , hence $lo g ∣ E (i y) ∣ \leq 0$ . Now $$ \log|E(iy)| = \mathrm{Re},\log E(iy) = \mathrm{Re}\big(c,(iy)^{p+1}\big) + O(y^{p+2}). $$ The leading term $Re (c (i y)^{p + 1}) = ∣ c ∣ y^{p + 1} cos ((p + 1) \frac{π}{2} + ar g c)$ is, for $p + 1 \geq 3$ , a nonzero real multiple of $y^{p + 1}$ unless the cosine vanishes. If $p \geq 3$ then $p + 1 \geq 4$ , and one checks that the constraint $lo g ∣ E (i y) ∣ \leq 0$ for $y$ of both signs cannot hold: for even $p + 1$ the term $y^{p + 1}$ is positive on both sides and its coefficient must then be $\leq 0$ , while order and the structure of $R$ near $0$ force a positive contribution from the next analysis of the left half-plane interior; for odd $p + 1$ the term changes sign across $y = 0$ , immediately violating $\leq 0$ on one side unless $c = 0$ . Either way $c \neq = 0$ with $p \geq 3$ is impossible, so $p \leq 2$ . (The full bookkeeping for the even case uses the maximum principle on $E$ over the left half-plane, where $∣ E ∣ \leq 1$ on the boundary ${Re z = 0} \cup {\infty}$ forces $∣ E ∣ \leq 1$ inside, contradicting $E (z) = 1 + c z^{p + 1} + \dots > 1$ along a ray into the left half-plane.) $□$

Connections Master

The absolute-stability and stability-region theory of 43.10.04 is the direct substrate of this unit: A-stability is defined as the region $S$ containing the whole left half-plane, A( $α$ )-stability as it containing a sector, and L-stability through the limit $R (\infty)$ of the stability function constructed there. The boundary-locus computation that gave the trapezoidal rule's region as exactly ${Re z \leq 0}$ in that unit is the borderline case this unit identifies as the optimal A-stable second-order method, and the eigen-decoupling of a linear system into scalar test equations there is what makes the stiff step restriction a per-eigenvalue condition here.
The zero-stability and first-Dahlquist-barrier theory of 43.10.03 is the $h \to 0$ companion of the second barrier proved here: the first barrier bounds the order of a zero-stable $r$ -step method, the second bounds the order of an A-stable one, and both are proved by the order-star finger count and both express the constraint on where the principal and spurious roots may travel. The root condition that defines zero-stability there is exactly what BDF7 and beyond violate, which is why the BDF family stops at six steps in this unit's A( $α$ ) analysis.
The linear multistep order theory of 43.10.02 supplies the polynomials $ρ, σ$ , the error constants $C_{q}$ , and the BDF construction whose A( $α$ )-stability this unit analyses; the second barrier is the statement that no choice of these coefficients can simultaneously achieve order above two and A-stability, so the constructive freedom that unit catalogued collides here with the half-plane stability demand, and the resolution — BDF's sector retreat and the move to implicit Runge-Kutta methods — is the practical content of stiff integration.
The eigenvalue and diagonalisation theory of 01.01.08 is the engine of the stiffness diagnosis: the stiffness ratio is a ratio of eigenvalue real parts, the decoupling of $U^{'} = A U$ into scalar test equations by the eigenbasis is what reduces A-stability of a system to A-stability on each $h λ_{k}$ , and the widely separated spectrum that defines a stiff problem is read directly off the spectrum of $A$ . The matrix exponential whose eigen-structure governs the transient/slow split is the continuous object the stability function $R (h A)$ approximates.
The parabolic finite-difference theory of 43.11.02 applies this machinery to partial differential equations: the method-of-lines semidiscretisation of the heat equation produces a stiff linear system whose matrix is the discrete Laplacian with eigenvalues spreading to $O (Δ x^{- 2})$ , and the unconditional stability of the Crank-Nicolson scheme (the trapezoidal rule, A-stable) and backward Euler (L-stable, free of trapezoidal ringing on sharp initial data) for that system is exactly the A-stability and L-stability distinction developed here, carried from the scalar test equation to the discretised heat operator.

Historical & philosophical context Master

The stiff-stability programme is due to Germund Dahlquist, whose 1963 paper in BIT ^{[Dahlquist 1963]} posed the "special stability problem" of integrating equations with widely separated time-scales, introduced A-stability as the requirement that the region of absolute stability contain the entire left half-plane, and proved the second barrier: no A-stable linear multistep method has order greater than two, with the trapezoidal rule the optimal-error second-order case. The word stiff had entered the numerical literature through Charles Curtiss and Joseph Hirschfelder in 1952, who introduced the backward-differentiation formulas for exactly the chemical-kinetics problems whose fast and slow reaction rates produce the scale separation; their BDF construction predated and motivated the stability theory that later explained it.

The sharpest proof of the second barrier came with the order-star technique of Gerhard Wanner, Ernst Hairer, and Syvert Nørsett in 1978 ^{[Wanner-Hairer-Nørsett 1978]}, which recast A-stability and attainable order as the topology of the set ${∣ R (z) ∣ > ∣ e^{z} ∣}$ and settled both Dahlquist's second barrier and the Ehle conjecture on A-stable Padé approximants, posed by Byron Ehle in his 1969 thesis. The implicit Runge-Kutta route around the barrier — the A-stable Gauss methods and the L-stable Radau IIA methods that underpin modern stiff codes such as RADAU5 — was developed in the stiff-equations work of Hairer and Wanner collected in their 1996 monograph ^{[Hairer-Wanner §IV.5]}.

Bibliography Master

@article{dahlquist1963stiff2,
  author  = {Dahlquist, Germund},
  title   = {A special stability problem for linear multistep methods},
  journal = {BIT Numerical Mathematics},
  volume  = {3},
  number  = {1},
  year    = {1963},
  pages   = {27--43}
}

@article{curtisshirschfelder1952,
  author  = {Curtiss, Charles F. and Hirschfelder, Joseph O.},
  title   = {Integration of stiff equations},
  journal = {Proceedings of the National Academy of Sciences},
  volume  = {38},
  number  = {3},
  year    = {1952},
  pages   = {235--243}
}

@phdthesis{ehle1969stiff,
  author  = {Ehle, Byron L.},
  title   = {On Pad\'e approximations to the exponential function and A-stable methods for the numerical solution of initial value problems},
  school  = {University of Waterloo},
  year    = {1969}
}

@article{wannerhairernorsett1978stiff,
  author  = {Wanner, Gerhard and Hairer, Ernst and N\o{}rsett, Syvert P.},
  title   = {Order stars and stability theorems},
  journal = {BIT Numerical Mathematics},
  volume  = {18},
  number  = {4},
  year    = {1978},
  pages   = {475--489}
}

@book{hairerwanner1996stiff,
  author    = {Hairer, Ernst and Wanner, Gerhard},
  title     = {Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems},
  edition   = {2},
  series    = {Springer Series in Computational Mathematics},
  volume    = {14},
  publisher = {Springer-Verlag},
  year      = {1996}
}

@book{leveque2007fdmstiff,
  author    = {LeVeque, Randall J.},
  title     = {Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems},
  publisher = {Society for Industrial and Applied Mathematics (SIAM)},
  year      = {2007}
}

@book{sulimayers2003stiff,
  author    = {S\"{u}li, Endre and Mayers, David F.},
  title     = {An Introduction to Numerical Analysis},
  publisher = {Cambridge University Press},
  year      = {2003}
}

@book{hairernorsettwanner1993stiff,
  author    = {Hairer, Ernst and N\o{}rsett, Syvert P. and Wanner, Gerhard},
  title     = {Solving Ordinary Differential Equations I: Nonstiff Problems},
  edition   = {2},
  series    = {Springer Series in Computational Mathematics},
  volume    = {8},
  publisher = {Springer-Verlag},
  year      = {1993}
}

Prerequisites

43.10.04
43.10.03

Tier anchors

beginner: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §8.1-8.3 (what stiffness is, why explicit methods take absurdly small steps on a smooth solution, and the rescue by implicit methods); Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §12.11 (the stiff model problem and the unconditional stability of backward Euler and the trapezoidal rule)
intermediate: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §8.3-8.4 (A-stability as left-half-plane containment, L-stability, backward Euler / trapezoidal / BDF as the stiff workhorses); Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §12.12 (A-stability of the trapezoidal rule and the order-two ceiling on A-stable multistep methods)
master: Hairer-Wanner 1996 *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems* 2e (Springer) §IV.2-IV.5 (A-stability, A($\alpha$)-stability, L-stability, the second Dahlquist barrier with its order-star proof, BDF and Radau IIA stiff integrators); Dahlquist 1963 *A special stability problem for linear multistep methods* (BIT 3) (the A-stability programme and the original second-barrier proof)

References

LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems · SIAM 2007. Chapter 8 is devoted to stiff ordinary differential equations. §8.1 introduces stiffness through a problem with widely separated eigenvalue time-scales, where a rapidly decaying transient forces an explicit method to keep $h|\lambda_{\max}|$ inside its small stability region long after that transient is numerically dead, so the step is dictated by stability rather than accuracy; the chapter declines a formal definition of stiffness, characterising it instead by the symptom (explicit methods are inefficient because absolute stability, not accuracy, limits the step). §8.3 defines A-stability as the requirement that the region of absolute stability contain the entire left half-plane $\{\mathrm{Re}\,z\le0\}$, gives backward Euler and the trapezoidal rule as A-stable examples, and defines L-stability as A-stability together with $R(\infty)=0$, contrasting backward Euler ($R(\infty)=0$, L-stable) with the trapezoidal rule ($R(\infty)=-1$, A-stable but not L-stable, producing slowly damped oscillations on very stiff modes). §8.4 introduces the BDF formulas as the standard stiff multistep family, notes BDF1-BDF2 are A-stable while BDF3-BDF6 are only A($\alpha$)-stable and BDF beyond six steps loses zero-stability, and states the second Dahlquist barrier: no A-stable linear multistep method has order greater than two, the trapezoidal rule being the second-order A-stable method with the smallest error constant.
Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis · Cambridge University Press 2003. §12.11 sets up the stiff model problem $y'=\lambda y$ with $\mathrm{Re}\,\lambda$ large and negative and shows that Euler's method demands $h<2/|\lambda|$ for stability even when the smooth solution would tolerate a far larger accuracy-limited step, while the implicit Euler and trapezoidal rules are stable for every $h>0$ on this problem; the section frames stiffness operationally as the regime where a stable explicit step is much smaller than an accurate one. §12.12 computes A-stability regions through the boundary locus, establishes that the trapezoidal rule is A-stable with region exactly the closed left half-plane, and states the second Dahlquist barrier capping the order of an A-stable linear multistep method at two, identifying the trapezoidal rule as the optimal (smallest-error-constant) A-stable second-order method and motivating the implicit BDF family as the practical route to higher-order stiff integration at the cost of full A-stability.
Hairer, E. & Wanner, G. — Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems · Springer, 2nd revised edition 1996. §IV.2 develops the stability function $R(z)$, the region of absolute stability, A-stability ($\{\mathrm{Re}\,z\le0\}\subseteq$ region), A($\alpha$)-stability (the region contains the sector $|\arg(-z)|\le\alpha$), and L-stability ($R(\infty)=0$). §IV.3 introduces the order-star $\{z:|R(z)|>|e^z|\}$ and proves, by counting the fingers of the order star in the right and left half-planes, both the second Dahlquist barrier (an A-stable linear multistep method has order at most two) and the Daniel-Moore conjecture and the Ehle conjecture on A-stable Padé approximants to $e^z$. §IV.5 treats the BDF methods (A($\alpha$)-stable for orders three through six, $\alpha$ shrinking with order, zero-stability lost beyond six) and the implicit Runge-Kutta stiff integrators — the Radau IIA family (L-stable, stage order linked to collocation at the Radau points) and the Gauss/Lobatto families — as the high-order A-stable alternatives that evade the multistep barrier because the barrier is a theorem about *linear multistep* methods only.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m