Initial-Value Integration and Stiffness

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

An adaptive ODE solver controls local error, but local accuracy alone does not decide whether an integration is efficient or stable. In a stiff system, rapidly decaying modes impose a tiny stability step on an explicit method even when the observable of interest evolves slowly. Choosing an integrator therefore requires a joint diagnosis of timescales, stability region, invariants, nonlinear-solve cost, and the finite interval on which accuracy is claimed.

Required background. Ordinary Differential Equations supplies initial-value problems, local error, and stability functions.

Helpful background. Formulations of General Relativity provides the relativistic systems that lead to multiscale reductions, while Partial Differential Equations and Green Functions helps interpret ODE systems obtained after spatial discretization or separation.

Stiffness is method dependent

For $y'=Ay$ , an explicit one-step method advances $y_{n+1}=R(hA)y_n$ . Each relevant $h\lambda_j(A)$ must lie in the method’s stability region; for nonnormal matrices, eigenvalues alone do not control transient amplification. If one eigenvalue has a large negative real part, an explicit method may require $h=O(1/\lvert\lambda_j\rvert)$ long after that mode has decayed below the requested error. An implicit method can remove that particular restriction, at the price of solving linear or nonlinear systems.

The scalar decay equation $z'=-\lambda z$ , with $\lambda\gt0$ , makes the distinction concrete. Forward Euler has amplification $1-h\lambda$ and needs $0\le h\lambda\le2$ for bounded steps. Backward Euler has amplification $(1+h\lambda)^{-1}$ , which tends to zero as $h\lambda\to\infty$ . The trapezoidal rule is also stable for every positive $h\lambda$ , but its amplification $(1-h\lambda/2)/(1+h\lambda/2)$ tends to $-1$ : unresolved fast modes alternate rather than disappear. A-stability includes the whole left half-plane in the stability region; L-stability additionally damps infinitely stiff modes to zero. This is why two implicit methods can behave differently on the same decaying transient.

This is not an intrinsic yes-or-no label attached to an equation. Its validity regime depends on scaling, interval length, tolerance, initial data, and observable. Hairer and Wanner develop the stability and implementation of stiff solvers in Solving Ordinary Differential Equations II, Chapter IV §§ 1–3, pp. 2–55. In a relativistic code, a poor variable scaling or nearly singular coordinate choice can create computational stiffness without introducing a new physical timescale.

Exact manufactured test

Before using a solver on a system without a closed form, test the scalar problem

y'=-\lambda(y-\cos t)-\sin t, \qquad y(0)=1,

whose exact solution is $y(t)=\cos t$ . With $\lambda=10^4$ on $0\le t\le10$ , the fast transient is absent from the exact initial data but still appears in the Jacobian. Compare an explicit Runge–Kutta method, an implicit Radau method, and a backward differentiation formula using the tolerance pairs

\begin{aligned} (\mathrm{rtol},\mathrm{atol})&=(10^{-6},10^{-9}),\\ &\quad(10^{-9},10^{-12}),\\ &\quad(10^{-12},10^{-15}). \end{aligned}

Absolute tolerance matters because the exact solution crosses zero. Record accepted and rejected steps, right-hand-side and Jacobian evaluations, linear solves, wall time, and

E_\infty=\max_n\lvert y_n-\cos t_n\rvert.

This benchmark separates accuracy from work. A method passes when error follows the tolerance until roundoff or implementation order limits it, and when tighter tolerances yield convergence rather than merely more steps. To check the boundary residual machinery later in the chapter, use the distinct second-order problem $z''=-\cos t$ with $z(0)=1$ and $z(10)=\cos10$ , whose exact solution is $z=\cos t$ . Imposing two endpoint values on the first-order equation above would overdetermine it.

Robertson kinetics as a stiffness stress test

The Robertson system is

\begin{aligned} y_1'&=-0.04y_1+10^4y_2y_3,\\ y_2'&=0.04y_1-10^4y_2y_3-3\times10^7y_2^2,\\ y_3'&=3\times10^7y_2^2, \end{aligned}

with $(y_1,y_2,y_3)(0)=(1,0,0)$ . It spans widely separated reaction times. There is no simple exact trajectory, but two discriminating checks are available:

y_1+y_2+y_3=1, \qquad y_i\ge0.

Integrate over $0\le t\le10^{11}$ at the same three relative tolerances, using componentwise

\boldsymbol{\mathrm{atol}} =\mathrm{rtol}\,(10^{-3},10^{-7},10^{-3}).

Sample 401 output points uniform in $\log(1+t)$ . Compare Radau and BDF with an independent 100-bit reference using $\mathrm{rtol}=10^{-15}$ and $\boldsymbol{\mathrm{atol}}=(10^{-18},10^{-22},10^{-18})$ . Define the reported trajectory discrepancy by

E_w=\max_{n,i} \frac{\lvert y_i(t_n)-y_i^{\rm ref}(t_n)\rvert} {\mathrm{atol}_i+\mathrm{rtol}\max(\lvert y_i(t_n)\rvert,\lvert y_i^{\rm ref}(t_n)\rvert)}.

The accepted comparison reports $E_w$ , the maximum conservation residual, minimum component, and work. Interpolate the reference with its solver’s dense output rather than linearly across the fast initial layer. Solver architectures that expose tolerances, Jacobians, and linear-solver choices are described in Hindmarsh et al. 2005, §§ 2–4.

Jacobians, structure, and failure modes

Implicit methods are only as reliable as their nonlinear solves. Supply either an analytic Jacobian or a finite-difference Jacobian with a documented perturbation scale, and compare both. A wrong Jacobian may lead to step rejection, false convergence of Newton iteration, or excessive work without obvious solution error. Sparse or band structure should be declared rather than inferred from one parameter point.

Structure-preserving methods can matter for Hamiltonian orbital systems: a symplectic scheme may control long-time energy drift better than a generic high-order method. That advantage does not make it suitable for dissipative stiffness. Match the method to the invariant and timescale actually being tested.

The environment record pins solver library and version, arithmetic precision, compiler or interpreter, method, tolerances, Jacobian mode, linear solver, event handling, maximum step, and output grid. A clean rerun executes the manufactured and Robertson inputs, emits machine-readable trajectories and work counters, and repeats tolerance and precision sweeps. Reject a result if conservation improves only because the solver projects the state silently, positivity fails, the reference uses the same defective Jacobian, or requested tolerances are below an observed saturation floor.

Check your understanding

An explicit solver and a Radau solver agree to $10^{-8}$ on the Robertson variables, but the explicit run uses a million steps while Radau uses a few hundred. Does the step-count difference by itself prove that Radau is more accurate?

Solution

No. It demonstrates a stability-cost advantage for the implicit method on this input, not greater accuracy. Compare both with an independent tighter reference, inspect conservation and positivity, and report error against work at several tolerances. The explicit solution may be equally accurate but inefficient because stability, rather than error, limits its step.

References

Hairer, E., and Wanner, G. (1996), Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd revised edition, Springer. DOI.
Hindmarsh, A. C., et al. (2005), “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers,” ACM Transactions on Mathematical Software 31, 363–396. DOI.