Stiff sets of differential equations

7.4 Stiff sets of differential equations

As soon as one deals with more than one first-order differential equation, the possibility of a stiff set of equations arises. Stiffness occurs in a problem where there are two or more very different scales of the independent variable on which the dependent quantities are varying. For example, consider the following set of equations:

\begin{equation*} \label{StiffEqua1} \begin{split} u'&=998u+1998v\\ v'&=-999u-1999v \end{split} \end{equation*}

(7.9)

with boundary conditions

\begin{equation*} u(0)=1, \quad v(0)=0. \end{equation*}

(7.10)

By means of the transformation

\begin{equation*} u=2y-z, \quad v=-y+z, \end{equation*}

(7.11)

one finds the solution

\begin{equation*} \begin{split} u&=2e^{-x}-e^{-1000x},\\ v&=-e^{-x}+e^{-1000x}. \end{split} \end{equation*}

(7.12)

If one integrated the system (7.9) with any of the methods given so far in this chapter, the presence of the \(e^{-1000x}\) term would require a step size \(h \ll 1/1000\) for the method to be stable—even though the \(e^{-1000x}\) term is completely negligible in determining the values of \(u\) and \(v\) as soon as one is away from the origin (see Fig. 7.3). This is the generic disease of stiff equations: one is required to follow the variation in the solution on the shortest length scale to maintain stability of the integration, even though accuracy requirements allow a much larger step size.

Figure 7.3: Example of an instability encountered in integrating a stiff equation (schematic). Here it is supposed that the equation has two solutions, shown as solid and short-dotted lines. Although the initial conditions are such as to give the solid solution, the stability of the integration (shown as the unstable dotted sequence of segments) is determined by the more rapidly varying dashed solution, even after that solution has effectively died away to zero.

To see how this problem might be cured, consider the single equation

\begin{equation*} y'=-cy \end{equation*}

(7.13)

where \(c\gt 0\) is a constant. The usual (or forward) Euler scheme for integrating this equation with step size \(h\) is

\begin{equation*} y_{n+1}=y_n+hy'_n=y_n+h(-cy_n)=(1-ch)y_n. \end{equation*}

(7.14)

The method is also called explicit because the new value \(y_{n+1}\) is given explicitly in terms of the old value \(y_n\). Clearly the method is unstable if \(h\gt 2/c\), for then \(\left|y_n\right|\rightarrow\infty\) as \(n\rightarrow\infty\).

The simplest cure is to resort to implicit differencing, where the right-hand side is evaluated at the new \(y\) location. In this case, one obtains the stable backward Euler scheme:

\begin{equation*} y_{n+1}=y_n+hy'_{n+1} \end{equation*}

(7.15)

\begin{equation*} y_{n+1}=\frac{y_n}{1+ch}. \end{equation*}

(7.16)