Stability Regions
A method's stability region is the set of step-times-eigenvalue values for which it stays bounded, explaining why stiff problems need implicit solvers.
The linear test equation
Stability of ODE solvers is analyzed on the scalar test equation y' = lambda y, whose exact solution decays when lambda has negative real part. Applying a method gives y_{n+1} = R(z) y_n with z = h lambda, where R is the method's amplification factor. The method is stable for a given z when |R(z)| <= 1.
The stability region
The set of complex z with |R(z)| <= 1 is the method's stability region. For a real system, h times each eigenvalue of the Jacobian must lie inside this region for the numerical solution to remain bounded. This turns the abstract idea of stability into a concrete constraint on the step size.
Explicit versus implicit shapes
- Explicit methods have bounded stability regions, so large negative eigenvalues force small h.
- Forward Euler's region is a unit disk centered at -1: h|lambda| must stay below 2.
- A-stable methods contain the entire left half-plane, allowing any step on decaying problems.
- Backward Euler and the trapezoidal rule are A-stable; that is why they suit stiff systems.
L-stability and beyond
A-stability alone can leave fast modes lightly damped, causing oscillatory residuals. An L-stable method additionally forces R(z) toward zero as z goes to negative infinity, aggressively killing the stiffest components. Backward Euler is L-stable; the trapezoidal rule is not, which is why the latter can ring on very stiff problems.
# amplification factors
R_forward_euler = lambda z: 1 + z
R_backward_euler = lambda z: 1/(1 - z)
R_trapezoid = lambda z: (1 + z/2)/(1 - z/2)Stability-region analysis guides solver and step choices for the stiff, multi-time-scale equations arising in breeder Hyperion plasma simulations.