# The wave equation

## Contents

using FundamentalsNumericalComputation
FNC.init_format()

┌ Info: verify download of index files...
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/MatrixDepot.jl:139
┌ Info: writing database
┌ Info: used remote sites are sparse.tamu.edu with MAT index and math.nist.gov with HTML index
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/MatrixDepot.jl:141


# 12.4. The wave equation¶

Closely related to the advection equation is the wave equation,

(12.4.1)$u_{tt} - c^2 u_{xx} = 0.$

This is our first PDE having a second derivative in time. As in the advection equation, $$u(x,t)=\phi(x-ct)$$ is a solution of (12.4.1), but now so is $$u(x,t)=\phi(x+c t)$$ for any twice-differentiable $$\phi$$ (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, \qquad t \ge 0,$

and two initial conditions,

(12.4.3)$\begin{split}u(x,0) &= f(x), \qquad 0 \le x \le 1, \\ u_t(x,0) &= g(x), \qquad 0 \le x \le 1. \end{split}$

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

$\frac{1}{\tau^2}(U_{i,j+1} - 2U_{i,j} + U_{i,j-1}) = \frac{c^2}{h^2} (U_{i+1,j} - 2U_{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{split}\begin{split} u_t &= y, \\ y_t &= c^2 u_{xx}. \end{split}\end{split}$

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

(12.4.5)$\begin{split} u_t &= z_x, \\ z_t &= c^2 u_{x}.\end{split}$

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{split}u(x,0) &= f(x), \qquad 0 \le x \le 1, \\ z(x,0) &= g(x), \qquad 0 \le x \le 1.\end{split}$

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{split} \begin{bmatrix} \mathbf{u}'(t) \\[2mm] \mathbf{z}'(t) \end{bmatrix} = \begin{bmatrix} \boldsymbol{0} & \mathbf{D}_x \\[2mm] c^2 \mathbf{D}_x & \boldsymbol{0} \end{bmatrix} \begin{bmatrix} \mathbf{u}(t) \\[2mm] \mathbf{z}(t) \end{bmatrix}.\end{split}$

The boundary conditions (12.4.2) suggest that we should remove both of the end values of $$\mathbf{u}$$ from the discretization, but retain all of the $$\mathbf{z}$$ values. We use $$\mathbf{w}(t)$$ to denote the vector of all the unknowns in the semidiscretization. When computing $$\mathbf{w}'(t)$$, we extract the $$\mathbf{u}$$ and $$\mathbf{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 $$\mathbf{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")

anim = @animate for t in range(0,2,length=120)
plot(x,extend(w(t)[1:m-1]),label=@sprintf("t=%.3f",t),
xaxis=(L"x"),yaxis=([-1,1],L"u(x,t)"),dpi=100,
title="Wave equation")
end
mp4(anim,"figures/wave-boundaries.mp4")

┌ Info: Saved animation to
└ @ Plots /Users/driscoll/.julia/packages/Plots/8K4be/src/animation.jl:114


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 $$\operatorname{diag}\bigl(c^2(x_0),\ldots,c^2(x_m)\bigr)$$.

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