All the methods for approximating the solution to initial-value problems have error terms that involve a higher derivative of the solution of the equation. If the derivative can be reasonably bounded, then the method will have a predictable error bound that can be used to estimate the accuracy of the approximation. Even if the derivative grows as the steps increase, the error can be kept in relative control, provided that the solution also grows in magnitude. Problems frequently arise, however, when the magnitude of the derivative increases but the solution does not. In this situation, the error can grow so large that it dominates the calculations. Initial-value problems for which this is likely to occur are called stiff equations.

Stiff differential equations are characterised as those whose exact solution has a term of the form \(e^{-ct}\), where \(c\) is a large positive constant. This is usually only a part of the solution, called the transient solution. The more important portion of the solution is called the steady-state solution. The transient portion of a stiff equation will rapidly decay to zero as t increases, but since the \(n^{th}\) derivative of this term has magnitude \(c^ne^{-ct}\), the derivative does not decay as quickly. In fact, since the derivative in the error term is evaluated not at \(t\) but at a number between zero and \(t\), the derivative terms can increase as \(t\) increases very rapidly.

Test Equations

Although stiffness is usually associated with systems of differential equations, the approximation characteristics of a particular numerical method applied to a stiff system can be predicted by examining the error produced when the method is applied to a simple test equation:

$$ \frac{dy}{dt} = \lambda y(t), \quad a \leq t \leq b, \quad y(a) = \alpha, \quad \lambda \lt 0 \\ \label{eq:stf1} $$

The solution to this equation is \(y(t) = \alpha e^{\lambda t}\), which contains the transient solution \(e^{\lambda t}\). The steady-state solution is zero, so the approximation characteristics of a method are easy to determine.

Consider a \(1^{st}\) order Taylor method applied to the test equation. The resulting difference-equation method may be written as:

\begin{equation} \begin{aligned} \omega_0 &= \alpha \\ \omega_{i+1} &= \omega_i + hf(t_i, \omega_i) = \omega_i + h \lambda \omega_i = (1 + h \lambda) \omega_i = (1 + h \lambda)^{i+1} \omega_0 = (1 + h \lambda)^{i+1} \alpha \quad i=1 \dots N-1 \\ \end{aligned} \label{eq:stf2} \end{equation}

Since the exact solution is \(y(t) = \alpha e^{\lambda t}\), the absolute error is:

$$ |y(t_i) - \omega_i| = |e^{ih \lambda} - (1 + h \lambda)^i||\alpha| = |(e^{\lambda h})^i - (1 + h \lambda)^i||\alpha| \\ \label{eq:stf3} $$

and the accuracy is determined by how well the term \((1 + h \lambda)\) approximates \(e^{\lambda h}\).

The exact solution decays to zero when \(\lambda \lt 0\) as \(i\) increases. The difference-equation approximation will have this property only if \(|1 + h \lambda| \lt 1\). This condition restricts the step size for the \(1^{st}\) order Taylor method to \(h \lt \frac{2}{|\lambda|}\).

In general, a function \(Q\) exists with the property that the difference method, when applied to the test equation, gives:

$$ \omega_{i+1} = Q(h \lambda) \omega_i \\ \label{eq:stf4} $$

The accuracy of the method depends upon how well \(Q(h \lambda)\) approximates \(e^{\lambda h}\), and the error will grow without bound if \(|Q(h \lambda)| \gt 1\). An \(n^{th}\) order Taylor method will have stability with regard to the absolute error, provided \(h\) is chosen to satisfy:

$$ \left|1 + h \lambda + \frac{1}{2}h^2 \lambda^2 + \dots + \frac{1}{n!}h^n \lambda^n \right| \lt 1\\ \label{eq:stf5} $$

When a multistep method is applied to the test equation, the resulting difference equation may be written as:

\begin{equation} \begin{aligned} \omega_{i+1} &= a_{m-1} \omega_i + \dots + a_0 \omega_{i+1-m} + h \lambda (b_m \omega_{i+1} + b_{m-1} \omega_i + \dots + b_0 \omega_{i+1-m}) \\ 0 &= (1 - h \lambda b_m) \omega_{i+1} - (a_{m-1} + h \lambda b_{m-1}) \omega_i - \dots - (a_0 + h \lambda b_0) \omega_{i+1-m} \\ \end{aligned} \label{eq:stf6} \end{equation}

for \(i=m-1 \dots N-1\). Associated with this homogeneous difference equation is a characteristic polynomial:

$$ Q(z,h \lambda) = (1 - h \lambda b_m) z^m - (a_{m-1} + h \lambda b_{m-1}) z^{m-1} - \dots - (a_0 + h \lambda b_0) \\ \label{eq:stf7} $$

Suppose \(\omega_0 \dots \omega_{m-1}\) are given, and, for fixed \(h \lambda\), let \(\beta_1 \dots \beta_m\) be zeros of the polynomial \(Q(z,h \lambda)\). Assuming \(\beta_1 \dots \beta_m\)are distinct, then \(c_1 \dots c_m\) exist with:

$$ \omega_i = \sum_{k=1}^{m} c_k (\beta)^i \quad i=1 \dots N \\ \label{eq:stf8} $$

If \(Q(z,h \lambda)\) has multiple zeros, \(\omega_i\) is similarly defined. If \(\omega_i\) is to accurately approximate \(y(t) = \alpha e^{\lambda t}\) then all zeros \(\beta_k\) must satisfy \(|\beta_k| \lt 1\) otherwise, certain choices of \(\alpha\) will result in \(c_k \ne 0\) and the term \(c_k(\beta_k)^i\) will not decay to zero.

The above analysis allows the region of absolute stability for a difference-equation method to be defined. The region \(R\) of absolute stability for a one step method is \(R=\{h \lambda \in C | \quad |Q(h \lambda)| \lt 1\}\) and for a multistep method, it is \(R=\{h \lambda \in C | \quad |\beta_k| \lt 1\}\).