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

u_t + c u_x = 0

(12.3.1)

be posed for $x \in [0,1]$ and subjected to periodic end conditions. If we use the central-difference matrix $\mathbf{D}_x$ defined in (12.1.9) to discretize the space derivative, we get the ODE system

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

(12.3.2)

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 12.3.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.3)

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.

Source

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):

Source

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.

Source

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.

Source

[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):

Source

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')

Source

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.

Source

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):

Source

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");

Source

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})$ .

12.3.1Advection–diffusion equations¶

The traffic flow equation (12.1.5) 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.4)

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.

Source

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.

Source

[x, Dx, Dxx] = diffper(40, [0, 1]);
clf
for ep = [0.001 0.01 0.05]
  lambda = eig(-Dx + ep*Dxx);
  str = sprintf("\\epsilon = %.3f", ep);
  scatter(real(lambda), imag(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.

Source

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

If we use a time stepping method with fixed step size τ on this system, then for any given ε and $m$ , there is a maximum value $\hat{\tau}$ of τ such that $\hat{\tau} \lambda$ fits inside the stability region of the method for all eigenvalues λ. Then the stability restriction for the fully discrete method is $\tau \le \hat{\tau}$ . (See Exercise 12.3.2.)

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, Example 12.2.7 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.5)

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.6)

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);

Source

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(λ))

-4.931967035822774

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);

Source

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}")

rightmost extent of eigenvalues: -4.93

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¶

Exercise 12.3.1

✍ (Similar to Exercise 11.3.7.) Let $\mathbf{D}_{x}$ be $m\times m$ and given by (12.2.1) for periodic end conditions. For any integer $k \in \{0,\ldots,m-1\}$ , define $\omega = \exp(2ik\pi/m)$ , and let $\mathbf{v}$ be the vector whose components are $v_j = \omega^j$ for $j=0,\ldots,m-1$ .

(a) Show that $\omega^m = 1$ .

(b) Let $\mathbf{v}' = \mathbf{D}_x \mathbf{v}$ . Show that for $j=1,\ldots,m-2$ ,

v_j' = \frac{1}{2h} \omega^{j} \left( \omega - \omega^{-1} \right).

(12.3.7)

(c) Show that the result of part (b) holds for $j=0$ and $j=m-1$ as well.

(d) Explain why the above results prove that $\mathbf{v}$ is an eigenvector of $\mathbf{D}_x$ with associated eigenvalue

\lambda = i\, m \sin \left( \frac{2k\pi}{m} \right).

(12.3.8)

Exercise 12.3.2

⌨ In Example 12.3.2 we saw that the eigenvalues of the semidiscretization of (12.3.1) for periodic end conditions lie in the left half of the complex plane. Suppose we want to apply the Euler time stepping formula. For a given eigenvalue λ, there is a value of τ such that $\zeta=\tau \lambda$ lies on the boudnary of the stability region.

(a) ✍ Show that if $\lambda = x + iy$ and $|\zeta + 1|^2 = 1$ , then

\tau = -\frac{2x}{x^2 + y^2}.

(12.3.9)

(b) ⌨ Use $c=1,$ $\epsilon=0.001,$ and $m=100$ to find the eigenvalues. For each eigenvalue, use (12.3.9) to find the value of τ that scales it to the stability region. Then find the minimum value of τ over all the eigenvalues. (This is the maximum stable time step $\hat{\tau}$ for Euler.)

Preface

Upwinding and stability

Preface

The wave equation