8) Hyperbolic Equations#
Last time:
Measuring errors
Stability
Consistency
Second order derivatives
Matrix representations and properties
Second order derivative with Dirichelet boundary conditions
Discrete Green’s functions
Interpolation by Vandermonde matrices
Today:
Introduction to Hyperbolic PDEs
1.1. The wave equation
1.2 A more general hyperbolic equationSystem of Hyperbolic Equations
Hyperbolic Equations with Variable Coefficients
3.1 Hyperbolic Systems with Variable CoefficientsBoundary Conditions
FD schemes for hyperbolic equations
5.1 FD schemes for the 1D wave equation
5.2 Example using Leapfrog
5.3 Example using Lax-FriedrichsStability of hyperbolic PDEs
6.1 The Courant-Friedrichs-Lewy ConditionConditional and unconditional stability
1. Introduction to Hyperbolic PDEs#
We begin with an overview of Hyperbolic PDEs; from the simplest model equation, to hyperbolic systems, and equations with variable coefficients.
We have already introduced central concepts in Numerical Analysis, such as convergence, consistency, and stability for Finite Difference (FD) schemes.
We saw how these concepts are related by the Lax Equivalence Theprem.
1.1 The wave equation#
The full second-order two-way wave equation yields solutions propagating both ways; it describes a standing wavefield resulting from superposition of two waves in opposite directions. By formally “factoring” the differential operator
it is clear that solutions to either
are solutions to the original equation.
These are known as one-way wave equations or transport equations or advection equations describing a physical transport mechanism (with propagation speed \(a\) [Units of Length /Units of time]).
The simplest prototype for Hyperbolic PDEs is the one-way wave equation
where \(a\) is a constant, \(t\in\mathbb{R}^{+}\) represents time, and \(x\in\mathbb{R}\) the spatial location. The initial state, \(u(0,x)\), must be specified.
Once the initial condition (IC) \(u(0,x)=u_0(x)\) is specified, the unique solution to the one-way wave equation for \(t>0\) is given by
The solution at time \(t\) is just a shift of the initial value, \(u_0(x)\).
When \(a>0\) it is a shift to the right and when \(a<0\) it is a shift to the left.
The solution depends only on the value of \(\xi=x-at\).
These are lines in the \((t,x)\)-plane that are called characteristic lines or characteristics, and
hence \(a\) is the propagation speed.
This is typical for Hyperbolic Equations:
Note
The solution propagates with finite speed along characteristics.
We note that the exact solution
requires no differentiability of \(u\) (or \(u_0\)), whereas the equation
appears to only make sense if \(u\) is differentiable.
This picture shows a volume with low pressure near the rear of the aircraft at high subsonic airspeeds. [Source: Transonic speed regime Wiki page].
Hyperbolic equations feature solutions that are discontinuous (worse than non-differentiable); e.g., the sonic boom produced by an aircraft exceeding the speed of sound (Mach-1, or \(\approx\) 750 miles per hour at sea level) is an example of this phenomena.
The discontinuity in pressure, temperature, and density of the medium, in a propagating disturbance that moves faster than the local speed of sound in the medium, creates a shock wave.
[Source: Phys.org]
Devising numerical schemes which allow for discontinuous solutions requires “a bit” of ingenuity.
1.2 A more general hyperbolic equation#
Where \(a\) and \(b\) are constants.
We can introduce the following change of variables (and its inverse):
With \(\tilde{u}(\tau,\xi)=u(t,x)\), we can transform the PDE to an ODE along the characteristics:
The exact solution is given by
which expressed in the original variables is
With some work this method can be extended to nonlinear equations of the form
From a numerical point of view, the key thing to note is that the solution evolves with finite speed along the characteristics.
2. System of Hyperbolic Equations#
Now consider systems of hyperbolic equations with constant coefficients in one space dimension; \({\mathbf{u}}\) is now a \(d\)-dimensional vector (containing various quantities, e.g., density (\(\rho\)), pressure (\(p\)), velocity (\(v\)), energy (\(E\)), and linear momentum (\(\rho v\)) of a fluid or gas).
Definition: Hyperbolic System
A system of the form
is hyperbolic if the matrix \(A\) is diagonalizable with real eigenvalues.
The matrix \(A\) is diagonalizable, if there exists a non-singular matrix \(P\) such that
is a diagonal matrix.
The eigenvalues \(\lambda_1,\dots,\lambda_d\) are the characteristic speeds of the system.
In the easiest case, \(B=0\), we get
under the change of variables \(\mathbf{w} = P \mathbf{u}\). This is a reduction to \(d\) independent scalar hyperbolic equations.
When \(B\not=0\), the resulting system is coupled, but only in undifferentiated terms. The lower order term \(B\mathbf{u}\) causes growth, decay, or oscillations in the solution but does not alter the primary feature of solutions propagating along characteristics.
Example 8.1:#
with \(u(0,x)\equiv u_0=1\) if \(|x|\le1\), and 0 otherwise; and \(v(0,x) \equiv v_0=0\).
The eigenvalues are \(\lambda=\{3,1\}\), and without too much difficulty we find
Hence, the solution is
using Plots
using LinearAlgebra
using LaTeXStrings
default(linewidth=3)
default(legendfontsize=12)
u0(x) = abs.(x) .<= 1
# Spatial grid
x = -2:0.1:12
ctr = 0
for t in 0:0.5:3
# Compute u and v
u = (u0(x .- 3t) .+ u0(x .- t)) ./ 2
v = (u0(x .- 3t) .- u0(x .- t)) ./ 2
# Create subplots
p1 = plot(
x, u,
color = :blue,
label = L"u(t,x)",
ylim = (0, 1.05),
xlim = (-2, 12),
title = "t=$t"
)
p2 = plot(
x, v,
color = :red,
label = L"v(t,x)",
ylim = (-1.05/2, 1.05/2),
xlim = (-2, 12)
)
p = plot(p1, p2, layout = (2, 1))
display(p)
ctr += 1
end
3. Hyperbolic Equations with Variable Coefficients#
What happens when the propagation speed is variable, e.g.,
In this example the solution is still constant along characteristics, but the characteristics are not straight lines.
Here, we get an ODE for the \(x\)-coordinate
If, e.g., \(a(\tau,x)=x\), then \(x(\tau) = \xi e^{\tau}\) (so that \(\xi = xe^{-t}\)), and we get
3.1 Hyperbolic Systems with Variable Coefficients#
We can extend the definition of hyperbolicity to systems:
Definition
A system of the form
is diagonal with real eigenvalues and the matrix norms of \(P(t,x)\) and \(P^{-1}(t,x)\) are bounded in \(x\) and \(t\) for \(x\in \mathbb{R}\), \(t\ge0\).
4. Boundary Conditions#
Let’s solve the 1D wave equation on a finite domain, e.g., \([0,1]\).
We saw that if \(a\) is positive then the information propagates to the right and if \(a\) is negative it propagates to the left.
When \(a>0\), in addition to the initial value \(u(0,x)\) \(0\le x \le 1\), we must also specify the boundary value \(u(t,0)\) for all \(t>0\), and when \(a<0\) we must specify \(u(t,1)\) for \(t>0\).
The problem of determining a solution when both initial and boundary data are present is known as an Initial-Boundary Value Problem (IBVP).
Consider the hyperbolic system (assume \(a>0\), \(b>0\))
on the interval \(0\le x \le 1\). The characteristic speeds are \((a+b)\) and \((a-b)\), so that with \(w=u+v\), and \(z=u-v\)
If \(b<a\), then both characteristic speeds are positive, but when \(b>a\), we get one positive and one negative speed.
In the above figure, in the left panel \(b<a\), so both characteristics propagate to the right. In the right panel \(b>a\), so the characteristics propagate in opposite directions.
In order for the IBVPs to be well-posed we must
[Left panel:] specify the initial condition and two boundary conditions at \(x=0\);
[Right panel:] the initial condition, a boundary condition at \(x=0\), and a boundary condition at \(x=1\).
Note that the specified boundary conditions must be linearly independent from the outgoing (leaving the domain) characteristic.
5. FD schemes for hyperbolic equations#
5.1 FD schemes for the 1D wave equation#
For the 1D one-way wave equation or transport equation
we can write down a number of finite difference approximations at \((t_n,x_m)\), e.g.,
The schemes presented so far can all be written expressing \(u_k^{n+1}\) as linear combinations of \(u_{\mu}^{\nu}\) at previous time-levels \(\nu \in \{n-1,\,n\}\).
The FTFS scheme can be written as
where \(r = \Delta t/\Delta x\) is the ratio of the time- and space grid sizes.
This scheme is a one-step scheme since it only involves information from one previous time-level.
The leapfrog scheme (CTCS) is a two-step (multi-step) scheme.
For \(n\)-step schemes, we need to allow for initialization of the first \(n\) time-levels. Before we can apply an \(n\)-step scheme we must define \(u_k^0\), \(\dots\), \(u_k^{n-1}\), \(\forall k\).
5.2 Example using Leapfrog#
with initial condition
and boundary condition
Let’s try with different \(r = \Delta t/\Delta x\).
# -------------------------------
# Shared FD utilities (plus `diff1_mat` from Module 1-7)
# -------------------------------
"""First-derivative differentiation matrix on a uniform grid.
- interior: 2nd-order centered
- boundaries: 1st-order one-sided
"""
# Dfferentiation matrix for a first-derivative
function diff1_mat(x)
n = length(x)
D = zeros(n, n)
h = x[2] - x[1]
D[1, 1:2] = [-1/h 1/h] # Use a first-order forward difference at the left boundary
for i in 2:n-1
D[i, i-1:i+1] = [-1/2h 0 1/2h] # In the interior points, use a second-order centered difference
end
D[n, n-1:n] = [-1/h 1/h] # Use a first-order backward difference at the right boundary
D
end
"""Apply homogeneous Dirichlet BCs in-place."""
apply_dirichlet0!(u) = (u[1] = 0.0; u[end] = 0.0; u)
# -------------------------------
# Hyperbolic solvers as functions
# -------------------------------
"""Convenience: uniform grid spacing."""
Δx(x) = x[2] - x[1]
"""Leapfrog (CTCS) for linear advection u_t + a u_x = 0.
Returns (t, x, U) where U[n, i] ≈ u(t[n], x[i]).
"""
function leapfrog_advection(u0, a, x, T; r=0.8)
h = Δx(x)
dt = r * h
nt = round(Int, T/dt) + 1
t = range(0.0, step=dt, length=nt)
D = diff1_mat(x) # re-used from Module 1-7
A = -a * D # semi-discrete operator u_t = A*u
U = zeros(nt, length(x))
# first I.C. (zeroth-step)
U[1, :] .= u0.(x)
# first-step starter (exact shift is fine for this linear test)
U[2, :] .= u0.(x .- a*dt)
for n in 2:nt-1
U[n+1, :] .= U[n-1, :] .+ 2dt .* (A * U[n, :])
apply_dirichlet0!(U[n+1, :])
end
t, x, U
end
"""Lax–Friedrichs scheme for linear advection u_t + a u_x = 0.
"""
function lax_friedrichs_advection(u0, a, x, T; r=0.8)
h = Δx(x)
dt = r * h
nt = round(Int, T/dt) + 1
t = range(0.0, step=dt, length=nt)
U = zeros(nt, length(x))
U[1, :] .= u0.(x)
for n in 1:nt-1
U[n+1, 2:end-1] .= 0.5 .* U[n, 1:end-2] .+ 0.5 .* U[n, 3:end] .-
(a*r/2) .* (U[n, 3:end] .- U[n, 1:end-2])
apply_dirichlet0!(U[n+1, :])
end
t, x, U
end
lax_friedrichs_advection
# Solution to the 1D wave equation (advection equation) u_t + a u_x = 0 using the Leapfrog (CTCS) scheme
# constant wave speed
a = 1.0
# final time
T = 2.0
# I.C., a triangle wave
u0(x) = (1 .- abs.(x)) .* (abs.(x) .<= 1)
# grid
x = collect(-2.0:0.1:4.0)
t, x, U = leapfrog_advection(u0, a, x, T; r=0.8)
animation = @animate for k in eachindex(t)
plot(x, U[k, :];
marker = :circle,
label = "leapfrog",
xlims = (-2, 4),
ylims = (-0.1, 1.1),
grid = true,
title = "t = $(round(t[k]; digits=3))",
)
plot!(x, u0.(x .- a*t[k]); label = "exact")
end
gif(animation, "leapfrog.gif", fps=10, show_msg=false)
# plot only the final time
plot(x, U[end, :];
marker = :circle,
label = "leapfrog",
xlims = (-2, 4),
ylims = (-0.1, 1.1),
grid = true,
title = "Final time T = $T",
)
plot!(x, u0.(x .- a*T); label = "exact")
When we try with \( r > 1\), clearly something “strange” happens
5.3 Example using Lax-Friedrichs#
Note that the Lax-Friedrichs scheme is a one-step scheme, whereas the leapfrog scheme is a 2-step scheme.
Example 8.2:#
Same example as above:
with initial condition
and boundary condition
# u_t + a u_x = 0 using the Lax–Friedrichs scheme
# constant wave speed
a = 1.0
# final time
T = 2.0
# I.C.
u0(x) = (1 .- abs.(x)) .* (abs.(x) .<= 1)
# store final-time solutions for each D
xs = Dict{Int, Vector{Float64}}()
ss = Dict{Int, Vector{Float64}}()
r = 0.8
for D in (5, 10, 20)
Δx_D = 1.0 / D
x = collect(-2.0:Δx_D:4.0)
t, x, U = lax_friedrichs_advection(u0, a, x, T; r=r)
animation = @animate for k in eachindex(t)
plot(x, U[k, :];
marker = :circle,
label = "Lax–Friedrichs",
xlims = (-2, 4),
ylims = (-0.1, 1.1),
grid = true,
title = "D = $D, t = $(round(t[k]; digits=3))",
)
plot!(x, u0.(x .- a*t[k]); label="exact")
end
gif(animation, "lax_friedrichs_D$(D).gif", fps=10, show_msg=false)
xs[D] = x
ss[D] = vec(U[end, :])
end
# final comparison plot (at time T)
x5, s5 = xs[5], ss[5]
x10, s10 = xs[10], ss[10]
x20, s20 = xs[20], ss[20]
x_exact = x20
s_exact = u0.(x_exact .- a*T)
pfinal = plot(x5, s5; marker=:circle, label="Δx=1/5")
plot!(pfinal, x10, s10; marker=:diamond, label="Δx=1/10")
plot!(pfinal, x20, s20; marker=:x, label="Δx=1/20")
plot!(pfinal, x_exact, s_exact; label="Exact")
xlims!(pfinal, (0, 4))
ylims!(pfinal, (0, 1.1))
title!(pfinal, "Lax–Friedrichs r=$r (T=$T)")
display(pfinal)
From the example above we can see that as \(\Delta x \to 0\), the approximate solution gets better and better.
Consistency alone is necessary, but not sufficient for a Finite Difference Scheme (FDS) to be convergent (Recall Lax Equivalence theorem).
Definition (from Strikwerda)
A finite difference scheme \(P_{k,h} v_{m}^{n}=0\) for a first-order equation is stable in a stability region \(\Lambda\) if there is an integer \(S\) such that for any positive time \(T\), there is a constant \(C_T\) such that
for \(0 \le n \Delta t \le T\), with \((\Delta t,\Delta x)\in\Lambda\).
Recall that the quantity
is the \(L^2\) norm of the grid function \(w\), and is a measure of the size (energy) of the solution. The multiplication by \(h\equiv \Delta x\) is needed so that the norm is not sensitive to grid refinements (the number of points increase as \(h\to 0\)).
With this notation, the inequality in the definition can be written
The inequality expresses a limit (in terms of energy) of how much the solution can grow.
Typically \(S=(n-1)\) for an \(n\)-step scheme.
6. Stability of hyperbolic PDEs#
We now turn our attention to the key stability criterion for hyperbolic PDEs.
In the example above, we saw some numerical evidence of the leapfrog scheme (applied to \(v_t+av_x=0\), \(a=1\)) breaking down when \( r>1\).
The condition \(|a r|<1\) is necessary for stability of many explicit FDS.
Consistency implies that the solution of the PDE, if it is smooth, is an approximate solution of the finite difference scheme (FDS).
Convergence means that a solution of the FDS approximates a solution of the PDE.
Consistency is necessary, but not sufficient for a FDS to be convergent.
We illustrate this with an example.
Consistency \(\not\Rightarrow\) Convergence#
We consider the one-way wave equation with constant \(a=1\) propagation speed, and apply the forward-time forward-space (F.T.F.S.) scheme
We rewrite this:
Let the initial condition be given by:
Hence the exact solution is \(v(t,x) = v_0(x-t)\), i.e., a “hump” of height and width one, traveling to the right with unitary speed.
The IC for the FDS is discretized as:
Let us see what region the numerical scheme depends on.
The figure above illustrates how the exact solution propagates; it is one in the pink band/region, and zero outside the band.
The IC of the FD scheme says that the values are zeros everywhere, except at the four points: \(u_0^0\), \(u_{-1}^0\), \(u_{-2}^0\), and \(u_{-3}^0\), where it is \(1\).
The figure above illustrates how the FD scheme solution propagates. In particular we have that \(u_k^n\equiv0\), \(\forall k>0,\ n\ge0\). Hence, \(u_k^n \not\rightarrow v(t^n,x_k)\) for \((t^n,x_k)\) in the part of the band strictly in the right half plane - \(v\) is one there, but \(u_k^n\) is zero, no matter how much we refine the grid.
If a scheme is convergent, then as \(u_k^n\) converges to \(v(t,x)\), then \(u_k^n\) is bounded in some sense; this is the essence of stability.
For almost all schemes there are restrictions on the way \(\Delta x\) and \(\Delta t\) can be chosen so that the particular scheme is stable.
Definition: Stability region
A stability region is any bounded non-empty region of the first quadrant of \(\mathbb{R}^2\) that has the origin as an accumulation point.
An explicit scheme is a scheme that can be written as
Implicit schemes, where the sum may contain terms with \(n'=n+1\) on the RHS, will be discussed in the future.
6.1 The Courant-Friedrichs-Lewy Condition#
The following result covers all one-step schemes we have seen so far.
Theorem
For an explicit scheme for the hyperbolic equation
of the form
with \(r = \Delta t/\Delta x\) held constant, a necessary condition for stability is the Courant-Friedrichs-Lewy (CFL)condition,
For systems of equations for which \(\mathbf{u}\) is a vector and \(\alpha\), \(\beta\), and \(\gamma\) are matrices, we must have \(|a_i r|\le1\) for all eigenvalues \(a_i\) of the matrix \(A\).
In the above figure, we see an illustration of the CFL condition, with \( r=1\) held fixed.
The green triangle shows the region of dependence, i.e., what region influences \(u^{n+1}_k\) (actually only the three points at the base of the triangle contribute).
The blue arrow corresponds to a characteristic with speed \(a=0.75\), which carries information inside the region of dependence.
The red arrow, corresponding to a characteristic speed of \(a=1.5\) carries information from outside the region of dependence; this information cannot be captured by the scheme.
Courant, Friedrichs and Lewy also proved the following theorem.
Theorem
There are no explicit, unconditionally stable, consistent finite difference schemes for hyperbolic systems of partial differential equations.
One way of thinking about these theorems is to define the numerical speed of propagation as \( r^{-1} = \Delta x/\Delta t\), and note that a necessary condition for the stability of a scheme is
This guarantees that the FDS can propagate information (energy) at least as fast as the PDE.
\(r^{-1}\) is the “speed limit” on the grid; which explains why (with \(a=1\)), we saw the breakdown of the leapfrog scheme when \(r>1 \Leftrightarrow r^{-1} < 1\).
7. Conditional and unconditional stability#
In general, even for other classes of problems (not necessarily hyperbolic PDEs) in the mathematics literature, we can find a variety of definitions of stability.
Let’s consider a two-level difference scheme in matrix form:
which will generally be a difference scheme for solving an initial-value problem (IVP) on \(\mathbb{R}\), for a homogeneous linear PDE.
A stronger definition of stability is that
Definition (from Thomas)
The difference scheme \(\mathbf{u}^{n+1} = Q \mathbf{u}^{n} , n \geq 0\) is said to be stable with respect to the norm \(\| \cdot \|\) if there exists positive constants \(\Delta x_0\) and \(\Delta t_0\), and non-negative constants \(K\) and \(\beta\) such that
for \(0 \leq t = (n+1) \Delta t\), \(0 < \Delta x \leq \Delta x_0\) and \(0 < \Delta t \leq \Delta t_0\).
Remark:
With this defintion, the solution is allowed to grow with time, but not with the number of time steps.
Another common definition that is used is one that does not allow for exponential growth, as in:
This latter definition of stability implies that the solutions to the difference equation must be bounded.
Proposition
The difference scheme \(\mathbf{u}^{n+1} = Q \mathbf{u}^{n} , n \geq 0\) is stable with respect to the \(\| \cdot \|\) norm if and only if there exists positive constants \(\Delta x_0\) and \(\Delta t_0\) and non-negative constants \(K\) and \(\beta\) such that
for \(0 \leq t = (n+1) \Delta t\), \(0 < \Delta x \leq \Delta x_0\), and \(0 < \Delta t \leq \Delta t_0\).
Remarks:
We notet that in the Proposition above the norm \(\| Q^{n+1} \|\) is a matrix norm, which is induced by the norm defined on the space in which we are working.
For \(p=2\), the Euclidean \(\ell_2\) norm of a matrix \(\| A \|\) is the spectral norm, that is, it is the largest singular value of \(A\), which is the square root of the largest eigenvalue of the matrix \(A A^*\) (where \(A^*\) denotes the conjugate gradient). Hence \( \| A \|_2 = \sqrt{ \lambda_{\textrm{max}} (A^* A)} = \sigma_{\textrm{max}} (A) \).
As with convergence, it is difficult to prove stability directly. (We will develop a method to find a necessary condition later in the course).
Every time we can find a condition (such as the \( |a r|\le 1\) condition above in the CFL condition), we say that the scheme is conditionally stable.
In the case where no restrictions on the relationship between \(\Delta t\) and \(\Delta x\) are needed for stability, we say that the scheme is stable or unconditionally stable.
In general, it is not easy to find eigenvalues of a matrix to have a sense of its norm. But Gershgorin Circle Theorem helps us as it gives an upper bound (so a worst-case scenario estimate) on the eigenvalues:
Gershgorin Circle Theorem
Let \(A\) be a \(n \times n\) matrix with entries \((a_{ij})\). For \(i \in \{ 1, 2, \ldots, n \}\) let \(s_i\) be the sum of the absolute values of the non-diagonal entries in the \(s\)-th row \(\rho_s = \sum_{\substack{j = 1 \\j \neq s}}^n |a_{sj}|\), then each eigenvalue of \(A\), \(\lambda\), is such that