Absolute stability - Fundamentals of Numerical Computation

In The method of lines we applied several different time stepping methods to a linear, constant coefficient problem in the form

\mathbf{u}'(t)=\mathbf{A}\mathbf{u}(t).

(11.3.1)

All of these methods are zero-stable in the sense of Zero-stability of multistep methods, in the limit as the time step size $\tau \to 0$ .^[1] Yet for some experiments with fixed τ, as in Demo 11.2.2, we have observed exponential growth in the different limit $n\to \infty$ .

Observe that if $\mathbf{A}$ has the eigenvalue decomposition $\mathbf{A}=\mathbf{V}\mathbf{D}\mathbf{V}^{-1}$ , then

\begin{align*} \mathbf{u}'&=(\mathbf{V}\mathbf{D}\mathbf{V}^{-1})\mathbf{u},\\ (\mathbf{V}^{-1} \mathbf{u}') &= \mathbf{D} (\mathbf{V}^{-1} \mathbf{u}), \\ \mathbf{y}' &= \mathbf{D} \mathbf{y}, \end{align*}

(11.3.2)

where $\mathbf{y}(t)=\mathbf{V}^{-1}\mathbf{u}(t)$ . Because $\mathbf{D}$ is diagonal, the dynamics of the components of $\mathbf{y}$ are completely decoupled: each row is a self-contained equation of the form $y_j'=\lambda_j y_j$ , where $\lambda_j$ is an eigenvalue of $\mathbf{A}$ .

The diagonalization argument suggests that we can look at the scalar problems

y' = \lambda y, \quad y(0)=1,

(11.3.3)

arising from the eigenvalues. These eigenvalues may not be real numbers, so in this section, $i$ stands for the imaginary unit, not an integer index. If we write λ in real and imaginary parts as $\lambda=\alpha + i\beta$ , then by Euler’s identity, the exact solution of (11.3.3) has magnitude

\bigl |e^{(\alpha+i\beta)t} \bigr| = \bigl |e^{\alpha t} \bigr| \cdot \bigl |e^{i \beta t} \bigr| = e^{\alpha t}.

(11.3.4)

We now consider the counterpart of this observation for the solution produced by a numerical IVP solver.

The fact that absolute stability depends only on the product $\zeta = \tau\lambda$ , and not independently on the individual factors, is a result of how the IVP solvers are defined, as we will see below. Since λ has units of inverse time according to (11.3.3), ζ is dimensionless.

11.3.1Stability regions¶

Each numerical IVP solver has its own collection of ζ values for which it is absolutely stable.

Example 11.3.1

Consider an Euler discretization of $y'=\lambda y$ :

y_{k+1} = y_k + \tau( \lambda y_k) = (1+ \zeta ) y_k.

(11.3.5)

Given that $y_0=1$ by (11.3.3), we easily deduce that $y_k = (1+\zeta)^k$ for all $k$ , and therefore

|y_k| = |1+\zeta|^k.

(11.3.6)

Hence $|y_k|$ remains bounded above as $k\to \infty$ if and only if $|1+\zeta| \le 1$ . Because ζ is a complex number, it’s easiest to interpret this condition geometrically:

|\zeta + 1 | = |\zeta - (-1) | \le 1.

(11.3.7)

That is, the distance in the plane from ζ to the point -1 is less than or equal to 1. This description defines a closed disk of radius 1 centered at $(-1,0)$ .

Stability regions for the most common IVP integrators are given in Figure 11.3.1 and Figure 11.3.2. Note that those for the implicit Adams-Moulton methods are larger than those for the explicit Adams-Bashforth methods of the same order. For the implicit backward differentiation methods, the exteriors of the curves provide large regions of stability, but significant portions of the imaginary axis may be excluded. Finally, while the single-step Runge-Kutta methods have smaller regions of stability, those of orders 3 and 4 do include significant portions of the imaginary axis.

Stability regions for Adams–Bashforth methods of order 1–4 (left) and Adams–Moulton methods of order 2–5 (right). The plots are in the complex \zeta-plane. — Figure 11.3.1:Stability regions for Adams–Bashforth methods of order 1–4 (left) and Adams–Moulton methods of order 2–5 (right). The plots are in the complex ζ-plane.

Stability regions for backward differentiation methods of order 1–4 (left, exteriors of curves) and Runge–Kutta methods of order 1–4 (right). The plots are in the complex \zeta-plane. — Figure 11.3.2:Stability regions for backward differentiation methods of order 1–4 (left, exteriors of curves) and Runge–Kutta methods of order 1–4 (right). The plots are in the complex ζ-plane.

For any particular method and value of λ in (11.3.3), we can use the stability region to deduce which, if any, values of the time step τ will give bounded solutions. Both the magnitude and the argument (angle) of λ play a role in determining such constraints.

Example 11.3.4

Suppose $\lambda=-4$ and Euler’s method is applied. Since the time step is always positive, $\zeta=-4\tau$ is always on the negative real axis. The only part of that line that lies within the stability region of Euler as derived in Example 11.3.1 is the real interval $[-2,0]$ . Hence we require $\zeta\ge -2$ , or $\tau \le 1/2$ . By contrast, the stability region of backward Euler includes the entire negative real axis, so absolute stability is unconditional, i.e., assured regardless of τ.

Now suppose instead that $\lambda=i$ , so that $\zeta=i\tau$ . Clearly ζ is always on the positive imaginary axis. But no part of this axis, aside from the origin, lies in the stability region of Euler’s method, so it is unconditionally unstable in this circumstance. The conclusion for backward Euler is the opposite; any value of τ will do, because the entire imaginary axis is within the stability region.

Example 11.3.4 does not contradict our earlier statements about the zero stability and convergence of Euler’s method in general, even for the case $\lambda=i$ . But those statements are based on the limit $\tau\to 0$ for $t$ in a finite interval $[a,b]$ . Both this limit and the limit $t\to \infty$ imply the number of steps $n$ goes to infinity, but the limits behave differently.

The fact that implicit methods have larger stability regions than their explicit counterparts is the primary justification for using them. While they have larger work requirements per step, they sometimes can take steps that are orders of magnitude larger than explicit methods and still remain stable.

When adaptive time stepping methods are used, as in most software for IVPs, the automatically determined time step is chosen to satisfy absolute stability requirements (otherwise errors grow exponentially). This phenomenon was manifested in Demo 11.2.5: in the explicit IVP method rk23, error control forced tiny step sizes compared to those used by Rodas4P, which is based on implicit methods.

11.3.2Heat equation¶

Now we return to the semidiscretization (11.2.5) of the heat equation, which was solved by Euler in Demo 11.2.2 and backward Euler in Demo 11.2.4.

Example 11.3.5 (Stability regions and the heat equation)

Julia

MATLAB

Python

Example 11.3.5

Euler and Backward Euler time-stepping methods were used to solve $\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}$ .

m = 40
_, _, Dₓₓ = FNC.diffper(m, [0, 1]);

The eigenvalues of this matrix are real and negative:

using Plots
λ = eigvals(Dₓₓ)
scatter(real(λ), imag(λ);
    title="Eigenvalues",
    frame=:zerolines,  aspect_ratio=1,
    xaxis=("Re λ"),  yaxis=("Im λ", (-1000, 1000)))

The Euler method is absolutely stable in the region $|\zeta+1| \le 1$ in the complex plane:

phi = 2π * (0:360) / 360
z = @. cis(phi) - 1;    # unit circle shifted to the left by 1

plot(Shape(real(z), imag(z));
    color=RGB(.8, .8, 1),
    xaxis=("Re ζ"),  yaxis=("Im ζ"),
    aspect_ratio=1,  frame=:zerolines,
    title="Stability region")

In order to get inside this region, we have to find τ such that $\lambda \tau > -2$ for all eigenvalues λ. This is an upper bound on τ.

λ_min = minimum(λ)
@show max_τ = -2 / λ_min;

Here we plot the resulting values of $\zeta=\lambda \tau$ .

ζ = λ * max_τ
scatter!(real(ζ), imag(ζ), title="Stability region and ζ values")

In backward Euler, the region is $|\zeta-1|\ge 1$ . Because they are all on the negative real axis, all of the ζ values will fit no matter what τ is chosen.

plot(Shape([-6, 6, 6, -6], [-6, -6, 6, 6]), color=RGB(.8, .8, 1))
z = @. cis(phi) + 1;   # unit circle shifted right by 1
plot!(Shape(real(z), imag(z)), color=:white)

scatter!(real(ζ), imag(ζ);
    xaxis=([-4, 2], "Re ζ"),  yaxis=([-3, 3], "Im ζ"),
    aspect_ratio=1,  frame=:zerolines,
    title="Stability region and ζ values")

Example 11.3.5

Euler and Backward Euler time-stepping methods were used to solve $\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}$ .

m = 40;  
[x, Dx, Dxx] = diffper(m, [0, 1]);

The eigenvalues of this matrix are real and negative:

lambda = eig(Dxx);
clf
plot(real(lambda), imag(lambda), 'o')
axis equal, grid on
xlabel('Re \lambda'),  ylabel('Im \lambda')
title('Eigenvalues')

The Euler method is absolutely stable in the region $|\zeta+1| \le 1$ in the complex plane:

phi = linspace(0, 2*pi, 361);
z = exp(1i*phi) - 1;   % unit circle shifted to the left by 1
fill(real(z), imag(z), [.8, .8, 1])
axis equal,  grid on
xlabel('Re \lambda'),  ylabel('Im \lambda')
title('Stability region')

In order to get inside this region, we have to find τ such that $\lambda \tau > -2$ for all eigenvalues λ. This is an upper bound on τ.

lambda_min = min(lambda)
max_tau = -2 / lambda_min

Here we plot the resulting values of $\zeta=\lambda \tau$ .

zeta = lambda * max_tau;
hold on
plot(real(zeta), imag(zeta), 'o')
title('Stability region and \zeta values')

In backward Euler, the region is $|\zeta-1|\ge 1$ . Because they are all on the negative real axis, all of the ζ values will fit no matter what τ is chosen.

clf,  fill([-6, 6, 6, -6],[-6, -6, 6, 6],[.8, .8, 1])
hold on
z = exp(1i*phi) + 1;   % unit circle shifted right by 1
fill(real(z), imag(z), 'w')
plot(real(zeta), imag(zeta), 'o')
axis equal,  axis([-4, 2, -3, 3]), grid on
xlabel('Re \lambda'),  ylabel('Im \lambda')
title('Stability region and \zeta values')

Example 11.3.5

Euler and Backward Euler time-stepping methods were used to solve $\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}$ .

m = 40
Dxx = FNC.diffper(m, [0, 1])[2]

The eigenvalues of this matrix are real and negative:

from scipy.linalg import eigvals
lamb = eigvals(Dxx)
plot(real(lamb), imag(lamb), "o")
xlabel("Re $\\lambda$")
ylabel("Im $\\lambda$")
title("Eigenvalues");

The Euler method is absolutely stable in the region $|\zeta+1| \le 1$ in the complex plane:

phi = 2 * pi * arange(361) / 360
z = exp(1j * phi) - 1  # unit circle shifted to the left by 1

fill(real(z), imag(z), color=(0.8, 0.8, 1))
xlabel("Re $\\lambda$")
ylabel("Im $\\lambda$")
axis("equal")
title("Stability region");

In order to get inside this region, we have to find τ such that $\lambda \tau > -2$ for all eigenvalues λ. This is an upper bound on τ.

lambda_min = min(real(lamb))
max_tau = -2 / lambda_min
print(f"predicted max time step is {max_tau:.3e}")

Here we plot the resulting values of $\zeta=\lambda \tau$ .

zeta = lamb * max_tau
fill(real(z), imag(z), color=(0.8, 0.8, 1))
plot(real(zeta), imag(zeta), "o")
xlabel("Re $\\lambda$")
ylabel("Im $\\lambda$")
axis("equal")
title("Stability region and $\zeta$ values");

In backward Euler, the region is $|\zeta-1|\ge 1$ . Because they are all on the negative real axis, all of the ζ values will fit no matter what τ is chosen.

fill([-6, 6, 6, -6], [-6, -6, 6, 6], color=(0.8, 0.8, 1))
z = exp(1j * phi) + 1
# unit circle shifted right by 1
fill(real(z), imag(z), color="w")

plot(real(zeta), imag(zeta), "o")
axis([-4, 2, -3, 3])
axis("equal")
xlabel("Re $\\lambda$")
ylabel("Im $\\lambda$")
title("Stability region and $\\zeta$ values");

The matrix $\mathbf{D}_{xx}$ occurring in (11.2.5) for semidiscretization of the periodic heat equation has eigenvalues that can be found explicitly. Assuming that $x\in[0,1)$ (with periodic boundary conditions), for which $h=1/m$ , then the eigenvalues are

\lambda_j = -4m^2 \sin^2 \left( \frac{j\pi}{m} \right), \qquad j = 0,\ldots,m-1.

(11.3.12)

This result agrees with the observation in Demo 11.3.5 that the eigenvalues are real and negative. Furthermore, they lie within the interval $[-4m^2,0]$ . In Euler time integration, this implies that $-4\tau m^2\ge -2$ , or $\tau\ge 1/(2m^2)=O(m^{-2})$ . For backward Euler, there is no time step restriction, and we say that backward Euler is unconditionally stable for this problem.

In summary, three things happen as $h\to 0$ :

The spatial discretization becomes more accurate like $O(h^2)$ .
The size of the matrix increases like $O(h^{-1})$ .
If we use an explicit time stepping method, then absolute stability requires $O(h^{-2})$ steps.

The last restriction becomes rather burdensome as $h\to 0$ , i.e., as we improve the spatial discretization, which is why implicit methods are preferred for diffusion. While any convergent IVP solver will get the right solution as $\tau\to 0$ , the results are exponentially large nonsense until τ is small enough to satisfy absolute stability.

11.3.3Exercises¶

✍ Use an eigenvalue decomposition to write the system
$\mathbf{u}'(t) = \begin{bmatrix} 0 & 4 \\ -4 & 0 \end{bmatrix} \mathbf{u}(t)$
(11.3.13)
as an equivalent diagonal system.
✍ For each system, state whether its solutions are bounded as $t\to \infty$ .
(a) $\mathbf{u}'(t) = \displaystyle \begin{bmatrix} 1 & 3 \\ 3 & 1 \end{bmatrix} \mathbf{u}(t)$
(b) $\mathbf{u}'(t) = \displaystyle \begin{bmatrix} -1 & 3 \\ -3 & -1 \end{bmatrix} \mathbf{u}(t)$
(c) $\mathbf{u}'(t) = \displaystyle \begin{bmatrix} 0 & 4 \\ -4 & 0 \end{bmatrix} \mathbf{u}(t)$
✍ Using Figure 11.3.1 and Figure 11.3.2, estimate the time step restriction (if any) for the system
$\mathbf{u}'(t) = \begin{bmatrix} -4 & 0 & 0 \\ 0 & -2 & 0 \\ 0 & 0 & -0.5 \end{bmatrix} \mathbf{u}(t)$
(11.3.14)
for the following IVP methods:
(a) RK4 $\qquad$ (b) AM4 $\qquad$ (c) AB2
✍ Using Figure 11.3.1 and Figure 11.3.2, find the time step restriction (if any) for the system
$\mathbf{u}'(t) = \begin{bmatrix} -1 & 0 & 0 \\ 0 & 0 & 4 \\ 0 & -4 & 0 \end{bmatrix} \mathbf{u}(t)$
(11.3.15)
for the following IVP methods:
(a) RK4 $\qquad$ (b) AM4 $\qquad$ (c) AB3
✍ Of the following methods, which would be unsuitable for a problem having eigenvalues on the imaginary axis? Justify your answer(s).
(a) AM2 $\qquad$ (b) AB2 $\qquad$ (c) RK2 $\qquad$ (d) RK3
✍ Of the following methods, which would have a time step restriction for a problem with real, negative eigenvalues? Justify your answer(s).
(a) AM2 $\qquad$ (b) AM4 $\qquad$ (c) BD4 $\qquad$ (d) RK4

✍ Let $\mathbf{D}_{xx}$ be $m\times m$ and given by (11.2.4). For any integer $k \in \{0,\ldots,m-1\}$ , define $\omega = \exp(2i k\pi/m)$ and $\mathbf{v} = \bigl[ 1,\; \omega,\; \omega^2,\; \ldots,\; \omega^{m-1} \bigr].$ Show that $\mathbf{v}$ is an eigenvector of $\mathbf{D}_{xx}$ , with eigenvalue
$\lambda = -4m^2 \sin^2 \left( \frac{k\pi}{m} \right).$
(11.3.16)
(This establishes that the eigenvalues all lie within the real interval $[-4m^2,0]$ .)
✍ (a) Derive an algebraic inequality equivalent to absolute stability for the AM2 (trapezoid) formula.
✍ (b) Argue that the inequality in part (a) is equivalent to the restriction $\operatorname{Re} \zeta\le 0$ . (Hint: Complex magnitude is equivalent to distance in the plane.)

Footnotes¶

In Chapter 6 we used $h$ rather than τ to denote the time step size, but now we reserve $h$ for spacing in the $x$ direction.
↩