4.5. A-stability

The plots of the region of absolute stability for the explicit RK4 method and the implicit Radau IA method are shown in Fig. 4.4 and Fig. 4.5. We can see that the region of absolute stability for the explicit method is bounded in the left-hand side of the complex plane, whereas the region of absolute stability for the implicit Radau IA method is unbounded. This means that the step length used for the RK4 method needs to be less than some value in order for the method to be stable, whereas for the Radua IA method the step length can be any value. This property is known as A-stability.

Fig. 4.4 Region of absolute stability for the RK4 method

Fig. 4.5 Region of absolute stability for the Radua IA method

Definition 4.8 (A-stability)

A method is said to be A-stable if its region of absolute stability satisfies

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

i.e., the method is stable for all points in the left-hand side of the complex plane.

Theorem 4.1 (Conditions for A-stability)

Given an implicit Runge-Kutta method with a stability function of the form

\[ R(z) = \frac{P(z)}{Q(z)} \]

and define the \(E\)-polynomial function

(4.8)\[ E(y)=Q(iy)Q(-iy)-P(iy)P(-iy), \]

then the method is A-stable if and only if the following are satisfied

• Criterion A: All roots of \(Q(z)\) have positive real parts;

• Criterion B: \(E(y)\geq 0\) for all \(y\in \mathbb{R}\).

Example 4.3

The Radua IA method is defined by the following Butcher tableau

\[\begin{split} \begin{array}{c|cc} 0 & \frac{1}{4} & -\frac{1}{4} \\ \frac{2}{3} & \frac{1}{4} & \frac{5}{12} \\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} \end{split}\]

Determine whether this method is A-stable and plot the region of absolute stability.

---

Solution

We saw in Example 4.2 that the stability function for the Radau IA method is

\[ \begin{align*} R(z) = \frac{1 + \frac{1}{3}z + \frac{1}{16}z^2}{1 - \frac{2}{3}z + \frac{1}{6}z^2}. \end{align*} \]

Here \(Q(z) = 1 - \frac{2}{3}z + \frac{1}{6}z^2\) which has roots at \(z = 2 \pm \sqrt{2}\) which both have positive real parts so criterion A for A-stability is satisfied. Using equation (4.8)

\[\begin{split} \begin{align*} E(y) &= \left( 1 - \frac{2}{3} i y - \frac{1}{6}y^2 \right) \left( 1 + \frac{2}{3} i y - \frac{1}{6}y^2 \right) \\ & \qquad - \left( 1 + \frac{1}{3} i y - \frac{1}{16} y^2 \right) \left( 1 - \frac{1}{3}i y - \frac{1}{16}y^2 \right) \\ &= \left( 1 + \frac{2}{3}iy - \frac{1}{6}y^2 - \frac{2}{3} i y + \frac{4}{9}y^2 + \frac{1}{9} i y^3 - \frac{1}{6}y^2 - \frac{1}{9} i y^2 + \frac{1}{36} y^4 \right) \\ & \qquad \left( 1 - \frac{1}{3} i y - \frac{1}{16}y^2 + \frac{1}{3} i y + \frac{1}{9} y^2 - \frac{1}{48} i y^2 - \frac{1}{16} y^2 + \frac{1}{48} i y^3 + \frac{1}{256} y^4 \right) \\ &= \left( 1 + \frac{1}{9} y^2 + \frac{1}{36} y^4 \right) - \left( 1 - \frac{1}{72} y^2 + \frac{1}{256} y^4 \right) \\ &= \frac{1}{8} y^2 + \frac{55}{2304} y^4. \end{align*} \end{split}\]

Since \(y_2\) and \(y_4\) are positive for all \(y\) then \(E(y)>0\) so criterion B for A-stability is satisfied. Since both criterion A and B are satisfied then we can say that this is an A-stable method.

---

4.5.1. Code

The Python and MATLAB code used to determine the stability function for the IRK method from Example 4.3 is given below.

Python

import sympy as sp

def P(z):
    return 1 + 1/3 * z + 1/16 * z ** 2

def Q(z):
    return 1 - 2/3 * z + 1/6 * z ** 2

# Check roots of Q(z)
sp.pprint(sp.solve(Q(z)))

# Check E(y)
y = sp.symbols('y')
sp.pprint(sp.simplify(Q(1j * y) * Q(-1j * y) - P(1j * y) * P(-1j * y)))

MATLAB

P = @(z) 1 + 1/3 * z + 1/16 * z ^ 2;
Q = @(z) 1 - 2/3 * z + 1/6 * z ^ 2;

% Check roots of Q(z)
solve(Q(z))

% Check E(y) >= 0
syms y
simplify(Q(1i * y) * Q(-1i * y) - P(1i * y) * P(-1i * y))