7.4. The Boundary Root Locus Method for Finding Regions of Absolute Stability

We illustrate this method by considering the following example.

Example 7.4

Find the region of absolute stability for the two-step explicit Adams method

\[ \begin{aligned} y_{j+1} ~=~ y_j + \frac{h}{2}(3f_j - f_{j-1}) \end{aligned} \]

when applied to \(y' = \lambda y\).

Solution

The characteristic equation is

\[ \begin{aligned} z^2 - (1 + \frac{3}{2}h\lambda)z + \frac{1}{2}h\lambda ~=~ 0 \end{aligned} \]

thus

\[ \begin{aligned} \rho(z) ~=~ z^2 - z \qquad\text{and}\qquad \sigma(z) ~=~ \frac{1}{2}(3z - 1) \end{aligned} \]

Since \(\qquad \L(z, h\lambda) ~=~ \rho(z) - h\lambda\sigma(z) ~=~ 0 \qquad\Rightarrow\qquad h\lambda ~=~\dfrac{\rho(z)}{\sigma(z)}\)

and we get \(\displaystyle\qquad h\lambda ~=~ \frac{z^2 - z}{\frac{1}{2}(3z - 1)} ~=~ \frac{2(z^2 - z)}{3z - 1}\,\).

Using the relation \(z = e^{\i\theta}\) for the locus of points on the boundary of the unit circle we get

\[ \begin{aligned} h\lambda ~=~ \frac{2(e^{\i 2\theta} - e^{i\theta})}{3e^{\i\theta} - 1}\,. \end{aligned} \]

The following derivations illustrate in detail the steps required for finding the \(x\), \(y\) coordinates of the points on the boundary locus of the region of absolute stability. We can use \(\,e^{\i n\theta} = \cos n\theta + \i\sin n\theta\,\) and get:

\[ \begin{aligned} h\lambda ~=~ \frac{2\,(\cos2\theta + \i\sin2\theta - \cos\theta - \i\sin\theta)}{3(\cos\theta + \i\sin\theta) - 1} \end{aligned} \]

Since \(\,h\lambda\) expression is complex, we can now separate the real and imaginary components as

\[\begin{split} \begin{aligned} h\lambda ~=~ x(\theta) + \i\,y(\theta) \qquad \text{where}&\qquad x(\theta) ~=~\frac{2(\cos\theta - 1)^2}{3\cos\theta - 5} \\[1ex] \text{and}&\qquad y(\theta) ~=~ \frac{2\sin\theta\,(\cos\theta - 2)}{3\cos\theta - 5}\,. \end{aligned} \end{split}\]

It can easily be seen that \(x(-\theta) = x(\theta)\) and \(y(-\theta) = -y(\theta)\), so that the locus is symmetric about the x-axis. Evaluating \(x\) and \(y\) at intervals of \(\theta\) of \(30^{\circ}\), we obtain the following table:

\(\theta\)	\(0^{\circ}\)	\(30^{\circ}\)	\(60^{\circ}\)	\(90^{\circ}\)	\(120^{\circ}\)	\(150^{\circ}\)	\(180^{\circ}\)
\(x\)	0.000	-0.015	-0.143	-0.400	-0.692	-0.917	-1.000
\(y\)	0.000	0.472	0.742	0.800	0.666	0.377	0.000

These points are plotted in the complex \(h\lambda\) plane, and the resulting region \(\Re\), \(\Re \subset \C\), of absolute stability is shown in the following figure:

Fig. 7.1 Stability region of AB 2-step method.

The region of absolute stability can be easily found by using a mathematical package. A worksheet solution for Example 3 above is attached which illustrates the computational steps involved when a Matlab program for finding the region of absolute stability is used.

The boundaries of the regions of absolute stability for various order Adams methods are shown

Fig. 7.2 Stability region of AB \(k\)-step method.

Fig. 7.3 Stability region of AM \(k\)-step method.

We learn from these diagrams that the region of absolute stability becomes smaller as the order of the method increases; for the formula of the same order the Adams-Moulton formula has a significantly larger region of absolute stability than the Adams-Bashforth formula. The size of these regions is usually quite acceptable for a practical point of view. For example, the real values of \(\,h\lambda\,\) in the region of absolute stability for the fourth-order Adams-Moulton formula are given by \(~-3 < h\lambda < 0\,\). This is not a serious restriction on \(h\) in most cases.

The Adams-family of formulae is very convenient for creating a variable order program, and their stability regions are quite acceptable. There will be difficulties with problems for which \(\lambda\) is negative and large in magnitude; these problems are best treated by other methods.

There are special methods for which the region of absolute stability consists of all complex values of \(\,h\lambda\,\) with \(\,Real(h\lambda) < 0\,\). These methods are called A-stable, and with them there is no restriction on \(h\) in order to have stability of the type we have been considering.

Definition 7.4 (\(A\)-stability)

The property that \(\displaystyle\,\lim_{j \to \infty} y_j = 0\;\) (\(y_j\,\) approximated solution) for all \(\,h>0\,\) and all complex \(\lambda\) with \(\,Real(\lambda)<0\,\) is called A-stability.

For certain applications, stiff differential equations, the A-stability property of the numerical method is an important consideration.

Example 7.5

Show that the backward Euler method (or first order Adams-Moulton formula)

\[ \begin{aligned} y_{j+1} ~=~ y_j + h\,f_{j+1} \end{aligned} \]

is A-stable.

Solution

Applying this formula to the test equation \(y' = \lambda y\) we find

\[\begin{split} \begin{aligned} y_{j+1} & =  y_j + hy_{j+1} \\ (1-h\lambda) y_{j+1} & = y_j \\ (1-h\lambda) y_{j+1} - y_j & =0 \end{aligned} \end{split}\]

This has a solution of the form \(y_j = A_1z_1^j\).

\[ \begin{aligned} && (1 - h\lambda)z - 1 ~&=~ 0 & \therefore\quad z ~&=~ \frac{1}{1 - h\lambda}\,. && \end{aligned} \]

Substituting gives:

\[ \begin{aligned} y_j ~&=~ A_1 \left( \frac{1}{1 - h\lambda} \right)^j,\quad \text{where $A_1$ is constant;} \end{aligned} \]

then \(y_j \to 0\) as \(x_j \to \infty\) if and only if

\[ \begin{aligned} \frac{1}{|1 - h\lambda|} \,<\, 1\,. \end{aligned} \]

This will be true for all \(h\lambda\) with \(Real(\lambda) < 0\), and thus the backward Euler formula is an A-stable method.

It would have been useful to have A-stable multistep methods of order greater than 2, but it has been shown that, Dahlquist (1963), there are no such methods.