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",
    right_margin=3Plots.mm)
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 
│   fn = /Users/driscoll/repos/fnc-julia/advection/figures/wave-boundaries.mp4
└ @ 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",
    right_margin=3Plots.mm)
Wave equation