Stiffness, Stability, and Implicit Methods

Definition An initial-value problem is stiff on an interval when some of its modes decay much faster than the solution being followed changes, so that stability, rather than the accuracy needed for that solution, forces an explicit method to take very small steps. For a linear system $X' = AX + g(t)$ with a constant matrix $A$, the fast modes belong to eigenvalues $\lambda$ of $A$ with $\operatorname{Re}\lambda < 0$ and $\lvert \operatorname{Re}\lambda \rvert$ large compared with the rates at which the solution itself varies. For a nonlinear system $X' = f(t, X)$, the eigenvalues of the Jacobian matrix $\partial f/\partial X$ along the solution play the same role as a local guide. Stiffness is a practical judgment that depends on the interval, the accuracy required and the method, not a property that one formula decides.

Definition The test equation is $y' = \lambda y$ with a constant $\lambda$, real or complex. Its nonzero solutions $y(t) = y(0)e^{\lambda t}$ decay to $0$ exactly when $\operatorname{Re}\lambda < 0$, and one step of length $h$ multiplies them by $e^{h\lambda}$. A Runge-Kutta method, including explicit and implicit Euler, applied to the test equation with step $h$ produces

where the function $R$ depends only on the method. The number $R(h\lambda)$ is the amplification factor of one step, and $z = h\lambda$ is the scaled eigenvalue. The test equation stands for one mode of a linear system. If $X' = AX$ with $A = VDV^{-1}$, where $D$ is diagonal with entries $\lambda_1, \dots, \lambda_n$, the substitution $X = Vw$ splits the system into $w_i' = \lambda_i w_i$, and these methods treat each $w_i$ exactly as they treat the test equation with $\lambda = \lambda_i$. For instance, Euler gives $X_{n+1} = (I + hA)X_n$, which multiplies each $w_i$ by $1 + h\lambda_i$.

Definition The region of absolute stability of a one-step method with amplification factor $R$ is the set of complex numbers $z$ with $\lvert R(z) \rvert < 1$. When $z = h\lambda$ lies in it, the numerical solution of the test equation tends to $0$ as $n \to \infty$ for every starting value. Some books also include the boundary points where $\lvert R(z) \rvert = 1$. For Euler’s method the region is the open disk

of radius $1$ centered at $-1$. For a real $\lambda < 0$ this means $0 < h < 2/\lvert \lambda \rvert$. For every explicit Runge-Kutta method, $R$ is a polynomial that is not constant, so $\lvert R(z) \rvert \to \infty$ as $\lvert z \rvert \to \infty$ and the region is bounded.

Definition For $y' = f(t, y)$ with step $h$ and $t_{n+1} = t_n + h$, the implicit Euler method, also called backward Euler, defines $y_{n+1}$ as a solution of

The slope is taken at the new point, so each step requires solving an equation for the unknown $y_{n+1}$. If $f$ is Lipschitz in $y$ with constant $L$ for all $y$ and $hL < 1$, this equation has exactly one solution, and simple fixed-point iteration converges to it. A stiff problem has a large $L$, and that condition would bring back the small steps, so stiff solvers usually solve the step equation by Newton’s method with the Jacobian $\partial f/\partial y$. On smooth problems the global error is of order $h$, the same order as explicit Euler. On the test equation the method gives $y_{n+1} = y_n + h\lambda y_{n+1}$, so for $h\lambda \ne 1$

Definition A method is A-stable if its region of absolute stability contains every complex $z$ with $\operatorname{Re} z < 0$. Then, for the test equation with any $\lambda$ satisfying $\operatorname{Re}\lambda < 0$, the numerical solution tends to $0$ for every step $h > 0$. The implicit Euler method is A-stable. If $\operatorname{Re} z < 0$, then

so $\lvert R(z) \rvert = 1/\lvert 1 - z \rvert < 1$. Its region is in fact the whole exterior $\lvert z - 1 \rvert > 1$ of the closed disk of radius $1$ centered at $1$. No explicit Runge-Kutta method is A-stable, because its region is bounded. A-stability is defined through the linear test equation, so for a nonlinear problem it describes the behavior near a solution through linearization and is not by itself a proof of stability.

Definition Physical stability is a property of the differential equation. A solution is stable if perturbations of it stay small, and asymptotically stable if they also decay to $0$. For the test equation the perturbations decay exactly when $\operatorname{Re}\lambda < 0$. Numerical stability is a property of a method used with a particular step on a particular problem. For the test equation it means $\lvert R(h\lambda) \rvert < 1$, so that perturbations of the computed values decay. Accuracy is a third question, how close the computed values are to the true solution, and it is judged from the method’s order, from error estimates and from refinement. A computation can pass any one of these tests and fail another.