The wave equation - Fundamentals of Numerical Computation

A close cousin of the advection equation is the wave equation.

The wave equation describes the propagation of waves, such as for sound or light. This is our first PDE having a second derivative in time. Consequently, we should expect to have two initial conditions:

\begin{split} u(x,0) &= f(x), \\ u_t(x,0) &= g(x). \end{split}

(12.4.2)

There are also two space derivatives, so if boundaries are present, two boundary conditions are required.

In the absence of boundaries, $u(x,t)=\phi(x-ct)$ is a solution of (12.4.1), just as it is for the advection equation. Now, though, so is $u(x,t)=\phi(x+c t)$ for any twice-differentiable ϕ (see Exercise 12.4.2).

12.4.1First-order system¶

In order to apply semidiscretization 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

\begin{split} u_t &= y, \\ y_t &= c^2 u_{xx}. \end{split}

(12.4.3)

This is a valid approach. However, the wave equation has another, less obvious option for transforming to a first-order system:

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

(12.4.4)

This second form is appealing in part because it’s equivalent to Maxwell’s equations for electromagnetism. In this form, we typically replace the velocity initial condition in (12.4.2) with a condition on $z$ ,

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

(12.4.5)

This alternative to (12.4.2) is useful in many electromagnetic applications, where $u$ and $z$ represent the electric and magnetic fields, respectively.

12.4.2Semidiscretization¶

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.4) leads to

\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}.

(12.4.6)

We now suppose that the domain is $[0,1]$ and impose the Dirichlet conditions

u(0,t) = u(1,t) = 0, \qquad t \ge 0.

(12.4.7)

The boundary conditions (12.4.7) suggest that we should remove both of the end values of $\mathbf{u}$ from the discretization, but retain all the $\mathbf{z}$ values. We use $\mathbf{w}(t)$ to denote the vector of all the unknowns in the semidiscretization. If the nodes are numbered $x_0,\ldots,x_m$ , then we have

\mathbf{w}(t) = \begin{bmatrix} u_1(t) \\ \vdots \\ u_{m-1}(t) \\ z_0(t) \\ \vdots \\ z_m(t) \end{bmatrix} = \begin{bmatrix} \mathbf{v}(t) \\[1mm] \mathbf{z}(t) \end{bmatrix} \in \mathbb{R}^{2m}.

(12.4.8)

When computing $\mathbf{w}'(t)$ , we extract the $\mathbf{v}$ and $\mathbf{z}$ components, and we define two helper functions: extend, which pads the $\mathbf{v}$ component with the zero end values, and chop, which deletes them from $\mathbf{u}$ to give $\mathbf{v}$ .

Example 12.4.1 (Wave equation with boundaries)

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

Julia

MATLAB

Python

Example 12.4.1

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]
    du_dt = Dₓ * z
    dz_dt = c^2 * (Dₓ * u)
    return [ chop(du_dt); dz_dt ]
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.

using OrdinaryDiffEq
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.

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

Source

anim = @animate for t in range(0 ,2, 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=150,  title="Wave equation")
end
mp4(anim, "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.

Example 12.4.1

c = 2;  m = 200;
[x, Dx] = 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.

chop = @(u) u(2:m);
extend = @(v) [0; v; 0];

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

f124wave.m

function dw_dt = wave(t, w, param)
    [c, m, Dx, chop, extend] = param{:};
    u = extend(w(1:m-1));
    z = w(m:2*m);
    du_dt = Dx * z;
    dz_dt = c.^2 .* (Dx * u);
    dw_dt = [ chop(du_dt); dz_dt ];
end

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

u_init = exp( -100*x.^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.

odefun = @(t, w) f124wave(t, w, {c, m, Dx, chop, extend});
t = linspace(0, 2, 101);
[t, W] = ode45(odefun, t, w_init);
W = W';

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.

n = length(t)-1;
U = [ zeros(1, n+1); W(1:m-1, :); zeros(1, n+1) ];

cmap = zeros(256, 3);
cmap(129:end, 1) = 2/3;
cmap(1:128, 2) = linspace(1, 0, 128);
cmap(129:256, 2) = linspace(0, 1, 128);
cmap(1:128, 3) = linspace(1, 0.5, 128);
cmap(129:256, 3) = linspace(0.5, 1, 128);
cmap = hsv2rgb(cmap);

clf,  contour(x, t, U', 24, linewidth=2)
colormap(cmap),  clim([-1, 1])
xlabel x,  ylabel t
title("Wave equation with boundaries")

Source

clf
plot(x, U(:, 1))
hold on
axis([-1, 1, -1.05, 1.05])
title("Wave equation with boundaries") 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/wave-boundaries.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = 1:length(t)
    cla, plot(x, U(:, frame))
    str = sprintf("t = %.2f", t(frame));
    text(-0.92, 0.85, str, fontsize=16);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

Example 12.4.1

m = 200
x, Dx, Dxx = FNC.diffmat2(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.

chop = lambda u: u[1:-1]
extend = lambda v: hstack([0, v, 0])

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

def dw_dt(t, w):
    u = extend(w[:m-1])
    z = w[m-1:]
    du_dt = Dx @ z
    dz_dt = c**2 * (Dx @ u)
    return hstack([chop(du_dt), dz_dt])

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

u_init = exp(-100 * x**2)
z_init = -u_init
w_init = hstack([chop(u_init), z_init])

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

c = 2
sol = solve_ivp(dw_dt, (0, 2), w_init, dense_output=True)
u = lambda t: extend(sol.sol(t)[:m-1])   # extract the u component

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 = linspace(0, 2, 80)
U = [u(tj) for tj in t]
contour(x, t, U, levels=24, cmap="RdBu", vmin=-1, vmax=1)
xlabel("$x$"),  ylabel("$t$")
title("Wave equation with boundaries");

Source

from matplotlib import animation
fig, ax = subplots()
curve = ax.plot(x, u_init)[0]
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ax.set_xlabel("$x$")
ax.set_ylabel("$u(x,t)$")
ax.set_ylim(-1.05, 1.05)
ax.set_title("Wave equation with boundaries")

def snapshot(t):
    curve.set_ydata(u(t))
    time_text.set_text(f"t = {t:.2f}")

anim = animation.FuncAnimation(
    fig, snapshot, frames=linspace(0, 2, 161)
    )
anim.save("figures/wave-boundaries.mp4", fps=30)
close()

12.4.3Variable speed¶

An interesting situation is when the wave speed $c$ in (12.4.1) changes discontinuously, as when light passes from one material into another. For this we must replace the term $c^2$ in (12.4.6) with the matrix $\operatorname{diag}\bigl(c^2(x_0),\ldots,c^2(x_m)\bigr)$ .

Example 12.4.2 (Wave equation with variable speed)

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

Julia

MATLAB

Python

Example 12.4.2

The ODE implementation has to change slightly.

ode = function(w,c,t)
    u = extend(w[1:m-1])
    z = w[m:2m]
    du_dt = Dₓ*z
    dz_dt = c.^2 .* (Dₓ*u)
    return [ chop(du_dt); dz_dt ]
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, 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
    )

Source

anim = @animate for t in range(0,5,181)
    plot(Shape([-1, 0, 0, -1], [-1, -1, 1, 1]), color=RGB(.8, .8, .8), l=0, label="")
    plot!(x, extend(w(t, idxs=1:m-1));
        label=@sprintf("t=%.2f", t), 
        xaxis=(L"x"),  yaxis=([-1, 1], L"u(x,t)"),
        dpi=150,  title="Wave equation, variable speed")
end
mp4(anim, "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.

Example 12.4.2

We now use a wave speed that is discontinuous at $x=0$ .

m = 120;
[x, Dx] = diffcheb(m, [-1, 1]);
c = 1 + (sign(x) + 1) / 2;
chop = @(u) u(2:m);
extend = @(v) [0; v; 0];

u_init = exp( -100*(x + 0.5).^2 );
z_init = -u_init;
w_init = [ chop(u_init); z_init ];  
odefun = @(t, w) f124wave(t, w, {c, m, Dx, chop, extend});
t = linspace(0, 2, 101);
[t, W] = ode45(odefun, t, w_init);
W = W';
U = [ zeros(1, n+1); W(1:m-1, :); zeros(1, n+1) ];

clf,  contour(x, t, U', 24, linewidth=2)
colormap(cmap),  clim([-1, 1])
xlabel x,  ylabel t
title("Wave equation with variable speed")

Source

clf
plot(x, U(:, 1))
hold on
axis([-1, 1, -1.05, 1.05])
title("Wave equation with variable speed") 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/wave-speed.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = 1:length(t)
    cla, plot(x, U(:, frame))
    str = sprintf("t = %.2f", t(frame));
    text(-0.92, 0.85, str, fontsize=16);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

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.

Example 12.4.2

The variable wave speed is set to be re-used

m = 120
x, Dx, Dxx = FNC.diffcheb(m, [-1, 1])
c = 1 + (sign(x) + 1) / 2
u_init = exp(-100 * x**2)
z_init = -u_init
w_init = hstack([chop(u_init), z_init])

sol = solve_ivp(dw_dt, (0, 5), w_init, dense_output=True, method="Radau")
u = lambda t: extend(sol.sol(t)[:m-1])

t = linspace(0, 5, 150)
U = [u(tj) for tj in t]
contour(x, t, U, levels=24, cmap="RdBu", vmin=-1, vmax=1)
xlabel("$x$"),  ylabel("$t$")
title("Wave equation with variable speed");

Source

fig, ax = subplots()
curve = ax.plot(x, u_init)[0]
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ax.set_xlabel("$x$")
ax.set_ylabel("$u(x,t)$")
ax.set_ylim(-1.05, 1.05)
ax.set_title("Wave equation with variable speed")

anim = animation.FuncAnimation(
    fig, snapshot, frames=linspace(0, 5, 251)
    )
anim.save("figures/wave-speed.mp4", fps=30)
close()

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.

12.4.4Exercises¶

Exercise 12.4.6

The vibrations of a stiff beam, such as a ruler, are governed by the PDE $u_{tt}=-u_{xxxx}$ , where $u(x,t)$ is the deflection of the beam at position $x$ and time $t$ . This is called the dynamic beam equation.

(a) ✍ Show that solutions of the first-order system

\begin{align*} u_t &= v_{xx}, \\ v_t &= -u_{xx}. \end{align*}

(12.4.11)

also satisfy the dynamic beam equation.

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