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. Stability functions

To examine the behaviour of the local truncation errors we use the following test equation

(4.2)\[ y' = \lambda y. \]

This test equation is chosen as it is a simple ODE which has the solution \(y=e^{\lambda t}\). Over a single step of a Runge-Kutta method, the values of \(y_{n+1}\) are updated using the values of \(y_n\), so the values of the local truncation error, \(e_{n+1}\), are also updated using \(e_n\) by the same method. This allows us to define a stability function for a method.

Definition 4.2 (Stability function)

The stability function of a method, denoted by \(R(z)\), is the rate of growth over a single step of the method when applied to calculate the solution of an ODE of the form \(y'=\lambda t\) where \(z = h \lambda\) and \(h\) is the step size

\[ y_{n+1} = R(z)y_n. \]

For example, if the Euler method is used to solve the test equation \(y'=\lambda y\) then the solution will be updated over one step using

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

Let \(z = h\lambda\) then

\[ y_{n+1} = y_n + z y_n =(1+z)y_n. \]

So the stability function of the Euler method is \(R(z) = 1 + z\).