Absolute stability - Fundamentals of Numerical Computation

The CFL criterion gives a necessary condition for convergence. It suggests, but cannot confirm, that a step size of $O(h)$ may be adequate in the advection equation. More details emerge when we adopt the semidiscretization point of view.

Let the advection equation (12.2.1) over $[0,1]$ be subjected to periodic end conditions. Suppose we use the central-difference matrix $\mathbf{D}_x$ defined in (12.1.7) to discretize the space derivative, leaving us with

\mathbf{u}' = -c \mathbf{D}_x \mathbf{u}.

(12.3.1)

To apply an IVP solver, we need to compare the stability region of the solver with the eigenvalues of $-c \mathbf{D}_x$ , as in Absolute stability. You can verify (see Exercise 1) that for $m$ points in $[0,1)$ , these are

\lambda_k = - i\, c m \sin \left( \frac{2\pi k}{m} \right), \qquad k = 0,\ldots,m-1.

(12.3.2)

Two things stand out about these eigenvalues: they are purely imaginary, which is consistent with conservation of magnitude, and they extend no farther than $O(m)=O(h^{-1})$ away from the origin. These characteristics suggest how to analyze the use of different time-stepping methods by referring to stability regions.

Example 12.3.1 (Eigenvalues for advection)

Julia

MATLAB

Python

Example 12.3.1

For $c=1$ we get purely imaginary eigenvalues.

using Plots
x, Dₓ = FNC.diffper(40, [0, 1])
λ = eigvals(Dₓ);
scatter(real(λ), imag(λ);
    aspect_ratio = 1,  frame=:zerolines,
    xlabel="Re λ",  ylabel="Im λ", 
    title="Eigenvalues for pure advection",  legend=:none)

Let’s choose a time step of $\tau=0.1$ and compare to the stability regions of the Euler and backward Euler time steppers (shown as shaded regions):

zc = @. cispi(2 * (0:360) / 360);     # points on |z|=1
z = zc .- 1;                          # shift left by 1
plot(Shape(real(z), imag(z)), color=RGB(.8, .8, 1))
ζ = 0.1 * λ
scatter!(real(ζ), imag(ζ);
    aspect_ratio=1,  frame=:zerolines,
    xaxis=("Re ζ", [-5, 5]),  yaxis=("Im ζ", [-5, 5]),
    title="Euler for advection")

In the Euler case it’s clear that no real value of $\tau>0$ is going to make ζ values fit within the stability region. Any method whose stability region includes none of the imaginary axis is an unsuitable choice for advection.

z = zc .+ 1;                        # shift circle right by 1
plot(Shape([-6, 6, 6, -6], [-6, -6, 6, 6]), color=RGB(.8, .8, 1))
plot!(Shape(real(z), imag(z)), color=:white)
scatter!(real(ζ), imag(ζ);
    aspect_ratio=1,  frame=:zerolines,
    xaxis=("Re ζ", [-5, 5]),  yaxis=("Im ζ", [-5, 5]),
    title="Backward Euler for advection")

The A-stable backward Euler time stepping tells the exact opposite story: it will be absolutely stable for any choice of the time step τ.

Example 12.3.1

For $c=1$ we get purely imaginary eigenvalues.

[x, Dx] = diffper(40, [0, 1]);
lambda = eig(Dx);
clf
scatter(real(lambda), imag(lambda))
axis equal,  grid on
title('Eigenvalues for pure advection')

Let’s choose a time step of $\tau=0.1$ and compare to the stability regions of the Euler and backward Euler time steppers (shown as shaded regions):

zc = exp( 1i * linspace(0,2*pi,361)' );    % points on |z|=1
z = zc - 1;    % shift circle left by 1
clf,  fill(real(z), imag(z), [.8, .8, 1])
hold on,  scatter(real(0.1*lambda), imag(0.1*lambda))
axis equal,  axis square
axis([-5, 5, -5, 5]),  grid on
title('Euler')

clf,  fill([-6, 6, 6, -6],[-6, -6, 6, 6], [.8, .8, 1])
z = zc + 1;   % shift circle right by 1
hold on,  scatter(real(0.1*lambda), imag(0.1*lambda))
fill(real(z), imag(z), 'w')
axis equal,  axis square
axis([-5 5 -5 5]),  grid on
title('backward Euler')

The A-stable backward Euler time stepping tells the exact opposite story: it will be absolutely stable for any choice of the time step τ.

Example 12.3.1

For $c=1$ we get purely imaginary eigenvalues.

from scipy.linalg import eigvals
x, Dx, Dxx = FNC.diffper(40, [0, 1])
lamb = eigvals(Dx)

plot(real(lamb), imag(lamb), "o")
axis([-40, 40, -40, 40])
axis("equal")
title("Eigenvalues for pure advection");

Let’s choose a time step of $\tau=0.1$ and compare to the stability regions of the Euler and backward Euler time steppers (shown as shaded regions):

zc = exp(2j * pi * arange(361) / 360)
# points on |z|=1

z = zc - 1    # shift circle left by 1
fill(real(z), imag(z), color=(0.8, 0.8, 1))
plot(real(0.1 * lamb), imag(0.1 * lamb), "o")
axis([-5, 5, -5, 5]),  axis("equal")
title("Euler");

z = zc + 1    # shift circle right by 1
fill([-6, 6, 6, -6], [-6, -6, 6, 6], color=(0.8, 0.8, 1))
fill(real(z), imag(z), color="w")
plot(real(0.1 * lamb), imag(0.1 * lamb), "o")
axis([-5, 5, -5, 5])
axis("equal")
title("Backward Euler");

The A-stable backward Euler time stepping tells the exact opposite story: it will be absolutely stable for any choice of the time step τ.

Many PDEs that conserve quantities will have imaginary eigenvalues, causing Euler and some other IVP methods to fail regardless of step size. Diffusion problems, in which the eigenvalues are negative and real, are compatible with a wider range of integrators, though possibly with onerous step size requirements due to stiffness.

The location of eigenvalues near $\pm ic/h$ also confirms what the CFL condition was suggesting. In order to use RK4, for example, whose stability region intersects the imaginary axis at around $\pm 2.8i$ , the time step stability restriction is $\tau c/h \le 2.8$ , or $\tau=O(h)$ . This is much more favorable than for diffusion, whose eigenvalues were as large as $O(h^{-2})$ , and it makes explicit IVP methods much more attractive for advection problems than for diffusion.

12.3.1Advection–diffusion equations¶

The traffic flow equation (12.1.4) combines a nonlinear advection with a diffusion term. The simplest linear problem with the same feature is the advection–diffusion equation

u_t+c u_x=\epsilon u_{xx}.

(12.3.3)

The parameter ε controls the relative strength between the two mechanisms, and the eigenvalues accordingly vary between the purely imaginary ones of advection and the negative real ones of diffusion.

Example 12.3.2 (Eigenvalues for advection–diffusion)

Julia

MATLAB

Python

Example 12.3.2

The eigenvalues of advection-diffusion are near-imaginary for $\epsilon\approx 0$ and get closer to the negative real axis as ε increases.

plt = plot(
    legend=:topleft,
    aspect_ratio=1,
    xlabel="Re ζ",  ylabel="Im ζ",
    title="Eigenvalues for advection-diffusion")
x, Dₓ, Dₓₓ = FNC.diffper(40, [0, 1]);
for ϵ in [0.001, 0.01, 0.05]
    λ = eigvals(-Dₓ + ϵ*Dₓₓ)
    scatter!(real(λ), imag(λ), m=:o, label="\\epsilon = $ϵ")
end
plt

Example 12.3.2

The eigenvalues of advection-diffusion are near-imaginary for $\epsilon\approx 0$ and get closer to the negative real axis as ε increases.

[x, Dx, Dxx] = diffper(40, [0, 1]);
tau = 0.1;
clf
for ep = [0.001 0.01 0.05]
  lambda = eig(-Dx + ep*Dxx);
  str = sprintf("\\epsilon = %.3f", ep);
  scatter(real(tau*lambda), imag(tau*lambda), displayname=str)
  hold on
end
axis equal,  grid on
legend(location='northwest')
title('Eigenvalues for advection-diffusion')

Example 12.3.2

The eigenvalues of advection-diffusion are near-imaginary for $\epsilon\approx 0$ and get closer to the negative real axis as ε increases.

x, Dx, Dxx = FNC.diffper(40, [0, 1])
tau = 0.1
for ep in [0.001, 0.01, 0.05]:
    lamb = eigvals(-Dx + ep * Dxx)
    plot(real(tau * lamb), imag(tau * lamb), "o", label=f"epsilon={ep:.1g}")
axis("equal")
legend()
title("Eigenvalues for advection-diffusion")

In a nonlinear problem, the eigenvalues come from the linearization about an exact solution, as in Stiffness.

12.3.2Boundary effects¶

Boundary conditions can have a dramatic effect on the eigenvalues of the semidiscretization. For instance, Demo 12.2.6 solves linear advection $u_t=u_x$ on $[0,1]$ with the homogeneous inflow condition $u(0,t)=0$ . Exclusion of the boundary node from the semidiscretization $\mathbf{u}$ to get the interior vector $\mathbf{v}$ is equivalent to

\mathbf{v} = \mathbf{E} \mathbf{u}, \quad \mathbf{u} = \begin{bmatrix} \mathbf{v} \\ 0 \end{bmatrix} = \mathbf{E}^T \mathbf{v},

(12.3.4)

where $\mathbf{E}$ is the $(m+1)\times (m+1)$ identity with the last row deleted. The ODE on the interior nodes is

\frac{d\mathbf{v}}{dt} = \mathbf{E} \left( \mathbf{D}_x \mathbf{u} \right) = \mathbf{E} \mathbf{D}_x \mathbf{E}^T \mathbf{v}.

(12.3.5)

As a result, we conclude that $\mathbf{A} = \mathbf{E} \mathbf{D}_x \mathbf{E}^T$ is the appropriate matrix for determining the eigenvalues of the semidiscretization. More simply, we can simply delete the last row and last column from $\mathbf{D}_x$ .

Example 12.3.3 (Eigenvalues for an inflow boundary)

Julia

MATLAB

Python

Example 12.3.3

Deleting the last row and column places all the eigenvalues of the discretization into the left half of the complex plane.

x, Dₓ, _ = FNC.diffcheb(40, [0, 1])
A = Dₓ[1:end-1, 1:end-1];     # delete last row and column
λ = eigvals(A);

scatter(real(λ), imag(λ);
    m=3,  aspect_ratio=1,
    legend=:none,  frame=:zerolines,
    xaxis=([-300, 100], "Re λ"),  yaxis=("Im λ"),
    title="Eigenvalues of advection with zero inflow")

Note that the rightmost eigenvalues have real part at most

maximum(real(λ))

Consequently all solutions decay exponentially to zero as $t\to\infty$ . This matches our observation of the solution: eventually, everything flows out of the domain.

Example 12.3.3

Deleting the last row and column places all the eigenvalues of the discretization into the left half of the complex plane.

[x, Dx, Dxx] = diffcheb(40, [0, 1]);
A = -Dx(2:end, 2:end);    % leave out first row and column
lambda = eig(A);

clf
scatter(real(lambda), imag(lambda))
axis equal,  grid on 
title('Eigenvalues of advection with zero inflow')

Note that the rightmost eigenvalues have real part at most

max(real(lambda))

Consequently all solutions decay exponentially to zero as $t\to\infty$ . This matches our observation of the solution: eventually, everything flows out of the domain.

Example 12.3.3

Deleting the last row and column places all the eigenvalues of the discretization into the left half of the complex plane.

from scipy.linalg import eigvals
x, Dx, Dxx = FNC.diffcheb(40, [0, 1])
A = -Dx[1:, 1:]  # leave out first row and column

lamb = eigvals(A)
plot(real(lamb), imag(lamb), "o")
xlim(-300, 100),  axis("equal"),  grid(True)
title("Eigenvalues of advection with zero inflow");

Note that the rightmost eigenvalues have real part at most

print(f"rightmost extent of eigenvalues: {max(real(lamb)):.3g}")

Consequently all solutions decay exponentially to zero as $t\to\infty$ . This matches our observation of the solution: eventually, everything flows out of the domain.

12.3.3Exercises¶

✍ Let $\mathbf{D}_{x}$ be $m\times m$ and given by (12.2.2). For any integer $k \in \{0,\ldots,m-1\}$ , define $\omega = \exp(2ik\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}_{x}$ , with eigenvalue
$\lambda = i\, m \sin \frac{2k\pi}{m}.$
(12.3.6)
(See also Exercise 11.3.7.)
⌨ Refer to the semidiscretization of the advection–diffusion equation on $x\in[0,1]$ with $c=1$ , $\epsilon=0.01$ , and periodic end conditions in Demo 12.3.2. For $m=100$ , find the upper bound on τ that gives absolute stability for Euler time stepping. (Hint: The stability region of Euler is the set of complex values whose distance from -1 is less than or equal to one. The effect of τ is to uniformly scale that distance for the eigenvalues.)
⌨ Refer to the semidiscretization in Demo 12.3.3. Find the upper bound on τ that gives absolute stability for Euler time stepping. (See the hint in the previous exercise.)
⌨ Modify Demo 12.3.3 so that it produces the eigenvalues of the problem $u_t+u_x=0$ with an outflow condition $u(1,t)=0$ . What is the behavior of solutions as $t\to\infty$ ?

Preface

Upwinding and stability

Preface

The wave equation