The wave equation

12.4. The wave equation#

Closely related to the advection equation is the wave equation,

(12.4.1)#

u_{t t} - c^{2} u_{x x} = 0.

This is our first PDE having a second derivative in time. As in the advection equation, $u (x, t) = ϕ (x - c t)$ is a solution of (12.4.1), but now so is $u (x, t) = ϕ (x + c t)$ for any twice-differentiable $ϕ$ (see Exercise 2). Thus, the wave equation supports advection in both directions simultaneously.

We will use $x \in [0, 1]$ and $t > 0$ as the domain. Because $u$ has two derivatives in $t$ and in $x$ , we need two boundary conditions. We will use the Dirichlet conditions

(12.4.2)#

u (0, t) = u (1, t) = 0, t \geq 0,

and two initial conditions,

(12.4.3)#

\begin{aligned} u (x, 0) & = f (x), 0 \leq x \leq 1, \\ u_{t} (x, 0) & = g (x), 0 \leq x \leq 1. \end{aligned}

One approach is to discretize both the $u_{t t}$ and $u_{x x}$ terms using finite differences:

\frac{1}{τ^{2}} (U_{i, j + 1} - 2 U_{i, j} + U_{i, j - 1}) = \frac{c^{2}}{h^{2}} (U_{i + 1, j} - 2 U_{i, j} + U_{i - 1, j}) .

This equation can be rearranged to solve for $U_{i, j + 1}$ in terms of values at time levels $j$ and $j - 1$ . Rather than pursue this method, however, we will turn to the method of lines.

First-order system#

In order to be compatible with the standard IVP solvers that we have encountered, we must recast (12.4.1) as a first-order system in time. Using our typical methodology, we would define $y = u_{t}$ and derive

(12.4.4)#

\begin{array}{r} \begin{aligned} u_{t} & = y, \\ y_{t} & = c^{2} u_{x x} . \end{aligned} \end{array}

However, there is another, less obvious option for reducing to a first-order system:

(12.4.5)#

\begin{aligned} u_{t} & = z_{x}, \\ z_{t} & = c^{2} u_{x} . \end{aligned}

This second form is appealing in part because it’s equivalent to Maxwell’s equations for electromagnetism. In the Maxwell form we typically replace the velocity initial condition in (12.4.3) with a condition on $z$ , which may be physically more relevant in some applications:

(12.4.6)#

\begin{aligned} u (x, 0) & = f (x), 0 \leq x \leq 1, \\ z (x, 0) & = g (x), 0 \leq x \leq 1. \end{aligned}

Because waves travel in both directions, there is no preferred upwind direction. This makes a centered finite difference in space appropriate. Before application of the boundary conditions, semidiscretization of (12.4.5) leads to

(12.4.7)#

\begin{array}{r} [\begin{array}{c} u^{'} (t) \\ z^{'} (t) \end{array}] = [\begin{array}{c} 0 & D_{x} \\ c^{2} D_{x} & 0 \end{array}] [\begin{array}{c} u (t) \\ z (t) \end{array}] . \end{array}

The boundary conditions (12.4.2) suggest that we should remove both of the end values of $u$ from the discretization, but retain all of the $z$ values. We use $w (t)$ to denote the vector of all the unknowns in the semidiscretization. When computing $w^{'} (t)$ , we extract the $u$ and $z$ components, and we use dedicated functions for padding with the zero end values or chopping off the zeros as necessary.

Demo 12.4.1

We solve the wave equation (12.4.5) with speed $c = 2$ , subject to (12.4.2) and initial conditions (12.4.6).

m = 200
x,Dₓ = FNC.diffcheb(m,[-1,1]);

The boundary values of $u$ are given to be zero, so they are not unknowns in the ODEs. Instead they are added or removed as necessary.

extend = v -> [0;v;0]
chop = u -> u[2:m];

The following function computes the time derivative of the system at interior points.

ode = function(w,c,t)
    u = extend(w[1:m-1])
    z = w[m:2m]
    dudt = Dₓ*z
    dzdt = c^2*(Dₓ*u)
    return [ chop(dudt); dzdt ]
end;

Our initial condition is a single hump for $u$ .

u_init = @. exp(-100*(x+0.5)^2)
z_init = -u_init
w_init = [ chop(u_init); z_init ];  

Because the wave equation is hyperbolic, we can use a nonstiff explicit solver.

IVP = ODEProblem(ode,w_init,(0.,2.),2)
w = solve(IVP,RK4());

We plot the results for the original $u$ variable only. Its interior values are at indices 1:m-1 of the composite $w$ variable.

t = range(0,2,length=80)
U = [extend(w(t)[1:m-1]) for t in t]
contour(x,t,hcat(U...)',color=:redsblues,clims=(-1,1),
    levels=24,xlabel=L"x",ylabel=L"t",title="Wave equation",
    right_margin=3Plots.mm)

[ Info: Saved animation to /Users/driscoll/Documents/fnc-julia/advection/figures/wave-boundaries.mp4

The original hump breaks into two pieces of different amplitudes, each traveling with speed $c = 2$ . They pass through one another without interference. When a hump encounters a boundary, it is perfectly reflected, but with inverted shape. At time $t = 2$ , the solution looks just like the initial condition.

Variable speed#

An interesting situation is when the wave speed $c$ changes discontinuously, as when light passes from one material into another. For this we must replace the term $c^{2}$ in (12.4.7) with the matrix $diag (c^{2} (x_{0}), \dots, c^{2} (x_{m}))$ .

Demo 12.4.2

We now use a wave speed that is discontinuous at $x = 0$ ; to the left, $c = 1$ , and to the right, $c = 2$ . The ODE implementation has to change slightly.

ode = function(w,c,t)
    u = extend(w[1:m-1])
    z = w[m:2m]
    dudt = Dₓ*z
    dzdt = c.^2 .* (Dₓ*u)
    return [ chop(dudt); dzdt ]
end;

The variable wave speed is passed as an extra parameter through the IVP solver.

c = @. 1 + (sign(x)+1)/2
IVP = ODEProblem(ode,w_init,(0.,5.),c)
w = solve(IVP,RK4());

t = range(0,5,length=80)
U = [extend(w(t)[1:m-1]) for t in t]
contour(x,t,hcat(U...)',color=:redsblues,clims=(-1,1),
    levels=24,xlabel=L"x",ylabel=L"t",title="Wave equation",
    right_margin=3Plots.mm)

[ Info: Saved animation to /Users/driscoll/Documents/fnc-julia/advection/figures/wave-speed.mp4

Each pass through the interface at $x = 0$ generates a reflected and transmitted wave. By conservation of energy, these are both smaller in amplitude than the incoming bump.

Exercises#

✍ Consider the Maxwell equations (12.4.5) with smooth solution $u (x, t)$ and $z (x, t)$ .

(a) Show that $u_{t t} = c^{2} u_{x x}$ .

(b) Show that $z_{t t} = c^{2} z_{x x}$ .
✍ Suppose that $ϕ (s)$ is any twice-differentiable function.

(a) Show that $u (x, t) = ϕ (x - c t)$ is a solution of $u_{t t} = c^{2} u_{x x}$ . (As in the advection equation, this is a traveling wave of velocity $c$ .)

(b) Show that $u (x, t) = ϕ (x + c t)$ is another solution of $u_{t t} = c^{2} u_{x x}$ . (This is a traveling wave of velocity $- c$ .)
✍ Show that the following is a solution to the wave equation $u_{t t} = c^{2} u_{x x}$ with initial and boundary conditions (12.4.2) and (12.4.3):

$u (x, t) = \frac{1}{2} [f (x - c t) + f (x + c t)] + \frac{1}{2 c} \int_{x - c t}^{x + c t} g (ξ) d ξ$

This is known as D’Alembert’s solution.
⌨ Suppose the wave equation has homogeneous Neumann conditions on $u$ at each boundary instead of Dirichlet conditions. Using the Maxwell formulation (12.4.5), we have $z_{t} = c^{2} u_{x}$ , so $z$ is constant in time at each boundary. Therefore, the endpoint values of $z$ can be taken from the initial condition and removed from the ODE, while the entire $u$ vector is now part of the ODE.

Modify Demo 12.4.1 to solve the PDE there with Neumann instead of Dirichlet conditions, and make an animation or space-time portrait of the solution. In what major way is it different from the Dirichlet case?
⌨ The equations $u_{t} = z_{x} - σ u$ , $z_{t} = c^{2} u_{x x}$ model electromagnetism in an imperfect conductor. Repeat Demo 12.4.2 with $σ (x) = 2 + 2 sign (x)$ . (This causes waves in the half-domain $x > 0$ to decay in time.)
The nonlinear sine–Gordon equation $u_{t t} - u_{x x} = \sin u$ has interesting solutions.

(a) ✍ Write the equation as a first-order system in the variables $u$ and $v = u_{t}$ .

(b) ⌨ Assume periodic end conditions on $[- 10, 10]$ and discretize at $m = 200$ points. Let $u (x, 0) = π e^{- x^{2}}$ and $u_{t} (x, 0) = 0$ . Solve the system using RK4 between $t = 0$ and $t = 50$ , and make a plot or animation of the solution.
The deflections of a stiff beam, such as a ruler, are governed by the PDE $u_{t t} = - u_{x x x x}$ .

(a) ✍ Show that the beam PDE is equivalent to the first-order system

$\begin{aligned} u_{t} & = v_{x x}, \\ v_{t} & = - u_{x x} . \end{aligned}$

(b) ⌨ Assuming periodic end conditions on $[- 1, 1]$ , use $m = 100$ , let $u (x, 0) = \exp (- 24 x^{2})$ , $v (x, 0) = 0$ , and simulate the solution of the beam equation for $0 \leq t \leq 1$ using Function 6.7.3 with $n = 100$ time steps. Make a plot or animation of the solution.

The wave equation

Contents

12.4. The wave equation#

First-order system#

Variable speed#

Exercises#