4. Stability

The stability of a numerical solver for ordinary differential equations refers to its ability to produce accurate and reliable solutions over a long interval, without the computed values of the solution becoming unbounded or oscillatory. In other words, a stable numerical solver for ODEs is one that does not produce solutions that blow up or oscillate wildly as time goes on.

For example, consider the Euler method solution to the following initial value problem using step lengths of \(h = 0.25\) and \(h = 0.125\)

\[ y' = -15y, \qquad t \in [0, 1], \qquad y(0) = 1. \]

The exact solution for this initial value problem is \(y = e^{-15t}\) and the numerical solutions are plotted against the exact solution in Fig. 4.1. Here we can see that the solution using \(h = 0.125\) remains bounded whereas the solution using \(h = 0.25\) diverges.

Fig. 4.1 Solutions to the initial value problem \(y' = -15y\), \(t\in [0, 1]\) and \(y(0) = 1\) using the Euler method with \(h=0.25\) and \(h=0.125\).

The solution using \(h = 0.125\) is an example of a stable solution where the numerical solution remains bounded whereas the solution using \(h = 0.25\) is an example of an unstable solution which is unbounded. The reason for the instability is due to the truncation of the Taylor series. We saw in the section on error analysis that the higher order terms that are omitted due to truncation causes a local truncation error, and the accumulation of these local truncation errors creates the global truncation error. If the local truncation error is too large then the global truncation error will eventually diverge as seen in Fig. 4.1 causing the solution to be unusable. This provides us with the definition of stability.

Definition 4.1 (Stability)

If \(e_n\) is the local truncation error of a numerical method defined by

\[\begin{align*} e_n = |y(t_{n}) - y_{n}|, \end{align*}\]

where \(y(t_n)\) is the exact solution and \(y_n\) the numerical approximation assuming that the solution at the previous step was exact. The method is considered stable if

(4.1)\[ |e_{n+1} - e_n | \leq 1, \]

for all steps of the method. If a method is not stable then it is unstable.

---

4.1. Stiffness

When solving IVPs numerically, some problems behave well and they can be solved accurately using large values of the step length \(h\). Other problems force us to use small step lengths use they become unstable. These problems as known as stiff problems.

An informal definition of stiffness is provided by Lambert [1991]

“If a numerical method with a finite region of absolute stability, applied to a system with any initial conditions, is forced to use in a certain interval of integration a step length which is excessively small in relation to the smoothness of the exact solution in that interval, then the system is said to be stiff in that interval.”

4.1.1. Linear systems

Given a linear system \(\mathbf{y}' = A \mathbf{y}\) then the stiffness ratio is defined as

\[ S= \frac{\max_i (Re(\lambda_i))}{\min_i (Re(\lambda_i))},\]

where \(\lambda_i\) are the eigenvalues of the matrix \(A\). If the value of \(S\) is large, then the system is stiff. For example, consider the system

\[\begin{split} \begin{align*} y_1' &= -1000y_1 + 999y_2, \\ y_2' &= -y_2. \end{align*} \end{split}\]

Here \(A = \begin{pmatrix} -1000 & 999 \\ 0 & -1 \end{pmatrix}\), so the eigenvalues satisfy

\[\begin{split} \begin{align*} 0 &= | A - \lambda I | \\ &= \begin{vmatrix} -1000 - \lambda & 999 \\ 0 & -1 - \lambda \end{vmatrix} \\ &= (-1000-\lambda) (-1 - \lambda). \end{align*} \end{split}\]

So the eigenvalues are \(\lambda_1 = -1000\) and \(\lambda_2 =-1\) and the stiffness ratio is

\[ S = \frac{\lambda_1}{\lambda_2} = \frac{-1000}{-1} = 1000.\]

Since \(S\) is large then this is a stiff system.

4.1.2. Nonlinear systems

Given a nonlinear system \(\mathbf{y}' = \mathbf{f}(t, \mathbf{y})\) then matrix \(A\) is replaced by the Jacobian matrix at a particular point \((t, \mathbf{y})\)

\[\begin{split} J(t, \mathbf{y}) = \frac{\partial \mathbf{f}}{\partial y} = \begin{pmatrix} \dfrac{\partial f_1}{\partial t_1} & \cdots & \dfrac{\partial f_1}{\partial t_N} \\ \vdots & \ddots & \vdots \\ \dfrac{\partial f_N}{\partial t_1} & \cdots & \dfrac{\partial f_N}{\partial t_N} \end{pmatrix}. \end{split}\]

the stiffness ratio is then defined locally by

\[ S(t, \mathbf{y}) = \frac{\max_i (Re(\lambda_i))}{\min_i (Re(\lambda_i))}.\]

For example, consider the nonlinear system

\[\begin{split} \begin{align*} y_1' &= -1000y_1, \\ y_2' &= -y_2^3 \end{align*} \end{split}\]

The Jacobian matrix is

\[\begin{split} J = \begin{pmatrix} -1000 & 0 \\ 0 & -3y_2^2 \end{pmatrix}. \end{split}\]

The eigenvalues satisfy

\[\begin{split} \begin{align*} 0 &= |J - \lambda I| \\ &= \begin{vmatrix} -1000 - \lambda & 0 \\ 0 & -3y_2^2 - \lambda \end{vmatrix} \\ &= (-1000 - \lambda)(-3y_2^2 - \lambda), \end{align*} \end{split}\]

so \(\lambda_1 = -1000\) and \(\lambda_2 = -3y_2^2\). The stiffness ratio is therefore

\[ S(y_2) = \frac{1000}{3y_2^2}. \]

Unlike the linear case, the stiffness ratio depends on the current value of \(y_2\), e.g.,

\[\begin{split} \begin{align*} S(y_2 = 1) &= 333.3, \\ S(y_2 = 0.1) &= 33333.3. \end{align*} \end{split}\]

As \(y_2\)​ approaches zero, the stiffness ratio increases dramatically. This shows that the stiffness of a nonlinear system can change as the solution evolves, so a system that is only mildly stiff at one point may become extremely stiff later.