4.2. 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 \(\tau_{n+1}\) are also updated using \(\tau_n\) by the same method. This allows us to define a stability function for a method.

Definition 4.3 (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\).

4.2.1. Absolute stability

We have seen that a necessary condition for stability of a method is that the local truncation errors must not grow from one step to the next. A method satisfying this basic condition is considered to be absolutely stable. Since the stability function \(R(z)\) is expressed using \(z=h\lambda\) then a method may be stable for some value of \(h\) and unstable for others. This provides the definition for absolute stability.

Definition 4.4 (Absolute stability)

A method is considered to be absolutely stable if \(|R(z)| \leq 1\) for \(z\in \mathbb{C}\).

Of course we require our methods to be stable so it is useful to know for what values of \(h\) we have a stable method. This gives the definition of the region of absolute stability.

Definition 4.5 (Region of absolute stability)

The region of absolute stability is the set of \(z\in \mathbb{C}\) for which a method is absolutely stable

\[ \{ z:z\in \mathbb{C},|R(z)|\leq 1 \} \]

The region of absolute stability for the Euler method is shown in Fig. 4.3.

Fig. 4.3 The region of absolute stability for the Euler method.

4.2.2. Plotting stability regions

We can plot the region of absolute stability by generate a set of points for \(z\) in the complex plane and plot the contour where \(|R(z)| = 1\) which is the boundary of the stability region. The code for producing a plot of the region of absolute stability of the Euler method using Python is shown below.

Python

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True  # allows use of LaTeX in labels

# Generate z values
xmin, xmax, ymin, ymax = -3, 1, -1.5, 1.5
X, Y = np.meshgrid(np.linspace(xmin, xmax, 200),np.linspace(ymin, ymax, 200))
Z = X + Y * 1j

# Define stability function
R = 1 + Z

# Plot stability region
fig = plt.figure()
contour = plt.contourf(X, Y, abs(R), levels=[0, 1], colors="#99ccff")  # Plot stability region
plt.contour(X, Y, abs(R), colors= "k", levels=[0, 1])                  # Add outline
plt.axhline(0, color="k", linewidth=1)                                 # Add x-axis line
plt.axvline(0, color="k", linewidth=1)                                 # Add y-axis line
plt.axis("equal")
plt.axis([xmin, xmax, ymin, ymax])
plt.xlabel("$\mathrm{Re}(z)$", fontsize=12)
plt.ylabel("$\mathrm{Im}(z)$", fontsize=12)
plt.show()

MATLAB

% Generate z values
xmin = -3;
xmax = 1;
ymin = -1.5;
ymax = 1.5;
[X, Y] = meshgrid(linspace(xmin, xmax, 200), linspace(ymin, ymax, 200));
Z = X + Y * 1i;

% Define stability function 
R = 1 + Z;

% Plot stability region
contourf(X, Y, abs(R), [0, 1], LineWidth=2)
xline(0, LineWidth=2)
yline(0, LineWidth=2)
colormap([153, 204, 255 ; 255, 255, 255] / 255)
axis equal
axis([xmin, xmax, ymin, ymax])
xlabel("$\mathrm{Re}(z)$", FontSize=12, Interpreter="latex")
ylabel("$\mathrm{Im}(z)$", FontSize=12, Interpreter="latex")

4.2.3. Interval of absolute stability

The choice of step length used in a method will depend on accuracy requirements, the computational resources available, the stability properties of the method and the ODE being solved. It is often necessary to use as large a value of the step length as possible permitted by the stability requirements to minimise the computational effort required to solve an ODE. The range values of the step length that can be chosen is governed by the stability region and provides use with the following definition.

Definition 4.6 (Interval of absolute stability)

The range of real values that the step length \(h\) of a method can take that ensures a method remains absolutely stable is known as the interval of absolute stability

\[ \begin{align*} \{ z : z \in \mathbb{R}, |R(z)| \leq 1 \} \end{align*} \]

The interval of absolute stability for the Euler method can be seen in Fig. 4.4.

Fig. 4.4 Interval of absolute stability for the Euler method.

The region of absolute stability for the Euler method shows that the interval of absolute stability for the test ODE \(y'=\lambda y\) is

\[ \begin{align*} z \in [-2,0], \end{align*} \]

Since \(z = h\lambda\) then

\[ \begin{align*} h \in \left[ -\frac{2}{\lambda},0 \right], \end{align*} \]

so we have the condition

\[ \begin{align*} h \leq -\frac{2}{\lambda}. \end{align*} \]

We saw in Fig. 4.1 that solving the ODE \(y' = -15y\) using \(h=0.25\) resulted in an unstable solution whereas using \(h=0.125\) resulted in a stable (albeit oscillatory) solution. This is because for this ODE \(\lambda = -15\) and the step length for the Euler method must satisfy

\[ \begin{align*} h \leq -\frac{2}{-15} = 0.1\dot{3}. \end{align*} \]

This is why the solution using \(h=0.25\) was unstable since \(0.25 > 0.1\dot{3}\) and the solution using \(h=0.125\) was stable since \(0.125 < 0.1\dot{3}\).