Stiffness: When Explicit Steps Collapse

A chemical kinetics network written as a system of rate equations behaves in two time scales at once. Some species react in microseconds and settle into equilibrium almost immediately; others evolve over seconds. The fast species are not the interesting part of the answer, but they determine the step size of every explicit integrator. This is stiffness, and it has a precise operational meaning: it is not about the equations being difficult to solve, but about the ratio between the fastest transient the dynamics can produce and the time scale over which the solution changes.

The operational definition

Take a linear constant-coefficient system $y’ = Ay$ with eigenvalues $\lambda_i$ of $A$. Each mode evolves as $e^{\lambda_i t}$. A mode with $\operatorname{Re}\lambda_i < 0$ decays; its rate is $|\operatorname{Re}\lambda_i|$. The stiffness ratio is

$$ S = \frac{\max_i |\operatorname{Re}\lambda_i|}{\min_i |\operatorname{Re}\lambda_i|} $$

over the modes whose amplitude is non-negligible, with the sign convention that decaying modes contribute. A problem is stiff when some decaying mode is fast and the solution of interest varies on a much longer scale. There is no sharp threshold. The practical cutoff depends on the method and the tolerance, and different authors put the transition anywhere from a ratio of a few tens to several thousand. The number is a property of the problem and the accuracy target together.

For nonlinear $y’ = f(t,y)$, the role of $A$ is taken by the Jacobian $J = \partial f/\partial y$, and the definition is applied locally, at points along the trajectory rather than once for the whole problem. A system can be stiff in one region of state space and non-stiff in another. Combustion kinetics is stiff where radical concentrations are near equilibrium and less so during ignition.

The stiffness does not have to come from the mathematics of the model. A grid discretization produces it. This is the essential point for anyone building a solver rather than studying one, and it connects to the fully discrete analysis below.

Why the eigenvalues bind an explicit method

For the test equation $y’ = \lambda y$ with $\lambda$ real and negative, forward Euler gives $y_{n+1} = (1 + h\lambda)y_n$. The mode is damped only when $|1 + h\lambda| \le 1$, that is, $-2 \le h\lambda \le 0$. With $\lambda = -1000$ the step cannot exceed $2 \times 10^{-3}$ no matter how slowly the solution of interest varies. If the quantity of interest changes on a scale of seconds, accuracy alone would allow a step of order 0.1; stability dictates a step two orders smaller.

Each extra digit of speed in the fastest mode divides the admissible step by the same factor while the solution of interest is untouched. The cost of an explicit integration over a fixed interval is therefore proportional to the stiffness ratio: the work is set by a transient that has long since died out, or, in the case of a mode that is never excited but is present in the Jacobian, by a component that contributes nothing at all.

Two mechanisms produce this in ordinary codes. A fine spatial grid discretizing a parabolic PDE has eigenvalues that scale as $-D/\Delta x^2$ with the diffusion coefficient $D$; refining the grid by two in each direction multiplies the fastest rate by four and forces a proportional reduction in the step, which is why explicit integration of diffusion is covered separately in explicit diffusion and the parabolic step limit. A reaction network with rates spanning five decades produces the same ratio from chemistry alone.

For complex eigenvalues the constraint is different in shape but not in kind. The explicit method’s region of absolute stability is a bounded set in the $h\lambda$ plane, and every eigenvalue must land inside it. The magnitude of the eigenvalues bounds the admissible step; their arrangement on the negative real axis versus off it decides which method is most economical.

Implicit steps and what they cost

An implicit method evaluates $f$ at the unknown state, which requires solving an equation. Backward Euler is

$$ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) $$

For the test equation this gives $y_{n+1} = y_n/(1 - h\lambda)$, which is damped for every negative $\lambda$ and, crucially, approaches zero as $h|\lambda| \to \infty$ rather than growing. The stability region contains the entire left half plane. That property is A-stability, and its consequences for the choice of method are the subject of A-stability and the Dahlquist barrier.

The price is a nonlinear solve per step. Newton’s method requires the Jacobian and the solution of a linear system with it:

/* One backward Euler step via Newton, dense n x n Jacobian J supplied. */
for (int iter = 0; iter < 10; iter++) {
    for (int i = 0; i < n; i++)
        r[i] = y_new[i] - y_old[i] - h * f(t_new, y_new, i);
    if (norm(r, n) < tol) break;               /* converged */
    for (int i = 0; i < n; i++)                /* M = I - h*J */
        for (int j = 0; j < n; j++)
            M[i][j] = (i == j) - h * J[i][j];
    solve_dense(M, r, n);                      /* r now holds the Newton update */
    for (int i = 0; i < n; i++) y_new[i] -= r[i];
}

The cost structure differs completely from the explicit case. Each step costs one or more Jacobian factorizations, $O(n^3)$ for a dense system of size $n$, plus triangular solves. The arithmetic is heavy and vectorizes poorly compared with a stencil sweep, so an implicit step on a three-dimensional grid may cost a hundred times an explicit one. The trade is favorable as soon as the stiffness ratio exceeds that factor, which for stiff kinetics is common.

The coupling inside the linear solve adds a second cost on parallel hardware. A triangular solve has a sequential dependency along the elimination order, while an explicit sweep over a grid has none, so the implicit step’s speedup from adding cores is limited by the factorization’s structure rather than by the arithmetic it performs. The degree to which the factorization can be reordered for parallelism is what makes the choice between a direct and an iterative inner solve a structural decision rather than a tuning parameter.

Sparsity is the real determinant

For a discretized PDE the Jacobian is sparse and the implicit step is not $O(n^3)$. A five-point stencil in two dimensions gives a Jacobian with five nonzero diagonals, and a banded or sparse factorization costs far less than a dense one. This is why implicit methods dominate PDE integration while explicit methods survive for problems whose Jacobians are dense but small, or whose stiffness is mild.

The structural fact is that the same stencil that produced the spatial operator reappears in the Jacobian. Building the Jacobian by finite differences of $f$ costs one function evaluation per column, which for a large stencil is optimistic; assembling it analytically from the known stencil is standard practice. The chapter on Newton iterations for the implicit spatial operator covers the assembly and the convergence conditions.

Is the answer an implicit discretization or a better model

Before reaching for an implicit solver, it is worth asking whether the fast mode should be in the model at all. If a species reaches equilibrium in microseconds and never departs, treating it as being in equilibrium at all times; setting its time derivative to zero and solving the resulting algebraic constraint; removes the mode from the equations. The result is a differential-algebraic system, smaller and non-stiff in the removed direction, and the step is limited by the remaining dynamics. This is the quasi-steady-state assumption, and its validity is a modelling question: it fails whenever the fast species is perturbed quickly, as happens during ignition or whenever the rates themselves change abruptly.

Similar reasoning applies to discretized PDEs. A grid that resolves a boundary layer much more finely than the interior dynamics require introduces eigenvalues far faster than the solution’s variation. A nonuniform grid or a different discretization can reduce the stiffness ratio substantially. Stiffness is sometimes a property of the model, sometimes of the discretization, and the distinction decides whether the fix belongs in the solver or in the setup.

There is also a hybrid worth knowing. IMEX schemes integrate the stiff linear or well-behaved part implicitly and the non-stiff nonlinear part explicitly, so the costly solve involves only the operator that needs it and the part that changes character along the solution stays cheap. The design of such partitions follows from the eigenvalue analysis above: the implicit part must cover the modes whose magnitudes force the step down.

The diagnostic that saves time

The failure signature of a non-stiff solver applied to a stiff problem is distinctive. The step controller reduces the step drastically for accuracy reasons while the step-size history shows the step oscillating: accepted, then rejected, while the solution’s variation per step stays small. If the solver also reports a step below what a smooth solution should require, or if shrinking the tolerance changes nothing about the accepted step, the Jacobian’s largest eigenvalue is the binding constraint and no amount of accuracy tuning will help. The measurement is straightforward when $f$ is available: estimate the spectral radius by a few power iterations on the Jacobian at a typical point, and compare it to the inverse of the current step. A ratio well above one is the definition of the trouble.