43.10.04 · numerical-analysis / numerical-odes

Absolute stability, stability regions, and the linear test equation

shipped3 tiersLean: none

Anchor (Master): Hairer-Nørsett-Wanner 1993 *Solving Ordinary Differential Equations I: Nonstiff Problems* 2e (Springer) §III.3 and Hairer-Wanner 1996 *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems* 2e (Springer) §IV.2-IV.3 (stability functions, stability regions, the boundary locus, order stars and the Ehle conjecture); Dahlquist 1963 *A special stability problem for linear multistep methods* (BIT 3) (the test-equation stability framework and the A-stability programme)

Intuition Beginner

Two questions decide whether a step-by-step method can be trusted. The first is whether it tracks the true curve well when you shrink the step toward zero — that is accuracy, and the earlier units handled it. The second is different and more practical: for the step size you actually plan to use, does a small error fade away or build up? That second question is absolute stability, and it is the one that tells you how big a step you are allowed to take before the computed answer turns to garbage.

The trick is to test every method on the same toy problem. Pick the equation whose slope is a fixed multiple of the current height: the height either decays smoothly toward zero or wobbles, depending on that multiple. The true answer fades when the multiple pulls downward. A good method, handed this fading problem, should also produce numbers that fade. A poor method, or a step that is too large, will instead produce numbers that grow — the exact opposite of the truth.

When you run one step of a method on this toy problem, the method multiplies the current value by a single number. Call that number the step's amplification factor. If its size is at most one, each step shrinks or holds the value, so errors fade and the method is stable. If its size is bigger than one, every step inflates the value and errors snowball. The amplification factor depends on the step size and the problem's multiplier only through their product, a single quantity usually written $z$ .

Because $z$ packs the step and the problem together, you can draw a map. Mark every value of $z$ for which the amplification factor has size at most one. That shaded set of points is the method's region of absolute stability. Choosing a step size for a real decaying problem becomes a geometry question: compute $z$ for your problem and step, and check whether it lands inside the shaded region. Forward Euler has a small region, so it forces tiny steps on fast-decaying problems; backward Euler has a huge region, so it stays calm with large steps. The shape of that region, not the method's accuracy, is what caps your step size.

Visual Beginner

Picture the plane of all possible values of the combined quantity $z$ . The horizontal direction measures decay or growth; the vertical direction measures wobble. Each method shades the patch of this plane where its amplification factor has size at most one. Land your problem's $z$ inside the patch and errors fade; land outside and they grow.

The table below shows where four standard methods place their shaded patch and what that means for the step size you may take on a fast-decaying problem. The pattern is the whole lesson: explicit methods carry small patches and force small steps, while implicit methods carry large patches and tolerate large steps.

method	shaded patch (its region)	step on a fast-decaying problem
forward Euler	small disc near the origin	forced very small
classical RK4	rounded bulb, bigger than Euler's disc	small but a few times larger
backward Euler	almost the whole plane	essentially unrestricted
trapezoidal rule	exactly the left half-plane	unrestricted for decay

Each row is a method's licence: the bigger the shaded patch, the larger the step it forgives before errors start to grow.

Worked example Beginner

Let us measure the largest stable step for forward Euler on a fast-decaying problem and watch what happens just past the edge. The toy problem has slope equal to $- 100$ times the height, so the true answer decays quickly toward zero. We start at height $1$ .

Forward Euler updates by taking one step of size $h$ along the current slope: the new height is the old height plus $h$ times the slope, which here is the old height times the factor $1 + h (- 100) = 1 - 100 h$ . So each step multiplies the height by the number $1 - 100 h$ . That number is the amplification factor, and stability needs its size to be at most $1$ .

Step 1. Solve for the boundary. The size of $1 - 100 h$ is at most $1$ exactly when $1 - 100 h$ sits between $- 1$ and $1$ . The upper end gives $1 - 100 h \leq 1$ , always true for positive $h$ . The lower end gives $1 - 100 h \geq - 1$ , that is $100 h \leq 2$ , so $h \leq 0.02$ . The largest stable step is $h = 0.02$ .

Step 2. Take a safe step, $h = 0.01$ . The factor is $1 - 100 (0.01) = 1 - 1 = 0$ . Starting at $1$ , the next value is $1 \times 0 = 0$ , then it stays at $0$ . The numbers fade, matching the true decay.

Step 3. Take an unsafe step, $h = 0.03$ . The factor is $1 - 100 (0.03) = 1 - 3 = - 2$ , whose size is $2$ . Starting at $1$ , the heights run $1, - 2, 4, - 8, 16$ , flipping sign and doubling in size. The computed answer explodes while the true answer quietly decays to zero.

Step 4. Read off the geometry. The combined quantity is $z = h \times (- 100)$ . The safe step $h = 0.01$ gives $z = - 1$ , inside forward Euler's region. The unsafe step $h = 0.03$ gives $z = - 3$ , outside it. The edge $h = 0.02$ gives $z = - 2$ , the left tip of the region.

What this tells us: accuracy did not fail here — the method is perfectly accurate for tiny steps. Stability failed. The step size was too large for forward Euler's small region, and that alone turned a decaying answer into a blow-up. A method with a larger region, such as backward Euler, would have stayed stable at every one of these steps.

Check your understanding Beginner

Exercise (easy, multiple choice).

Backward Euler is run on a decaying problem with a very large step. Its amplification factor for the combined quantity $z$ is $1/ (1 - z)$ , and the decaying problem has $z$ far out on the negative axis. What happens to the computed values?

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

Hint

For $z$ a large negative number, compute the size of $1/ (1 - z)$ . Is it bigger or smaller than one?

Answer

B. They fade toward zero. For a large negative $z$ , the denominator $1 - z$ is large, so $1/ (1 - z)$ has size well below one and each step shrinks the value. Feedback-correct: backward Euler is stable for every decaying problem at any step, which is the practical reason it survives where forward Euler blows up. Feedback-wrong: large steps do not destabilise backward Euler, because its region covers the entire decaying side of the plane.

Formal definition Intermediate+

Fix the Dahlquist linear test equation $$ y'(t) = \lambda y(t), \qquad \lambda \in \mathbb{C}, \qquad y(0) = y_0, $$ whose exact solution $y (t) = y_{0} e^{λ t}$ decays when $Re λ < 0$ and oscillates without growth when $Re λ = 0$ . The complex multiplier $λ$ models one eigenvalue of a linearised system; the plane in which $λ$ lives, and the rescaled variable $z = hλ \in C$ , are the holomorphic setting of 06.01.01.

Definition (stability function of a one-step method). A one-step method applied to the test equation produces the recurrence $u^{n + 1} = R (z) u^{n}$ with $z = hλ$ , where the rational or polynomial function $R$ is the method's stability function. For the basic methods, $$ R_{\text{FE}}(z) = 1 + z, \quad R_{\text{BE}}(z) = \frac{1}{1 - z}, \quad R_{\text{TR}}(z) = \frac{1 + z/2}{1 - z/2}, \quad R_{\text{RK4}}(z) = 1 + z + \tfrac{z^2}{2} + \tfrac{z^3}{6} + \tfrac{z^4}{24}, $$ for forward Euler, backward Euler, the trapezoidal rule, and classical fourth-order Runge-Kutta. For a Runge-Kutta method with Butcher tableau $(A, b)$ and $1 = (1, \dots, 1)^{⊤}$ , $R (z) = 1 + z b^{⊤} (I - z A)^{- 1} 1$ , a rational function whose Taylor expansion matches $e^{z}$ to the order of the method.

Definition (region of absolute stability). The region of absolute stability of a one-step method is $$ \mathcal{S} = {, z \in \mathbb{C} : |R(z)| \le 1 ,}, $$ the set of $z = hλ$ for which the test-equation iteration neither grows. The method is absolutely stable at a given $z$ when $z \in S$ . Its intersection with the negative real axis is the interval of absolute stability.

Definition (characteristic polynomial and region for a linear multistep method). For an $r$ -step linear multistep method 43.10.02 with first and second characteristic polynomials $ρ (ζ) = \sum_{j} α_{j} ζ^{j}$ and $σ (ζ) = \sum_{j} β_{j} ζ^{j}$ , applying the method to $y^{'} = λ y$ gives the homogeneous recurrence with the stability polynomial $$ \pi(\zeta; z) = \rho(\zeta) - z,\sigma(\zeta), \qquad z = h\lambda. $$ The region of absolute stability is $$ \mathcal{S} = {, z \in \mathbb{C} : \text{every root } \zeta_i(z) \text{ of } \pi(\cdot,; z) \text{ has } |\zeta_i| \le 1,\ \text{with boundary roots simple} ,}. $$ A single-root one-step method recovers the previous definition, with $π (ζ; z) = ζ - R (z)$ and the lone root $ζ = R (z)$ lying in the closed unit disc iff $∣ R (z) ∣ \leq 1$ . The notation $z = hλ$ , the stability function $R (z)$ , the stability polynomial $π (ζ; z)$ , the region $S$ , and the boundary-locus curve $z (θ)$ are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

"Absolute stability is the same as zero-stability." They are different limits. Zero-stability 43.10.03 is the $z = 0$ statement that the roots of $ρ$ alone satisfy the root condition as $h \to 0$ . Absolute stability fixes $h$ and asks the root condition of $π (ζ; z) = ρ (ζ) - z σ (ζ)$ for a particular $z = hλ \neq = 0$ . A zero-stable method can have an empty interior or a small region; leapfrog is zero-stable yet has no region of absolute stability in the open left half-plane.
"Inside the region means accurate." The region only controls whether errors decay; it says nothing about accuracy. A step well inside $S$ can still be far too large to resolve the solution. Stability is necessary, not sufficient, for a good answer.
"The boundary locus $z (θ) = ρ (e^{i θ}) / σ (e^{i θ})$ is always the boundary of the region." The locus is the set where some root sits exactly on the unit circle, so $\partial S \subseteq z ([0, 2 π))$ ; but the curve can self-intersect and enclose several lobes, and a winding-number/root count is needed to decide which side is the stable region. The locus contains the boundary; it is not literally equal to it where the curve crosses itself.

Key theorem with proof Intermediate+

The signature result is the boundary-locus characterisation: the boundary of a linear multistep method's region of absolute stability is traced by the simple formula $z = ρ (ζ) / σ (ζ)$ as $ζ$ runs over the unit circle. This converts an implicit condition on all roots of a $z$ -dependent polynomial into an explicit parametric curve one can plot directly, and it is the computational heart of the whole stability-region machinery ^{[LeVeque §7.3; Süli-Mayers §12.12]}.

Theorem (boundary locus). Let an $r$ -step linear multistep method have characteristic polynomials $ρ, σ$ with $σ (ζ) \neq = 0$ on $∣ ζ ∣ = 1$ , and stability polynomial $π (ζ; z) = ρ (ζ) - z σ (ζ)$ . A point $z_{0} \in C$ lies on the boundary $\partial S$ of the region of absolute stability only if some root of $π (\cdot; z_{0})$ has modulus exactly one; equivalently, $\partial S$ is contained in the image of the unit circle under the boundary-locus map $$ z(\theta) = \frac{\rho(e^{i\theta})}{\sigma(e^{i\theta})}, \qquad \theta \in [0, 2\pi). $$ For forward Euler ( $ρ = ζ - 1$ , $σ = 1$ ) this is the circle $z (θ) = e^{i θ} - 1$ of radius one centred at $- 1$ ; for backward Euler ( $σ = ζ$ ) it is $z (θ) = 1 - e^{- i θ}$ , the circle of radius one centred at $+ 1$ , whose exterior together with that disc's complement is the (large) region.

Proof. Membership $z \in S$ is governed by the location of the roots $ζ_{i} (z)$ of $π (\cdot; z)$ relative to the closed unit disc. As $z$ varies, the roots move continuously: $π (ζ; z)$ is a polynomial in $ζ$ whose coefficients depend continuously (indeed analytically) on $z$ , and the roots of a polynomial depend continuously on its coefficients. The number of roots of $π (\cdot; z)$ strictly inside the open unit disc is therefore a locally constant integer-valued function of $z$ on any open set where no root sits on the circle $∣ ζ ∣ = 1$ . Consequently this count can change only across values of $z$ at which some root crosses $∣ ζ ∣ = 1$ . The region $S$ is, up to the boundary-simplicity stipulation, a union of connected components on which the inside-count is maximal, so its topological boundary $\partial S$ is contained in the locus where a root lies on the unit circle.

A root of modulus one is a value $ζ = e^{i θ}$ with $π (e^{i θ}; z) = 0$ , that is $ρ (e^{i θ}) - z σ (e^{i θ}) = 0$ . Since $σ (e^{i θ}) \neq = 0$ by hypothesis, solve for $z$ : $$ z = \frac{\rho(e^{i\theta})}{\sigma(e^{i\theta})} = z(\theta). $$ Hence every boundary point of $S$ has the form $z (θ)$ for some $θ$ , which is the claimed inclusion $\partial S \subseteq z ([0, 2 π))$ . For forward Euler, $z (θ) = (e^{i θ} - 1) /1$ parametrises the unit circle centred at $- 1$ ; for backward Euler, $z (θ) = (e^{i θ} - 1) / e^{i θ} = 1 - e^{- i θ}$ parametrises the unit circle centred at $+ 1$ . To finish identifying which side is $S$ , evaluate $∣ R (z) ∣$ at one interior test point: for forward Euler $R (z) = 1 + z$ , and $z = - 1$ (the centre) gives $∣ R (- 1) ∣ = 0 \leq 1$ , so the disc about $- 1$ is the region; for backward Euler $R (z) = 1/ (1 - z)$ , and $z = - 1$ gives $∣ R (- 1) ∣ = 1/2 \leq 1$ , placing the entire exterior of the disc about $+ 1$ inside $S$ . $□$

Bridge. The boundary locus is the foundational reason a method's stability region is computable at all: the implicit demand that all roots of $π (\cdot; z)$ stay in the disc collapses, on the boundary, to the explicit one-line map $z = ρ / σ$ evaluated on the unit circle, and this is exactly the same companion-matrix root-tracking used for zero-stability in 43.10.03, now run at fixed $z \neq = 0$ rather than at $z = 0$ . This builds toward the A-stability theory of 43.10.05, where the question becomes whether the region contains the entire left half-plane, and the trapezoidal rule's region — proved below to be exactly ${Re z \leq 0}$ — is the borderline case that putting these together with the order theory yields the Dahlquist second barrier. The construction generalises the scalar one-step picture, since for a single-root method $z (θ)$ degenerates to the level curve $∣ R (z) ∣ = 1$ ; it is dual to the eigenvalue decoupling of 01.01.08, where each eigenvalue $λ_{k}$ of a linear system contributes one test equation and the step is admissible exactly when every $h λ_{k}$ lands in $S$ . The locus appears again in 43.11.02, where the eigenvalues of the discrete Laplacian are fed into a stability region to fix the parabolic step restriction, and the central insight is that one parametric curve, the image of the unit circle under $ρ / σ$ , encodes the entire fixed-step stability budget of the method.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Derive the boundary locus of the backward Euler method ( $ρ (ζ) = ζ - 1$ , $σ (ζ) = ζ$ ) and identify the region of absolute stability.

Hint

Compute $z (θ) = ρ (e^{i θ}) / σ (e^{i θ}) = (e^{i θ} - 1) / e^{i θ}$ , then simplify and test an interior point.

Answer

$z (θ) = (e^{i θ} - 1) / e^{i θ} = 1 - e^{- i θ}$ , which traces the circle of radius $1$ centred at $z = 1$ as $θ$ runs over $[0, 2 π)$ . The plane is split into the open disc $∣ z - 1∣ < 1$ and its exterior. Testing $z = - 1$ (exterior): $R (z) = 1/ (1 - z)$ gives $∣ R (- 1) ∣ = 1/2 \leq 1$ , so the exterior of the disc about $+ 1$ is the region of absolute stability. It contains the whole left half-plane and most of the right half-plane, so backward Euler is A-stable and unconditionally stable on decaying problems. Rubric: full credit for the simplified locus, the circle identification, and the correct stable side via an interior test.

Exercise 4 (medium, symbolic).

The two-step Adams-Bashforth method has $ρ (ζ) = ζ^{2} - ζ$ and $σ (ζ) = \frac{3}{2} ζ - \frac{1}{2}$ . Write the boundary-locus map $z (θ)$ and evaluate it at $θ = π$ to find where the locus meets the negative real axis.

Hint

$z (θ) = (e^{2 i θ} - e^{i θ}) / (\frac{3}{2} e^{i θ} - \frac{1}{2})$ . At $θ = π$ , $e^{i θ} = - 1$ and $e^{2 i θ} = 1$ .

Answer

$z (θ) = \frac{e ^{2 i θ} - e ^{i θ}}{\frac{3}{2} e ^{i θ} - \frac{1}{2}}$ . At $θ = π$ : numerator $= 1 - (- 1) = 2$ , denominator $= \frac{3}{2} (- 1) - \frac{1}{2} = - 2$ , so $z (π) = 2/ (- 2) = - 1$ . At $θ = 0$ , $z (0) = (1 - 1) / (\frac{3}{2} - \frac{1}{2}) = 0$ , the principal-root point. So the AB2 boundary locus on the negative real axis runs from $0$ to $- 1$ , and the interval of absolute stability is $(- 1, 0)$ — a smaller interval than forward Euler's $(- 2, 0)$ , reflecting the price of higher order in an explicit method. Rubric: full credit for the map, $z (π) = - 1$ , and the resulting interval.

Exercise 5 (medium, numeric).

A linear system $U^{'} = A U$ has eigenvalues $λ_{1} = - 1$ and $λ_{2} = - 1000$ . You integrate with forward Euler. What is the largest step $h$ that keeps the computation absolutely stable, and which eigenvalue forces it?

Hint

Each eigenvalue must satisfy $h λ_{k}$ in forward Euler's interval $[- 2, 0]$ . The binding constraint is the most negative eigenvalue.

Answer

$h = 0.002$ , forced by $λ_{2} = - 1000$ . Stability needs $h λ_{k} \geq - 2$ for each $k$ , that is $h \leq 2/∣ λ_{k} ∣$ . For $λ_{1}$ : $h \leq 2$ . For $λ_{2}$ : $h \leq 2/1000 = 0.002$ . The smaller bound wins, so $h = 0.002$ . The fast-decaying mode $λ_{2}$ dictates a step a thousand times smaller than accuracy on the slow mode would need — the signature of stiffness 43.10.05, where the stability region, not accuracy, caps the step. Rubric: full credit for the per-eigenvalue bounds and identifying $λ_{2}$ as binding.

Exercise 6 (medium, symbolic).

Show that the classical RK4 stability function $R (z) = 1 + z + z^{2} /2 + z^{3} /6 + z^{4} /24$ satisfies $R (z) = e^{z} + O (z^{5})$ , and explain why this forces the RK4 interval of absolute stability to include a neighbourhood of $0$ on the negative real axis.

Hint

Expand $e^{z} = \sum_{k \geq 0} z^{k} / k!$ and compare termwise through degree $4$ . Then note $R (0) = 1$ and consider $∣ R (z) ∣$ for small negative real $z$ .

Answer

The Taylor series $e^{z} = 1 + z + z^{2} /2 + z^{3} /6 + z^{4} /24 + z^{5} /120 + \dots$ agrees with $R (z)$ through the $z^{4}$ term, so $R (z) - e^{z} = - z^{5} /120 + O (z^{6}) = O (z^{5})$ ; the agreement to order $4$ is the order of RK4. On the negative real axis, $e^{z} < 1$ for $z < 0$ , and since $R (z)$ approximates $e^{z}$ with $R (0) = 1$ and $R^{'} (0) = 1$ , for small $z < 0$ we have $R (z) \approx e^{z} < 1$ , so $∣ R (z) ∣ < 1$ there. Hence a neighbourhood of $0$ on the negative axis lies in the region; numerically the RK4 real interval extends to about $z \approx - 2.79$ . Rubric: full credit for the order-4 match, the $O (z^{5})$ remainder, and the $e^{z} < 1$ argument for small negative $z$ .

Exercise 7 (hard, symbolic).

Prove that the trapezoidal rule's region of absolute stability is exactly the closed left half-plane ${z : Re z \leq 0}$ , and show that its stability function satisfies $R (z) \to - 1$ as $z \to - \infty$ along the real axis. Comment on what the second fact implies for very stiff decaying modes.

Hint

Use $R (z) = (1 + z /2) / (1 - z /2)$ . For the region, compare $∣1 + z /2 ∣^{2}$ with $∣1 - z /2 ∣^{2}$ . For the limit, divide numerator and denominator by $z$ .

Answer

Write $z = x + i y$ . Then $∣1 + z /2 ∣^{2} - ∣1 - z /2 ∣^{2} = [(1 + x /2)^{2} + (y /2)^{2}] - [(1 - x /2)^{2} + (y /2)^{2}] = 2 x$ . So $∣ R (z) ∣ \leq 1 ⟺ ∣1 + z /2 ∣^{2} \leq ∣1 - z /2 ∣^{2} ⟺ 2 x \leq 0 ⟺ Re z \leq 0$ , with equality $∣ R ∣ = 1$ exactly on the imaginary axis. Hence $S = {Re z \leq 0}$ , the closed left half-plane: the trapezoidal rule is A-stable. For the limit, $R (z) = (1 + z /2) / (1 - z /2) = (2/ z + 1) / (2/ z - 1) \to 1/ (- 1) = - 1$ as $∣ z ∣ \to \infty$ . So $R (\infty) = - 1$ , not $0$ : on a very stiff mode ( $z = hλ$ with $λ$ hugely negative) the amplification factor approaches $- 1$ in size, so the fast transient is damped only marginally and oscillates in sign rather than being annihilated. This is why the trapezoidal rule is A-stable but not L-stable, and why backward Euler ( $R (\infty) = 0$ ) is sometimes preferred for the very stiffest modes 43.10.05. Rubric: full credit for the half-plane proof, the $R (\infty) = - 1$ limit, and the stiff-mode interpretation.

Exercise 8 (hard, symbolic).

For a constant-coefficient linear system $U^{'} = A U$ with $A$ diagonalisable, $A = V Λ V^{- 1}$ , show that applying a one-step method with stability function $R$ gives $U^{n} = V R (h Λ) V^{- 1} U^{0}$ , and deduce that the numerical solution stays bounded for all $n$ if and only if every eigenvalue $λ_{k}$ of $A$ satisfies $h λ_{k} \in S$ .

Hint

In the eigenbasis the system decouples into scalar test equations $w_{k}^{'} = λ_{k} w_{k}$ . The method acts on each as $w_{k}^{n + 1} = R (h λ_{k}) w_{k}^{n}$ . Reassemble.

Answer

Set $W^{n} = V^{- 1} U^{n}$ , the coordinates in the eigenbasis of $A$ 01.01.08. Since $A = V Λ V^{- 1}$ , the system $U^{'} = A U$ becomes $W^{'} = Λ W$ , which is decoupled: $w_{k}^{'} = λ_{k} w_{k}$ for each $k$ . A one-step method depends on $A$ only through $h A$ , and on each scalar test equation acts as $w_{k}^{n + 1} = R (h λ_{k}) w_{k}^{n}$ , so $W^{n + 1} = R (h Λ) W^{n}$ with $R (h Λ) = diag (R (h λ_{k}))$ . Thus $W^{n} = R (h Λ)^{n} W^{0}$ and, undoing the change of basis, $U^{n} = V R (h Λ)^{n} V^{- 1} U^{0}$ , i.e. $U^{n} = V diag (R (h λ_{k})^{n}) V^{- 1} U^{0}$ . The $k$ -th component stays bounded as $n \to \infty$ iff $∣ R (h λ_{k}) ∣ \leq 1$ , that is $h λ_{k} \in S$ ; the whole solution is bounded iff this holds for every $k$ (a finite intersection of conditions, $V$ being a fixed invertible change of basis). Hence the largest stable step is $h_{m a x} = sup {h : h λ_{k} \in S \forall k}$ , governed by the eigenvalue farthest from the origin in the unstable direction. Rubric: full credit for the eigen-decoupling, the $R (h Λ)^{n}$ formula, and the per-eigenvalue boundedness equivalence.

Advanced results Master

The boundary locus gives the picture; the structure behind it is the rational-approximation theory of $e^{z}$ , the half-plane and exterior-of-disc geometries that distinguish the stable families, and the sharp order-versus-stability obstructions that the next unit turns into the Dahlquist second barrier.

Theorem 1 (stability function as a rational approximation to $e^{z}$ ). For a one-step method of order $p$ , the stability function satisfies $R (z) = e^{z} + C z^{p + 1} + O (z^{p + 2})$ with $C \neq = 0$ , so $R$ is a rational (for Runge-Kutta, $R (z) = 1 + z b^{⊤} (I - z A)^{- 1} 1$ ) approximation to the exponential of order $p$ . The exact propagator over one step of the test equation is $e^{z}$ ; the method replaces it by $R (z)$ , and the entire stability and order theory of one-step methods is the theory of how well a rational function of prescribed degree can both approximate $e^{z}$ near $z = 0$ and keep $∣ R ∣ \leq 1$ on the left half-plane. The Padé table of $e^{z}$ supplies the optimal such approximations: the $(j, k)$ Padé approximant has order $j + k$ , and the diagonal and first two subdiagonals are A-stable ^{[Hairer-Nørsett-Wanner Vol. II §IV.3]}.

Theorem 2 (A-stability is a half-plane containment, characterised by the boundary on the imaginary axis). A method is A-stable when ${z : Re z \leq 0} \subseteq S$ . For a one-step method this holds iff $∣ R (i y) ∣ \leq 1$ for all real $y$ together with $R$ having no pole in the open left half-plane — by the maximum-modulus principle 06.01.01, if $R$ is analytic on ${Re z < 0}$ and bounded by one on its boundary ${Re z = 0} \cup {\infty}$ , then $∣ R ∣ \leq 1$ throughout. The trapezoidal rule realises the extreme case $S = {Re z \leq 0}$ exactly, with $∣ R (i y) ∣ = 1$ for all $y$ ; backward Euler realises the strictly larger exterior-of-disc region with $∣ R (i y) ∣ < 1$ for $y \neq = 0$ . The further refinement $R (\infty) = 0$ , L-stability, holds for backward Euler but not the trapezoidal rule, whose $R (\infty) = - 1$ damps the stiffest transients only marginally ^{[Hairer-Nørsett-Wanner Vol. II §IV.3]}.

Theorem 3 (the multistep region need not be connected or convex, and shrinks with order in the explicit families). For the Adams-Bashforth family the region of absolute stability shrinks rapidly as the order rises: the AB1 (forward Euler) real interval is $(- 2, 0)$ , AB2 is $(- 1, 0)$ , AB3 is $(- 6/11, 0)$ , and higher-order explicit Adams methods have still smaller regions, with portions of the boundary locus producing extra lobes. The implicit Adams-Moulton regions are larger but, beyond AM2 (the trapezoidal rule), are bounded and not A-stable. This monotone shrinking of explicit regions against the unbounded growth of $∣ λ ∣$ in stiff problems is the structural reason explicit multistep methods are unusable on stiff systems, independently of their order ^{[Süli-Mayers §12.12]}.

Theorem 4 (BDF methods and $A (α)$ -stability). The backward-differentiation formulas are the standard stiff multistep family precisely because their regions of absolute stability contain a large sector of the left half-plane. BDF1 (backward Euler) and BDF2 are A-stable; BDF3 through BDF6 are not A-stable but are $A (α)$ -stable, meaning $S$ contains the sector ${z : ∣ ar g (- z) ∣ \leq α}$ with $α$ decreasing from about $8 6^{\circ}$ (BDF3) to about $1 8^{\circ}$ (BDF6); BDF7 and beyond are not even zero-stable 43.10.03, so the family stops at six steps. The boundary locus $z (θ) = ρ (e^{i θ}) / σ (e^{i θ})$ for BDF methods opens leftward and, for the stable members, leaves the entire far left half-plane inside the region ^{[Hairer-Nørsett-Wanner Vol. II §IV.2]}.

Theorem 5 (relative stability and the order-star bridge to the barriers). The map $z \mapsto R (z)$ versus $z \mapsto e^{z}$ is compared by the order star ${z : ∣ R (z) ∣ > ∣ e^{z} ∣}$ , whose finger geometry near $z = 0$ encodes both the order (the number of fingers meeting at the origin) and the A-stability (whether any finger reaches into the left half-plane). This single picture proves the Ehle conjecture — that the $(j, k)$ Padé approximant to $e^{z}$ is A-stable iff $k \leq j \leq k + 2$ — and it is the technique that proves the Dahlquist second barrier for A-stable multistep methods, the subject of 43.10.05 ^{[Hairer-Nørsett-Wanner Vol. II §IV.3]}.

Synthesis. The region of absolute stability is the foundational reason the usable step on a stiff or oscillatory problem is set by geometry rather than accuracy: each eigenvalue $λ_{k}$ of the linearised system contributes one scaled point $h λ_{k}$ , and the step is admissible exactly when every such point lies in $S$ , so the spectral spread of the problem and the shape of the region together fix the step.

This is exactly the fixed- $h$ counterpart of the $h \to 0$ zero-stability of 43.10.03: there the root condition was asked of $ρ$ alone at $z = 0$ , here it is asked of $π (ζ; z) = ρ (ζ) - z σ (ζ)$ at $z = hλ \neq = 0$ , and zero-stability is precisely the $z = 0$ boundary point of $S$ . The boundary-locus map $z = ρ / σ$ is dual to the eigenvalue decoupling of 01.01.08, in which the companion-matrix spectrum that governed zero-stability is replaced by the spectrum of the system matrix feeding the test equation. Putting these together, the stability function $R (z)$ is the central insight that unifies the one-step and multistep pictures — a rational approximation to $e^{z}$ whose order is the method's order and whose left-half-plane sublevel set is its A-stability — and the order-star geometry of ${∣ R ∣ > ∣ e^{z} ∣}$ generalises the boundary locus into the tool that yields the optimal-order results. The bridge to the next unit is the recognition that demanding the region contain the entire left half-plane is A-stability, and combining that demand with the order theory of 43.10.02 forces the Dahlquist second barrier of 43.10.05; the same region machinery appears again in 43.11.02, where the eigenvalues of the discrete Laplacian are dropped into $S$ to read off the parabolic step restriction.

Full proof set Master

Proposition 1 (forward Euler region is the unit disc about $- 1$ ). The region of absolute stability of forward Euler is the closed disc ${z : ∣ z + 1∣ \leq 1}$ , and its real interval is $[- 2, 0]$ .

Proof. Forward Euler on $y^{'} = λ y$ gives $u^{n + 1} = u^{n} + hλ u^{n} = (1 + z) u^{n}$ with $z = hλ$ , so $R (z) = 1 + z$ and $S = {z : ∣1 + z ∣ \leq 1}$ . The condition $∣1 + z ∣ \leq 1$ is precisely $∣ z - (- 1) ∣ \leq 1$ , the closed disc of radius $1$ centred at $- 1$ . Intersecting with the real axis, $z = x$ real gives $∣1 + x ∣ \leq 1 ⟺ - 2 \leq x \leq 0$ , the interval $[- 2, 0]$ . $□$

Proposition 2 (trapezoidal region is the closed left half-plane). The trapezoidal rule has $R (z) = (1 + z /2) / (1 - z /2)$ and region of absolute stability $S = {z : Re z \leq 0}$ ; moreover $∣ R (i y) ∣ = 1$ for every real $y$ , so the boundary of $S$ is the imaginary axis.

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)$ (the value $z = 2$ excluded as a pole). 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. $$ Hence $∣ R (z) ∣^{2} = ∣1 + z /2 ∣^{2} /∣1 - z /2 ∣^{2} \leq 1 ⟺ 2 x \leq 0 ⟺ Re z \leq 0$ , so $S = {Re z \leq 0}$ , and the difference vanishing at $x = 0$ gives $∣ R (i y) ∣ = 1$ for all real $y$ , so the imaginary axis is exactly the boundary. $□$

Proposition 3 (boundary locus passes through the origin tangent to the imaginary axis for any consistent method). For a consistent linear multistep method, the boundary locus $z (θ) = ρ (e^{i θ}) / σ (e^{i θ})$ satisfies $z (0) = 0$ and, to first order in $θ$ , $z (θ) = i θ + O (θ^{2})$ , so the locus leaves the origin tangent to the imaginary axis.

Proof. Consistency gives $ρ (1) = 0$ and $ρ^{'} (1) = σ (1)$ 43.10.02, with $σ (1) \neq = 0$ (else the method would be degenerate). At $θ = 0$ , $ζ = e^{i θ} = 1$ , so $z (0) = ρ (1) / σ (1) = 0/ σ (1) = 0$ . Differentiate $z (θ) = ρ (e^{i θ}) / σ (e^{i θ})$ in $θ$ . Writing $ζ = e^{i θ}$ so $d ζ / d θ = i ζ$ , the quotient rule gives $$ z'(\theta) = \frac{\rho'(\zeta),i\zeta,\sigma(\zeta) - \rho(\zeta),\sigma'(\zeta),i\zeta}{\sigma(\zeta)^2}. $$ At $θ = 0$ , $ζ = 1$ , $ρ (1) = 0$ , so the second term drops and $$ z'(0) = \frac{\rho'(1),i,\sigma(1)}{\sigma(1)^2} = i,\frac{\rho'(1)}{\sigma(1)} = i,\frac{\sigma(1)}{\sigma(1)} = i, $$ using $ρ^{'} (1) = σ (1)$ . Hence $z (θ) = z (0) + z^{'} (0) θ + O (θ^{2}) = i θ + O (θ^{2})$ , a curve through the origin with tangent direction $i$ — the imaginary axis. This is the structural reason consistent methods all have the imaginary axis tangent to their stability boundary at the origin: near $z = 0$ the test problem is a pure oscillation, and a consistent method must keep marginal oscillations marginal. $□$

Proposition 4 (eigen-decoupling of the step restriction). Let $U^{'} = A U$ with $A \in C^{d \times d}$ diagonalisable, $A = V Λ V^{- 1}$ , $Λ = diag (λ_{1}, \dots, λ_{d})$ . A one-step method with stability function $R$ produces $U^{n} = V R (h Λ)^{n} V^{- 1} U^{0}$ , and $sup_{n} ∥ U^{n} ∥ < \infty$ for all initial data if and only if $h λ_{k} \in S$ for every $k$ .

Proof. Put $W^{n} = V^{- 1} U^{n}$ . A one-step method applied to $U^{'} = A U$ advances by a matrix function $U^{n + 1} = R (h A) U^{n}$ , where for a Runge-Kutta method $R (h A) = I + h A b^{⊤} (I - h AA)^{- 1} 1$ and in general $R (h A)$ is the same rational/polynomial expression in $h A$ that defines $R$ on scalars; this uses only that the method is linear and depends on $A$ through $h A$ . Since $A = V Λ V^{- 1}$ and $R$ is a rational function, $R (h A) = V R (h Λ) V^{- 1}$ with $R (h Λ) = diag (R (h λ_{k}))$ (any rational function commutes with the similarity, applied entrywise on the diagonal). Therefore $U^{n} = R (h A)^{n} U^{0} = V R (h Λ)^{n} V^{- 1} U^{0}$ and $W^{n} = R (h Λ)^{n} W^{0}$ , decoupled into scalar recurrences $w_{k}^{n} = R (h λ_{k})^{n} w_{k}^{0}$ . Each is bounded in $n$ iff $∣ R (h λ_{k}) ∣ \leq 1$ , i.e. $h λ_{k} \in S$ . As $V$ is a fixed invertible matrix, $∥ U^{n} ∥ \leq ∥ V ∥ ∥ V^{- 1} ∥ max_{k} ∣ R (h λ_{k}) ∣^{n} ∥ W^{0} ∥$ is bounded over $n$ and all data iff each factor is, i.e. iff $h λ_{k} \in S$ for every $k$ . $□$

Connections Master

The zero-stability and Dahlquist-equivalence theory of 43.10.03 is the $z = 0$ edge of this unit: zero-stability asks the root condition of $ρ$ alone as $h \to 0$ , while absolute stability asks the same of $π (ζ; z) = ρ (ζ) - z σ (ζ)$ at a fixed $z = hλ \neq = 0$ . The companion-matrix root-tracking that proved the equivalence theorem there is exactly the continuity-of-roots argument that gives the boundary locus here, and the strong-versus-weak-stability classification of spurious roots from that unit is the input that decides which lobe of the boundary locus is the stable region.
The linear multistep order theory of 43.10.02 supplies the polynomials $ρ, σ$ , the consistency conditions $ρ (1) = 0$ , $ρ^{'} (1) = σ (1)$ , and the construction of the Adams and BDF families whose regions are computed here. Consistency is precisely what makes the boundary locus pass through the origin tangent to the imaginary axis (Proposition 3), and the order of a method controls how closely its stability function $R (z)$ approximates $e^{z}$ , linking the order constants of that unit to the stability geometry of this one.
The eigenvalue and diagonalisation theory of 01.01.08 is the engine of the eigen-decoupling (Proposition 4): a constant-coefficient linear system is reduced to independent scalar test equations by the eigenbasis of its matrix, so the abstract region $S$ acquires its operational meaning — the step is admissible exactly when every $h λ_{k}$ lands in the region. The spectral radius and Jordan structure that governed power-boundedness for zero-stability reappear as the per-eigenvalue boundedness condition here.
The holomorphic-function and complex-plane theory of 06.01.01 is the setting in which the region $S \subseteq C$ , the rational stability function $R (z)$ , and the maximum-modulus characterisation of A-stability all live: $R$ is a rational function of the complex variable $z = hλ$ , its poles locate the implicit-method singularities, and the proof that a half-plane-bounded analytic $R$ with $∣ R ∣ \leq 1$ on the imaginary axis is A-stable is a maximum-modulus argument.
The parabolic finite-difference theory of 43.11.02 applies this machinery to PDEs: the method-of-lines semidiscretisation of the heat equation yields a stiff linear system whose matrix is the discrete Laplacian, and the parabolic step restriction $Δ t \leq \frac{1}{2} Δ x^{2}$ for explicit time-stepping is read off by dropping the eigenvalues of that discrete Laplacian into forward Euler's region $S$ — the eigen-decoupling of Proposition 4 carried to the PDE setting.

Historical & philosophical context Master

The absolute-stability programme is due to Germund Dahlquist, whose 1963 paper in BIT ^{[Dahlquist 1963]} posed the "special stability problem" of integrating stiff equations and introduced A-stability as the requirement that the region of absolute stability contain the entire left half-plane; the same paper proved the second Dahlquist barrier, that no A-stable linear multistep method exceeds order two. The test equation $y^{'} = λ y$ and the practice of reading stability off a single complex parameter $z = hλ$ were already standard in the difference-equation analyses of the 1950s, but Dahlquist's framing turned the question into the geometry of a region in the complex plane. The stability function $R (z)$ of Runge-Kutta methods as a rational approximation to $e^{z}$ was developed by Byron Ehle in his 1969 thesis, where he conjectured the now-proven characterisation of which Padé approximants of $e^{z}$ are A-stable.

The boundary-locus method — tracing $z (θ) = ρ (e^{i θ}) / σ (e^{i θ})$ to plot a region — was systematised in the multistep-method literature of the 1960s and appears in Peter Henrici's 1962 monograph and Jack Lambert's subsequent texts. The sharpest modern proofs of the order-versus-stability barriers came with the order-star technique introduced by 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 the Ehle conjecture together with the Daniel-Moore conjecture on the maximal order of A-stable methods.

Bibliography Master

@article{dahlquist1963stiff,
  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}
}

@phdthesis{ehle1969,
  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{wannerhairernorsett1978stab,
  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{hairerwanner1996,
  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{hnw1993stab,
  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}
}

@book{lambert1991,
  author    = {Lambert, John D.},
  title     = {Numerical Methods for Ordinary Differential Systems: The Initial Value Problem},
  publisher = {John Wiley \& Sons},
  year      = {1991}
}

@book{leveque2007fdmstab,
  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{sulimayers2003stab,
  author    = {S\"{u}li, Endre and Mayers, David F.},
  title     = {An Introduction to Numerical Analysis},
  publisher = {Cambridge University Press},
  year      = {2003}
}

Prerequisites

43.10.03
43.10.02
01.01.08
06.01.01

Tier anchors

beginner: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §7.1-7.3 (the test equation $y'=\lambda y$, the amplification factor of a step, and the stability region drawn in the $z=h\lambda$ plane); Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §12.11 (absolute stability of Euler and the trapezoidal rule, the interval of absolute stability)
intermediate: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §7.1-7.4 (the stability function $R(z)$ of one-step methods, the characteristic equation $\pi(\zeta;z)=\rho(\zeta)-z\sigma(\zeta)$ for multistep methods, regions of absolute stability and the boundary-locus method); Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §12.11-12.12 (absolute stability regions for Euler, trapezoidal, Adams, and BDF methods)
master: Hairer-Nørsett-Wanner 1993 *Solving Ordinary Differential Equations I: Nonstiff Problems* 2e (Springer) §III.3 and Hairer-Wanner 1996 *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems* 2e (Springer) §IV.2-IV.3 (stability functions, stability regions, the boundary locus, order stars and the Ehle conjecture); Dahlquist 1963 *A special stability problem for linear multistep methods* (BIT 3) (the test-equation stability framework and the A-stability programme)

References

LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems · SIAM 2007. Chapter 7 develops absolute stability from the Dahlquist test equation $y'=\lambda y$. §7.1 applies a one-step method to the test problem to obtain $u^{n+1}=R(z)u^n$ with $z=h\lambda$ and the *stability function* $R(z)$ ($R(z)=1+z$ for forward Euler, $1/(1-z)$ for backward Euler, $(1+z/2)/(1-z/2)$ for the trapezoidal rule, and the degree-4 Taylor polynomial of $e^z$ for classical RK4), and defines the *region of absolute stability* $\{z\in\mathbb C:|R(z)|\le1\}$. §7.3 extends this to linear multistep methods through the *characteristic equation* $\pi(\zeta;z)=\rho(\zeta)-z\sigma(\zeta)=0$: the region is $\{z:\text{all roots }\zeta_i(z)\text{ of }\pi(\cdot;z)\text{ satisfy }|\zeta_i|\le1,\text{ boundary roots simple}\}$, and the *boundary-locus method* traces $\partial$ of the region as $z(\theta)=\rho(e^{i\theta})/\sigma(e^{i\theta})$, $\theta\in[0,2\pi)$. §7.4 treats coupled linear systems by diagonalising $U'=AU$ so that each eigenvalue $\lambda_k$ of $A$ must satisfy $h\lambda_k\in$ the region, and uses the region to bound the usable step on a decaying or oscillatory problem.
Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis · Cambridge University Press 2003. §12.11 defines absolute stability of a numerical method through the model problem $y'=\lambda y$ with $\mathrm{Re}\,\lambda<0$, requiring $|u^n|\to0$ as $n\to\infty$, and computes the *interval of absolute stability* on the negative real axis for Euler's method ($-2<h\lambda<0$) and the unconditional stability of the implicit Euler and trapezoidal rules; §12.12 computes the absolute stability regions of the Adams-Bashforth, Adams-Moulton, and backward-differentiation families through the boundary locus $h\lambda=\rho(e^{i\theta})/\sigma(e^{i\theta})$ and contrasts the shrinking explicit regions with the large implicit ones, motivating implicit methods for stiff problems.
Hairer, E., Nørsett, S. P. & Wanner, G. — Solving Ordinary Differential Equations I and II · Springer, 2nd revised editions 1993/1996. Vol. I §III.3 introduces the stability function $R(z)$ of one-step and Runge-Kutta methods as $R(z)=1+z b^\top(I-zA)^{-1}\mathbf 1$ for a Butcher tableau $(A,b)$, identifies it as a rational approximation to $e^z$ of order equal to the method order, and draws the stability domains. Vol. II §IV.2 defines the region of absolute stability, the boundary locus, $A$-stability ($\{\mathrm{Re}\,z\le0\}\subseteq$ region), $A(\alpha)$-stability and $L$-stability ($R(\infty)=0$); §IV.3 develops the order-star theory of Wanner-Hairer-Nørsett that proves the optimal-order results (the Ehle conjecture: the $(j,k)$ Padé approximant to $e^z$ is $A$-stable iff $k\le j\le k+2$) and connects the geometry of $\{|R(z)|=|e^z|\}$ to attainable order under $A$-stability.

Estimated time

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