Unit 7: Absolute Stability Analysis, Stiff Differential Equations & Linear Multi-Step Methods

1. Germund Dahlquist's Test Problem

To analyze whether numerical solutions remain bounded as $t \to \infty$ for fixed step size $h > 0$, Germund Dahlquist (1956) introduced the scalar complex test equation:

The exact solution is $y(t) = y_0 e^{\lambda t}$. Since $\text{Re}(\lambda) < 0$, the exact solution asymptotically decays to zero:

Applying a numerical method with step size $h$ yields a discrete recurrence relation:

where $R(z): \mathbb{C} \to \mathbb{C}$ is the stability function (or amplification factor). After $n$ steps:

Definition: Region of Absolute Stability

The region of absolute stability $\mathcal{S} \subset \mathbb{C}$ is the set of complex parameters $z = h\lambda$ for which the numerical solution decays or remains bounded:

Its boundary is the locus $\partial \mathcal{S} = \{ z \in \mathbb{C} : |R(z)| = 1 \}$.

2. Stability Functions for Classic Methods

1. Forward Euler:

Stability condition: $|1 + z| < 1$, which is an open disk of radius 1 centered at $z = -1$. On the negative real axis ($\lambda \in \mathbb{R}^-$), stable step sizes require:

2. Backward Euler (Implicit):

Stability condition: $|1 - z| > 1$, which is the exterior of the disk of radius 1 centered at $z = 1$. Notice that $\mathcal{S}$ includes the entire left half-plane $\mathbb{C}^- = \{z : \text{Re}(z) < 0\}$.

3. Implicit Trapezoidal Rule (Crank-Nicolson):

For $\text{Re}(z) < 0$, $|1 + z/2| < |1 - z/2|$, so $|R(z)| < 1$ holds for all $z \in \mathbb{C}^-$.

4. Explicit Runge-Kutta Methods of Order $p \le 4$:

For any explicit $p$-stage, $p$-th order RK method ($p = 1, 2, 3, 4$), the stability function is the truncated Taylor polynomial of $e^z$:

For classical RK4:

The interval of absolute stability on the real axis is approximately $(-2.78, 0)$, and it contains a segment of the imaginary axis $[-2\sqrt{2}i, 2\sqrt{2}i]$, making RK4 capable of integrating undamped wave oscillators.

1. The Phenomenon of Stiffness

Consider a linear system of differential equations:

with eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_d$ such that $\text{Re}(\lambda_i) < 0$.

Definition: Stiffness Ratio

The system is called stiff if:

1. $\text{Re}(\lambda_i) < 0$ for all $i = 1, \dots, d$.
2. The stiffness ratio:

Physical Reality: Stiff systems describe multi-scale dynamics—such as chemical kinetics, combustion, or structural mechanics—where fast transient modes decay almost instantaneously while slow master modes dominate the physical trajectory.

The Failure of Explicit Methods:

To keep explicit methods from experiencing exponential divergence, the step size $h$ is severely constrained by the fastest decaying eigenvalue:

Even after the fast transient has decayed to zero ($10^{-10}$), an explicit integrator is forced to take millions of tiny steps purely to satisfy numerical stability, even though accuracy would permit a step size $10^6$ times larger!

2. A-Stability and L-Stability

To solve stiff problems efficiently, we require methods whose stability region encompasses the entire physics of decaying systems.

Definition 1: A-Stability (Dahlquist, 1963)

A numerical method is A-stable if its stability region contains the entire open left half-plane:

For an A-stable method, any step size $h > 0$ produces bounded solutions, eliminating the stability bottleneck completely!

Definition 2: L-Stability (Stiff Decay)

While the Trapezoidal rule is A-stable, its stability function satisfies:

As a consequence, very stiff modes do not decay to zero; they oscillate wildly between $\pm y_0$!

A method is L-stable if:

1. It is A-stable.
2. Its stability function vanishes at infinity:

Backward Euler is L-stable since $\lim_{z \to -\infty} \frac{1}{1 - z} = 0$. Stiff solvers in production (e.g., Radau IIA, BDF) are L-stable, ensuring that stiff transients are damped out instantly.

1. General Linear Multi-Step Framework

A linear $k$-step method uses numerical values and function derivatives from the $k$ preceding points $(y_{n+j}, f_{n+j})$ to compute $y_{n+k}$:

where $\alpha_k \ne 0$ and $|\alpha_0| + |\beta_0| > 0$. By standard convention, normalize $\alpha_k = 1$.

• If $\beta_k = 0$, the method is explicit ($y_{n+k}$ appears only on the left).
• If $\beta_k \ne 0$, the method is implicit ($y_{n+k}$ appears on both sides inside $f(t_{n+k}, y_{n+k})$).

2. The Adams Family of Integrators

Integrating $y' = f(t, y)$ from $t_{n+k-1}$ to $t_{n+k}$:

Replace $f(t, y(t))$ by an interpolating polynomial $P(t)$:

1. Adams-Bashforth Methods (Explicit):

$P(t)$ interpolates past values $\{f_{n}, f_{n+1}, \dots, f_{n+k-1}\}$ at $k$ points:

• AB1 (Euler): $y_{n+1} = y_n + h f_n$
• AB2: $y_{n+2} = y_{n+1} + \frac{h}{2} (3f_{n+1} - f_n)$
• AB3: $y_{n+3} = y_{n+2} + \frac{h}{12} (23f_{n+2} - 16f_{n+1} + 5f_n)$
• AB4: $y_{n+4} = y_{n+3} + \frac{h}{24} (55f_{n+3} - 59f_{n+2} + 37f_{n+1} - 9f_n)$

2. Adams-Moulton Methods (Implicit):

$P(t)$ interpolates $\{f_{n}, \dots, f_{n+k-1}, f_{n+k}\}$ at $k+1$ points:

• AM1 (Backward Euler): $y_{n+1} = y_n + h f_{n+1}$
• AM2 (Trapezoidal): $y_{n+1} = y_n + \frac{h}{2} (f_{n+1} + f_n)$
• AM3: $y_{n+2} = y_{n+1} + \frac{h}{12} (5f_{n+2} + 8f_{n+1} - f_n)$
• AM4: $y_{n+3} = y_{n+2} + \frac{h}{24} (9f_{n+3} + 19f_{n+2} - 5f_{n+1} + f_n)$

3. Predictor-Corrector Architecture (PECE)

To avoid solving nonlinear systems for implicit Adams-Moulton methods, modern solvers use a Predict-Evaluate-Correct-Evaluate (PECE) sequence:

1. P (Predict): Compute initial guess $y_{n+1}^{(0)}$ using explicit Adams-Bashforth:

2. E (Evaluate): Evaluate $f_{n+1}^{(0)} = f(t_{n+1}, y_{n+1}^{(0)})$.

3. C (Correct): Refine using implicit Adams-Moulton:

4. E (Evaluate): Final update $f_{n+1} = f(t_{n+1}, y_{n+1})$ for the next time step.

This achieves order 4 accuracy with only two function evaluations per step, compared to RK4's four evaluations!

1. Characteristic Polynomials

Associated with the linear multi-step method $\sum_{j=0}^k \alpha_j y_{n+j} = h \sum_{j=0}^k \beta_j f_{n+j}$ are the two characteristic polynomials:

Consistency Conditions:

A linear multi-step method is consistent if and only if:

2. Zero-Stability and the Dahlquist Root Condition

Applying the method to the trivial differential equation $y' = 0$ yields the difference equation:

whose general solution is a linear combination of $r_m^n$, where $r_m$ are roots of $\rho(r) = 0$.

Definition (Dahlquist Root Condition / Zero-Stability): A linear multi-step method is zero-stable if all roots $r$ of the first characteristic polynomial $\rho(r) = 0$ satisfy:

1. $|r| \le 1$ (all roots lie within or on the closed unit disk).
2. Any root lying on the unit circle ($|r| = 1$) is simple (multiplicity 1).

Theorem (Dahlquist Equivalence Theorem): For any linear multi-step method:

3. The Dahlquist Stability Barriers

Germund Dahlquist established two profound mathematical bounds that dictate what numerical integrators can and cannot do:

Dahlquist's First Barrier (1956):

The order of convergence $p$ of a zero-stable $k$-step linear multi-step method cannot exceed:

Dahlquist's Second Barrier (1963):

1. No explicit linear multi-step method can be A-stable.
2. An A-stable linear multi-step method cannot have order of accuracy exceeding $p = 2$.
3. Among all second-order A-stable linear multi-step methods, the Implicit Trapezoidal Rule has the smallest asymptotic error constant:

This theorem proves that high-order A-stable methods cannot be achieved with linear multi-step formulas, motivating the development of Backward Differentiation Formulas (BDF) (which sacrifice A-stability for $A(\alpha)$-stability at orders $k \le 6$) and Implicit Runge-Kutta methods (which overcome the barrier by utilizing multiple internal stages).