Three Methods, Three Regions
For the test equation dy/dx = λy, three explicit ODE methods have the following stability regions in the complex hλ-plane:
Euler's method (first-order): stability region is the disk |1 + hλ| ≤ 1, a circle of radius 1 centered at (-1, 0). Real negative hλ must lie in [-2, 0].
Runge-Kutta 2 (midpoint method) (second-order): stability region is |1 + hλ + (hλ)²/2| ≤ 1. Larger than Euler's disk, but still bounded.
Runge-Kutta 4 (fourth-order): stability region satisfies |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Real negative hλ extends to approximately -2.785. The region is substantially larger than Euler's.
Backward Euler (implicit): stability region is the entire complex plane except the disk |1 - hλ|⁻¹ > 1, equivalently |1/(1-hλ)| ≤ 1. For λ in the left half-plane (Re(λ) < 0), this is unconditionally stable — no constraint on h from stability.
The Amplification Function
For any Runge-Kutta method, the amplification factor per step R(hλ) is a polynomial approximation to e^(hλ):
- Euler: R(z) = 1 + z (truncated at degree 1)
- RK2: R(z) = 1 + z + z²/2 (truncated at degree 2)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (truncated at degree 4)
The stability region is {z ∈ ℂ : |R(z)| ≤ 1}. The true solution amplification: |e^z| = e^(Re(z)). For Re(z) < 0 (stable ODE), the true solution decays. The numerical method is stable if |R(z)| ≤ 1 — matching the decaying behavior.
Pure Imaginary Eigenvalues: Oscillatory Systems
Many physical systems have purely imaginary eigenvalues: λ = iω (oscillations without damping). The spring-mass system, orbital mechanics, pendulum dynamics.
For λ = iω: hλ = ihω lies on the imaginary axis.
Euler's stability on the imaginary axis: |1 + ihω|² = 1 + (hω)² > 1 for any h > 0. Euler is unstable for any step size on purely imaginary eigenvalues. The computed 'oscillation' grows without bound.
RK4's stability on the imaginary axis: the stability region extends to approximately |hω| ≤ 2.83 on the imaginary axis. For small enough h, RK4 handles undamped oscillations. Euler cannot.
This geometry is why Euler fails on conservative systems (spring-mass, orbits, wave equations) even with small h, while RK4 handles them well.
The Geometry of Stiff Problems
A stiff ODE system has eigenvalues with very different magnitudes. The stiffness ratio: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Why stiffness is expensive for explicit solvers:
Stability requires h·max|λᵢ| ≤ C (where C depends on the method). The most negative eigenvalue sets the bound.
Accuracy for the slow dynamics requires h·min|λᵢ| ≥ ε (resolve the slowest mode adequately).
If κ is large, these two requirements force a tiny h: small enough for stability of the fast mode, large enough to sample the slow mode. The number of steps scales with κ.
Geometric picture in the eigenvalue spectrum: the eigenvalues of the Jacobian ∂f/∂y form a set of points in the complex plane. An explicit solver's stability region must contain all points h·λᵢ. If eigenvalues span from -1 to -1000, the stability region must cover a range of 1000 along the real axis — requiring very small h.
Implicit solvers: backward Euler's stability region covers the entire left half-plane. All eigenvalues with Re(λ) < 0 are automatically inside the stability region regardless of h. The constraint on h comes only from accuracy, not stability.
Stiffness Ratio & Cost
Consider a chemical reaction network with fast reactions (time scale 10⁻⁶ s) and slow reactions (time scale 1 s).
Stiffness ratio: κ = 10⁶ / 1 = 10⁶.
With RK4 (stability limit h·|λ_max| ≤ 2.785): h_max = 2.785 / 10⁶ ≈ 2.8 × 10⁻⁶ s.
To integrate over 10 s of reaction time: steps = 10 / (2.8 × 10⁻⁶) ≈ 3.6 × 10⁶.
With backward Euler (unconditionally stable): h can be chosen for accuracy of the slow reactions. h = 10⁻² s (100 samples over 1 s scale). Steps = 10 / 10⁻² = 1000.
Cost ratio: explicit 3.6 million steps vs implicit 1000 steps — a factor of 3600. Each implicit step requires solving a linear system (cost per step is higher), but the total cost is much lower for very stiff problems.
Why n-Dimensional Tubes Are Not What You Think
In 2D, a 'tube' of radius ε around a curve C is the set of points within distance ε of C. The cross-section is a circle of radius ε. The volume of the tube grows proportionally with its length.
In n dimensions, the tube geometry changes fundamentally, due to the phenomenon from Chapter 9:
The n-dimensional corner paradox: in n-dimensional space, almost all the volume of an n-dimensional hypercube lies in the corners — not in the central region. As n increases, the fraction of volume within distance ε of the center goes to zero for any fixed ε.
Applied to ODE solution tubes:
In 2D: if the true solution passes through the center of a tube, most nearby points are close to the curve. Small perturbations keep you near the true solution.
In high dimensions: most points within the tube's bounding box are actually far from the true solution curve. The tube's 'volume' is dominated by the corners — regions that are far from the center in multiple dimensions simultaneously.
Consequence for simulation: with 28 coupled ODEs (Hamming's Navy intercept problem), a perturbation of size ε in each dimension can produce a total displacement of ε√28 ≈ 5.3ε from the true solution. The tube must be understood in terms of the L2 norm across all dimensions, not just the maximum displacement in any one dimension.
Stability in high dimensions: a system where each component decays independently (each eigenvalue has negative real part) may still show large combined displacements because the components' errors add in L2 norm. The 28-dimensional tube is not just 28 independent 1-dimensional tubes — the geometry couples them.
From Geometry to Design
The geometric insights of Chapters 18-20 come together as a set of design principles for numerical simulation:
Step size selection: h must place h·λ inside the stability region for every eigenvalue. For stiff systems, implicit methods remove the stability constraint, leaving only accuracy requirements.
Error accumulation in high dimensions: the global error is a vector in n-dimensional space. Its norm grows as √n times the per-component error. High-dimensional simulations need tighter per-step accuracy requirements.
Feedback as stabilizer: if the simulation incorporates feedback (the computed output influences subsequent inputs, as in a guidance system), convergent feedback damps errors. The simulation can tolerate imprecise inputs for quantities inside the feedback loop.
The instability as signal: for problems with divergent direction fields, instability can be exploited: the direction of divergence carries information about the initial condition error, enabling corrective adjustment.