8.2. Stiff Systems of Differential Equations

Although there are situations when a single differential equation is stiff, in real life problems, stiffness is mainly a phenomenon related to systems of differential equations. A stiff system of differential equation, similarly, has a transient solution as well as a slowly varying or steady-state solution.

Consider a general system of \(\,n\,\) ordinary differential equations of the form

\[\begin{aligned} \mathbf{y}' ~=~ \mathbf{f} (x,\mathbf{y}) \end{aligned}\]

given by

\[\begin{aligned} {\bf y^\prime} ~=~ {\bf \Lambda\,y} + {\bf g}(x). \end{aligned}\]

where \(\,\bf y\,\) and \(\,\bf g\,\) are vectors with \(\,n\,\) components and \(\,\bf \Lambda\,\) is the Jacobian matrix, J, of \(\,\bf f\,\) (i.e. \(\,{\bf f}_y(x_j,y_j)\,\) at some local point, or at some step, \(\,j\,\)). For example, for a system of three equations:

\[\begin{split}\begin{aligned} y'_1 ~&=~ f_1(x,y_1,y_2,y_3) \\[1ex] y'_2 ~&=~ f_2(x,y_1,y_2,y_3) \\[1ex] y'_3 ~&=~ f_3(x,y_1,y_2,y_3) \end{aligned}\end{split}\]

the associated Jacobian matrix is

\[\begin{split}\begin{aligned} \mathbf{J} = \begin{pmatrix} \dfrac{\partial f_1}{\partial y_1} & \dfrac{\partial f_1}{\partial y_2} & \dfrac{\partial f_1}{\partial y_3}\\[4pt] \dfrac{\partial f_2}{\partial y_1} & \dfrac{\partial f_2}{\partial y_2} & \dfrac{\partial f_2}{\partial y_3}\\[4pt] \dfrac{\partial f_3}{\partial y_1} & \dfrac{\partial f_3}{\partial y_2} & \dfrac{\partial f_3}{\partial y_3}\\[4pt] \end{pmatrix} \end{aligned}\end{split}\]

Here, the eigenvalues \(\,\lambda_i\,\) of the Jacobian matrix play the part of \(\,\lambda\,\). The system of equations is stiff if there is an eigenvalue, \(\,\lambda_i\,\), of \(\mathbf{f}\), with a negative real part of a very large magnitude. In a system of equations there are several \(\,\lambda's\,\) with corresponding time constants,

\[\begin{aligned} && \tau_i ~&=~ \dfrac{1}{|Real\,(\lambda_i)|} & i &= 1,2,...,n && \end{aligned}\]

In a system of equations the larger the magnitude difference between the eigenvalues of the corresponding Jacobian matrix, the higher is the degree of stiffness. The degree of stiffness can then be measured by calculating the stiffness ratio, \(\,S(x)\,\), given by

\[\begin{aligned} && S(x) ~&=~ \dfrac{\max{|Real\,(\lambda_i)|}}{\min{|Real\,(\lambda_i)|}} & i &= 1,2,...,n && \end{aligned}\]

The higher the stiffness ratio the higher the degree of stiffness of the system.

Note

The above system is a linear system and the terms in the Jacobian, \(\dfrac{\partial \mathbf{f}}{\partial \mathbf{y}}\,\), and thus the corresponding eigenvalues are constant throughout the integration. For non-linear systems the terms in the Jacobian, \(\dfrac{\partial \mathbf{f}}{\partial \mathbf{y}}\,\), and thus the corresponding eigenvalues can be a function of some variables in the system and they are calculated locally as the integration proceeds.

Example 8.3

Consider the non-linear system

\[\begin{split}\begin{aligned} f_1 ~&=~ -0.013y_1 + 10^3y_1y_3\\ f_2 ~&=~ +2.5\times10^3y_2y_3\\ f_3 ~&=~ 0.013y_1 - 10^3y_1y_3 - 2.5\times10^3y_2y_3 \end{aligned}\end{split}\]

with initial conditions

\[\begin{split}\begin{aligned} \mathbf{y}(0) ~=~ \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}. \end{aligned}\end{split}\]

This is an example from chemical reaction kinetics and \(y_3\) represents the concentration of a very reactive species which is an intermediate in the course of reaction and always stays small. \(y_1\) and \(y_2\) are monotonically decreasing and increasing respectively. \(y_3\) increases to a maximum and thereafter is monotonically decreasing. It can be shown that \(\,y_3<1.3\times10^{-5}\). The Jacobian of the system is given by

\[\begin{split}\begin{aligned} \mathbf{J} ~=~ \begin{pmatrix} -0.013+10^3y_3 & 0 & \phantom{-}10^3y_1\phantom{~-2.5\times10^3y_2} \\ 0 & \phantom{-}2.5 \times10^3y_3 & \phantom{-10^3y_1-~}2.5\times10^3y_2 \\ \phantom{-}0.013-10^3y_3 & -2.5 \times10^3 y_3 & -10^3y_1-2.5 \times10^3 y_2 \end{pmatrix} \end{aligned}\end{split}\]

Since \(y_3\) is small, the last column has elements which are much larger than the remainder - the large magnitude differences between the columns is a characteristic of stiff systems. Note that at any time step \(j\), the calculated values of \(y_1\), \(y_2\), and \(y_3\) can be used for calculating the corresponding values of \(\lambda_1\), \(\lambda_2\), and \(\lambda_3\).

Example 8.4

Consider the differential equation system

\[\begin{split}\begin{aligned} && u' ~&=~ \phantom{-}998u + 1998v & u(0) &= 1 && \\ && v' ~&=~ -999u - 1999v & v(0) &= 0 && \end{aligned}\end{split}\]

The analytic solution is

\[\begin{split}\begin{aligned} u ~&=~ 2e^{-t} - e^{-1000t}\\ v ~&=~ -e^{-t} + e^{-1000t}\\ \end{aligned}\end{split}\]

where \(t\) is the independent variable.

Analytical Solution to Linear differential equation system

For a homogeneous (linear) system of differential equations, we can use the eigenvalues and eigenvectors of the system to obtains its analytical solution.

Supposing we have a linear system of differential equations given by

(8.1)\[ \mathbf{U}'=\mathbf{A}\mathbf{U}, \]

and assuming its solution can be expressed as

(8.2)\[ \mathbf{U}=\mathbf{v}e^{\lambda t}, \]

where \(\mathbf{A}\) is the matrix of the linear system, and \(\mathbf{v}\) is a vector. Differentiating Eq. (8.2), we can obtain

(8.3)\[ \mathbf{U}'=\lambda\mathbf{v}e^{\lambda t}. \]

Substituting Eqs. (8.2) and (8.3) into Eq. (8.1), we can obtain

(8.4)\[ \mathbf{A}\mathbf{v}e^{\lambda t} = \lambda \mathbf{v}e^{\lambda t}, \]

dividing \(e^{\lambda t}\) gives:

\[ \mathbf{A}\mathbf{v} = \lambda \mathbf{v}. \]

This means \(\lambda\) and \(\mathbf{v}\) are the eigenvalue and eigenvector of the system, respectively. Once we obtain the \(\lambda\) and \(\mathbf{v}\), we can obtain the analytical solution to the system (8.1) given by

\[ \mathbf{U}=\sum_{i=1}^n \alpha_i \mathbf{v}_i e^{\lambda_i t}, \]

where \(n\) is the size of the system, \(\alpha_i (i=1,2,\ldots,n)\) are coefficients need to be determined by using the initial condition.

Regarding Example 8.4, we have

\[\begin{split} \mathbf{A}=\mathbf{J}= \begin{bmatrix} 998 & 1998 \\ -999 & -1999 \end{bmatrix}, \end{split}\]

its eigenvalues are

\[ \lambda_1=-1, \quad \lambda_2=-1000, \]

its eigenvectors are given by

\[\begin{split} \mathbf{v}_1=\begin{bmatrix} 2 \\ -1 \end{bmatrix}, \quad \mathbf{v}_2=\begin{bmatrix} 1 \\ -1 \end{bmatrix}. \end{split}\]

Now the solution to the system can be expressed as

\[\begin{split} \begin{aligned} \mathbf{U} & = \alpha_1 \mathbf{v}_1 e^{\lambda_1 t} + \alpha_2 \mathbf{v}_2 e^{\lambda_2 t} \\ & = \alpha_1 \begin{bmatrix} 2 \\ -1 \end{bmatrix} e^{-t}+ \alpha_2 \begin{bmatrix} 1 \\ -1 \end{bmatrix} e^{-1000t}. \end{aligned} \end{split}\]

At \(t=0\), \(\mathbf{U}(0)=[1, 0]^{\mathrm{T}}\), so we obtain

\[ \alpha_1 = 1, \quad \alpha_2=-1. \]

Therefore the analytical solution is given by

\[\begin{split} \mathbf{U} = \begin{bmatrix} 2 \\ -1 \end{bmatrix} e^{-t} - \begin{bmatrix} 1 \\ -1 \end{bmatrix} e^{-1000t}= \begin{bmatrix} 2e^{-t}- e^{-1000t} \\ -e^{-t}+ e^{-1000t} \end{bmatrix}. \end{split}\]

The Jacobian \(\mathbf{J}\) of this system is:

\[\begin{split}\begin{aligned} \begin{pmatrix} \phantom{-}998 & \phantom{-}1998 \\ -999 & -1999 \end{pmatrix} \end{aligned}\end{split}\]

The corresponding eigenvalues are \(\,\lambda_1 = -1\,\) and \(\,\lambda_2 = -1000\,\). Thus, the stiffness ratio is \(1000\) which is not too high; in practical problems a ratio of the order of \(10^{6\pm2}\) is not unusual. The stiffness phenomenon of this problem can also be observed from the analytical solution; both dependent variables contain both fast and slow components. Depending which numerical method is chosen for the solution of this system there will be restrictions on \(h\) even when the rapidly changing component is insignificant.

Note

Note that the above problem is a linear system and the terms of the Jacobian, \(\,\dfrac{\partial \mathbf{f}}{\partial \mathbf{y}}\,\), and thus the corresponding eigenvalues are constant throughout the integration!