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, \(\tau\), 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 \(\tau_n\) is the local truncation error of a numerical method defined by

\[\begin{align*} \tau_n = |y_n - \tilde{y}_n|, \end{align*}\]

where \(\tilde{y}_n\) is the exact solution and \(y_n\) the numerical approximation then the method is considered stable if

(4.1)\[ |\tau_{n+1} - \tau_n | \leq 1, \]

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

4.1. Stiffness

The stability of a numerical solver for ODEs is closely related to the stiffness of the ODEs being solved. Stiff ODEs are those where the solution varies over several orders of magnitude or changes rapidly over a small interval, compared to the domain of the ODEs. In such cases, standard numerical solvers may produce inaccurate or unstable solutions if the time step used in the computation is not small enough to capture the fast dynamics of the problem.

There is no mathematical definition of what makes an ODE stiff, however, a useful description 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. Stiffness ratio

Stiffness is not strictly defined but we can attempt to define a metric to help us determine whether an ODE is stiff. Consider a system of linear ODEs for the form

\[ \vec{y}' = A \vec{y}, \]

where \(\vec{y} = (y_1, y_2, \ldots, y_n)^T\) and \(A\) is an \(n \times n\) matrix of constants. Let \(V = (\vec{v}_1, \vec{v}_2, \ldots, \vec{v}_n)\) be a matrix containing the eigenvalues of \(A\) then

\[V = (A \vec{v}_1, A \vec{v}_2, \ldots, A \vec{v}_n).\]

Since \(A \vec{v}_i = \lambda_i \vec{v}_1\) where \(\lambda_i\) is the eigenvalue corresponding to the eigenvector \(\vec{v}_i\) then

\[ AV = (\lambda_1 \vec{v}_1, \lambda_2 \vec{v}_2, \ldots, \lambda_n \vec{v}_n) = V \Lambda \]

where \(\Lambda\) is a diagonal matrix of the eigenvalues \(\lambda_i\) (take care to make the distinction between \(\Lambda\) which is the the uppercase \(\lambda\) character and A). Multiplying both sides by \(V^{-1}\) gives

\[ V^{-1} A V = \Lambda, \]

which means that \(V^{-1} A\) is a linear transformation that that diagonalises the matrix \(A\). Let \(\vec{u} = V^{-1} A \vec{y}\) then we have the transformed system \(\vec{u}' = \Lambda \vec{u}\) which is

\[\begin{split} \begin{align*} \begin{pmatrix} u_1' \\ u_2' \\ \vdots \\ u_n' \end{pmatrix} = \begin{pmatrix} \lambda_1 \\ & \lambda_2 \\ && \ddots \\ &&& \lambda_n \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{pmatrix} = \begin{pmatrix} \lambda_1 u_1 \\ \lambda_2 u_2 \\ \vdots \\ \lambda_n u_n \end{pmatrix}. \end{align*} \end{split}\]

Each equation in the transformed system has the solution \(u_i = e^{\lambda_i t}\) so the fastest growing solution will be for the equation with the largest value of \(\lambda\) and the slowest growing solution will be for the equation with the smallest value of \(\lambda\). Stiffness occurs when there is a large variation in the solutions of equations in a system so we can gauge this by calculating the ratio between the largest and smallest eigenvalues of \(A\).

Definition 4.2 (Stiffness ratio)

The stiffness ratio for a system of ODEs of the form \(\vec{y}' = A \vec{y}\) is

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

where \(\lambda_i\) are the eigenvalues of \(A\).

If the stiffness ratio of a system of ODEs is large[1] then the system is considered to be stiff.

Example 4.1

Determine the stiffness ratio of the following ODE

\[ y'' + 1001 y' + 1000 y = 0.\]

Solution (click to show)

Rewriting the second-order ODE as a system of two first-order ODEs

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

which can be written as the matrix equation

\[\begin{align*} \begin{pmatrix} y_1' \\ y_2' \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1000 & -1001 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}. \end{align*}\]

Calculating the eigenvalues of the coefficient matrix

\[\begin{align*} 0 &= \begin{vmatrix} -\lambda & 1 \\ -1000 & -1001 - \lambda \end{vmatrix} = \lambda^2 + 1001 \lambda + 1000 = (\lambda + 1000) (\lambda + 1), \end{align*}\]

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

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

Since \(S = 1000\) is not too large then this system is not considered to be stiff. For practical problems the stiffness ration of \(10^{6}\) can be common.

4.1.2. Solving stiff systems

Different ODE solvers have different stability properties and are stable over larger values of \(h\). For example the implicit Euler method is an A-stable method meaning that it is stable for any value of \(h\).

\[ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) \]

A comparison of the Euler and implicit Euler methods used to solve the same initial value problem as before is shown in Fig. 4.2 where the implicit Euler method solutions remain bounded despite using the same step length.

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

So if we were constrained with using a step length of \(h=0.25\) then we could choose the implicit Euler method to solve this initial value problem.

---

[1]

What constitutes “large” is vague and subjective. As a rule of thumb if \(S\) is in the hundreds or greater we can consider the system to be stiff.