using FundamentalsNumericalComputation
┌ Info: verify download of index files...
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/MatrixDepot.jl:139
┌ Info: reading database
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:23
┌ Info: adding metadata...
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:67
┌ Info: adding svd data...
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:69
┌ Info: writing database
└ @ MatrixDepot /Users/driscoll/.julia/packages/MatrixDepot/GEDc3/src/download.jl:74
┌ Info: used remote sites are with MAT index and 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 ]

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]
    levels=24,xlabel=L"x",ylabel=L"t",title="Wave equation")
Wave equation
anim = @animate for t in range(0,2,length=120)
        title="Wave equation")
┌ 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 ]

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]
    levels=24,xlabel=L"x",ylabel=L"t",title="Wave equation")
Wave equation