Numerical Methods for Waves

Learning Objectives

By the end of this section you will be able to:

1. Replace the continuous wave PDE $\,u_{tt} = c^{2} u_{xx}\,$ with its discrete cousin on a space-time grid.
2. Derive the explicit five-point stencil $\,u_i^{n+1} = 2u_i^{n} - u_i^{n-1} + r^{2}\bigl(u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n}\bigr)\,$ and read it as "curvature in space drives curvature in time."
3. Prove and feel the CFL bound $\,r = c\,\Delta t/\Delta x \le 1\,$ — the one inequality that decides whether your simulation is a wave or an explosion.
4. Handle Dirichlet, Neumann, and absorbing boundary conditions, and explain physically what each one does to a reflected pulse.
5. Implement the scheme in plain Python first, then translate it to vectorized PyTorch so it runs on a GPU and plugs into machine-learning pipelines.

Why Go Numerical?

"Analytic solutions are gifts. Numerical solutions are tools."

The previous four sections solved the wave equation by hand — d'Alembert's travelling waves, separation of variables, Fourier series, Bessel functions. Those methods are beautiful, but they only work on clean geometries with clean initial data. Real problems rarely cooperate:

The cure is to discretize: replace the continuous space-time with a finite grid, replace each derivative with a finite difference, and march the solution forward in time. The result is the Finite Difference Method (FDM) — the oldest, simplest, and arguably most pedagogically valuable way to solve a PDE on a computer.

Intuition: imagine the membrane of a drum painted with a grid. At every intersection we track a single number, the height of that pinpoint at a single instant. We then step time forward in tiny ticks. To know tomorrow's height of a point, we only need today's heights of that point and its immediate neighbors. Local rule, global wave. That is the whole game.

Discretizing Space and Time

Pick a spatial step $\Delta x$ and a temporal step $\Delta t$. Replace the continuous function $u(x, t)$ by samples on a lattice:

The superscript $n$ is the time row, the subscript $i$ is the spatial column. Think of$\,u_i^n\,$ as a giant matrix where each column is one moving point of the string and each row is one snapshot in time.

Approximating the derivatives

Taylor expansion of $u$ around a grid point gives, after the linear terms cancel, the famous second-order central differences:

Each formula has the same shape: the second derivative of a smooth function at a point is, up to an error of order $h^{2}$, the difference between the value there and the average of its two neighbors, scaled by $1/h^{2}$. That is geometry: it measures how much the curve bulges away from the straight line through its neighbors. That bulge is curvature, and curvature is what drives the wave.

The discrete wave equation

Plug the two central differences into $u_{tt} = c^{2} u_{xx}$ at the grid point $(i, n)$:

Solve for the only unknown, $u_i^{n+1}$, and introduce the CFL number $r = c\,\Delta t / \Delta x$:

This is the only formula you have to remember. Everything else — stability, accuracy, boundary handling — spins out of this one line.

The Five-Point Stencil

The update rule reads from five grid points: the three neighbors at the current time row, the same point one step in the past, and writes one new value in the future. Drawn on a space-time grid it looks like a plus sign:

Two intuitions sit inside this stencil:

• Spatial part — curvature drives acceleration. The bracket $u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n}$ is the discrete Laplacian: positive when the point is below its neighbors, negative when above. The tension wants to pull the membrane back toward the average, and the wave equation says this force is proportional to acceleration.
• Temporal part — momentum. The combination $2u_i^{n} - u_i^{n-1}$ is exactly where the point would be next if there were no force at all (constant velocity extrapolation). It is the discrete version of inertia.

Add inertia to the curvature-driven kick and you get the next time step. That is Newton's second law on a grid.

The First-Step Trick

The stencil needs two previous time rows, but at $n = 0$ we only have one. We need a special recipe for the first step.

Use the initial-velocity data $u_t(x, 0) = g(x)$ and a Taylor expansion in time:

Substitute the wave equation $u_{tt}(x, 0) = c^{2}\,u_{xx}(x, 0)$ and replace $u_{xx}$ by central differences:

The half in front of $r^{2}$ is the signature of this special first step — miss it and the whole scheme silently drops to first-order accuracy. Released-from-rest data means $g(x) \equiv 0$ and the middle term vanishes.

Worked Example: Step Through By Hand

Take the classic plucked string setup: $u_{tt} = u_{xx}$ on $[0, 1]$ with $u(x, 0) = \sin(\pi x)$, $u_t(x, 0) = 0$ and fixed ends. The exact solution is $u(x, t) = \sin(\pi x)\cos(\pi t)$ — the fundamental mode breathing at frequency one. We will compute three time steps by hand and compare to the exact answer.

The CFL Stability Condition

Why $r \le 1$? The physical intuition is beautiful: information in the wave equation travels at speed $c$. In one numerical time step $\Delta t$, the physical wave moves a distance $c\,\Delta t$. The numerical stencil gathers information from one grid spacing $\Delta x$ away.

Courant's rule: the numerical "sight" of one cell per step must reach at least as far as the physical wave front travels. Otherwise the scheme is trying to predict a value the wave has already passed through — information arrives that the stencil has not even looked at yet, and the simulation goes berserk.

Mathematically the same condition pops out of von Neumann stability analysis. Substitute the Fourier mode $u_j^{n} = \xi^{n} e^{ikj\Delta x}$ into the update rule. You get a quadratic in the growth factor $\xi$:

For every wavenumber $k$, both roots have modulus $|\xi| = 1$ exactly when the coefficient stays in $[-2, 2]$, which forces $r \le 1$. Equality, the so-called magic CFL number, is even better: it makes the FDM update exactly equal to d'Alembert's solution — no dispersion at all.

Push the slider just past $r = 1$ and the Gaussian pulse explodes within a few dozen steps. Pull it back and watch the wave snap right back into clean propagation. The transition is razor sharp because the underlying inequality is binary — you are either inside the unit circle or you are not.

Boundary Conditions

The interior stencil has the same shape everywhere, but at $i = 0$ and $i = N$ we need to say something about what is beyond the edge. Three flavors cover essentially every classical problem.

The physics: a fixed end is a wall — force the displacement to zero, the wave hands its energy back inverted. A free end is a ring on a frictionless rod — no transverse force, so the slope must be flat. An absorbing end is an imaginary infinite extension — we model what a far-field observer would see at the truncation point.

Interactive 1D Wave Simulator

Time to play. The next panel runs the full FDM scheme live in your browser. Switch initial conditions, slide the CFL number into the danger zone, swap boundary conditions, and observe what each knob does to the physics.

Things worth trying:

• With sine mode + fixed BC, turn on "show exact". At $r = 1$ the cyan curve sits exactly on top of the yellow dashed line forever — the magic CFL number.
• Launch a Gaussian with absorbing BCs. The pulse escapes through the boundaries instead of bouncing.
• With two pulses + fixed BCs, watch them collide, pass through each other (linearity!), and rebound off the walls inverted.
• Push the CFL slider past 1. Watch the wave turn red and detonate.

Plain Python Implementation

Now that you have the rule in muscle memory, here is the whole solver as roughly thirty lines of NumPy. Read the explanations and feel how each line maps to a physical idea.

np.array([1,2,3]) + np.array([4,5,6]) -> [5,7,9]

solve_wave_1d(c=2.0, r=0.8, boundary='absorbing')

L=1, nx=5 -> dx = 0.25

r=1.0 means dt = dx/c (sharpest stable scheme)

np.linspace(0,1,5) -> [0, 0.25, 0.5, 0.75, 1.0]

If g(x)=0 (released from rest), the dt·g term drops out.

If u_curr = [0, 0.67, 0.95, 0.67, 0] and u_prev = [0, 0.71, 1, 0.71, 0], with r²=0.16, the interior of u_next is [0.578, 0.817, 0.578]

When r=1 the formula collapses to u_next[N-1] = u_curr[N-2]: the right boundary just copies its inner neighbor (perfectly transparent).

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_wave_1d(c=1.0, L=1.0, T=2.0, nx=101, r=0.9,
5                  f=lambda x: np.sin(np.pi * x),
6                  g=lambda x: np.zeros_like(x),
7                  boundary="fixed"):
8    """Finite-difference solver for u_tt = c^2 u_xx on [0, L] x [0, T]."""
9    dx = L / (nx - 1)
10    dt = r * dx / c                    # CFL: r = c*dt/dx <= 1
11    nt = int(np.ceil(T / dt)) + 1
12    x  = np.linspace(0.0, L, nx)
13    r2 = r * r
14
15    u_prev = f(x).astype(float)        # u at time level n-1
16    u_curr = np.empty_like(u_prev)     # u at time level n
17    u_next = np.empty_like(u_prev)     # u at time level n+1
18
19    # Special first step uses the initial velocity g(x)
20    u_curr[1:-1] = (u_prev[1:-1]
21        + dt * g(x[1:-1])
22        + 0.5 * r2 * (u_prev[2:] - 2 * u_prev[1:-1] + u_prev[:-2]))
23    u_curr[0] = u_curr[-1] = 0.0
24
25    history = [u_prev.copy(), u_curr.copy()]
26
27    for n in range(1, nt - 1):
28        u_next[1:-1] = (2 * u_curr[1:-1] - u_prev[1:-1]
29            + r2 * (u_curr[2:] - 2 * u_curr[1:-1] + u_curr[:-2]))
30
31        if boundary == "fixed":        # Dirichlet
32            u_next[0]  = 0.0
33            u_next[-1] = 0.0
34        elif boundary == "free":       # Neumann (zero-derivative)
35            u_next[0]  = u_next[1]
36            u_next[-1] = u_next[-2]
37        elif boundary == "absorbing":  # Mur 1st-order open BC
38            u_next[0]  = u_curr[1]  + ((r - 1) / (r + 1)) * (u_next[1]  - u_curr[0])
39            u_next[-1] = u_curr[-2] + ((r - 1) / (r + 1)) * (u_next[-2] - u_curr[-1])
40
41        u_prev, u_curr, u_next = u_curr, u_next, u_prev
42        history.append(u_curr.copy())
43
44    return x, np.array(history), dt
45
46x, U, dt = solve_wave_1d(nx=101, r=0.9)
47print(f"steps={U.shape[0]}, dt={dt:.4f}")
48print(f"max|u| over time = {np.max(np.abs(U)):.4f}")

Reading the algorithm in code form makes the trade-off explicit: we traded a clean closed-form Fourier series for a vectorized 3-line update we can run on anything — arbitrary c(x), arbitrary domain shapes, arbitrary forcing.

Vectorized PyTorch Implementation

The NumPy version maxes out a single CPU core. For 2D / 3D wave problems, batch simulation, or for plugging the solver into a neural network, we rewrite the same code in PyTorch. The pleasant surprise: only the import line really changes.

Drop the decorator if you want gradients of T(u₀) for adjoint optimization.

u_prev = exp(-200·(x-0.5)²) for a pulse instead of a sine

torch.nn.functional.conv1d(u.view(1,1,-1), torch.tensor([[[1,-2,1]]]))

1import torch
2
3@torch.no_grad()
4def solve_wave_1d_torch(c=1.0, L=1.0, T=2.0, nx=2049, r=0.9,
5                         device="cuda" if torch.cuda.is_available() else "cpu",
6                         dtype=torch.float32):
7    """Same FDM scheme as the NumPy version, but as a GPU-friendly torch.Tensor."""
8    dx = L / (nx - 1)
9    dt = r * dx / c
10    nt = int((T / dt)) + 1
11    r2 = r * r
12
13    x = torch.linspace(0.0, L, nx, device=device, dtype=dtype)
14    u_prev = torch.sin(torch.pi * x).clone()
15    u_curr = torch.empty_like(u_prev)
16    u_next = torch.empty_like(u_prev)
17
18    # First step (g(x)=0 here, so the dt*g term vanishes)
19    u_curr[1:-1] = (u_prev[1:-1]
20                    + 0.5 * r2 * (u_prev[2:] - 2 * u_prev[1:-1] + u_prev[:-2]))
21    u_curr[0] = u_curr[-1] = 0.0
22
23    for _ in range(1, nt - 1):
24        u_next[1:-1] = (2 * u_curr[1:-1] - u_prev[1:-1]
25                        + r2 * (u_curr[2:] - 2 * u_curr[1:-1] + u_curr[:-2]))
26        u_next[0] = u_next[-1] = 0.0
27        u_prev, u_curr, u_next = u_curr, u_next, u_prev
28
29    return x.cpu(), u_curr.cpu()
30
31# Conv1d view — the Laplacian is literally a kernel convolution
32def laplacian_kernel(device, dtype):
33    return torch.tensor([[[1.0, -2.0, 1.0]]], device=device, dtype=dtype)
34
35def step_conv(u_prev, u_curr, kernel, r2):
36    """One time step expressed as a 1-D convolution. Drop-in batched version."""
37    lap = torch.nn.functional.conv1d(u_curr.view(1, 1, -1), kernel).view(-1)
38    u_next = 2 * u_curr[1:-1] - u_prev[1:-1] + r2 * lap
39    out = torch.zeros_like(u_curr)
40    out[1:-1] = u_next
41    return out

That last function is the bridge to modern research. Once your FDM time step is a $\texttt{conv1d}$, the boundary between "classical PDE solver" and "learnable neural network" collapses. Differentiable solvers, neural operators, and physics-informed networks all start exactly here.

Accuracy, Dispersion, and Error

The scheme is second-order accurate in both $\Delta x$ and $\Delta t$: halving the grid roughly quarters the error.

Numerical dispersion

A subtler issue is numerical dispersion. The exact wave equation has every wavenumber propagating at the same speed $c$. The discrete scheme has a slightly wavenumber-dependent speed:

Long-wavelength components (small $k$) travel at almost $c$; short ones lag. The net effect is that a sharp Gaussian pulse spreads out slightly with time — energy at high $k$ gets left behind. At the magic value $r = 1$ the formula collapses to $c_{\text{num}} = c$ identically — no dispersion, full agreement with d'Alembert.

How to keep the error small in practice

1. Pick the finest $\Delta x$ your memory allows.
2. Pick the largest $\Delta t$ that still satisfies $r \le 1$. Smaller $\Delta t$ is not always better — it actually increases dispersion.
3. Resolve the shortest physical wavelength with at least 10–20 grid points.
4. If you need very long-time integration, switch to a higher-order scheme (Crank–Nicolson, spectral, or DG).

Beyond Basic FDM

Every one of these methods sits on the same conceptual chassis as the five-point stencil: replace continuous operators with discrete proxies, respect stability, and march in time. Master the simplest one and the rest will read like dialect.

Connection to Machine Learning

Modern deep-learning research has rediscovered FDM in three flavours, all of which lean on the equation you just implemented:

Common Pitfalls

• Violating CFL by a hair. Set $r = 1.0001$ by accident in a test and the simulation looks fine for 50 steps before exploding. Always pin $r$ to at most 0.99 in production code.
• Forgetting the first-step factor of one-half. Using the standard rule at $n = 0$ drops the scheme to first-order overall and the phase error grows visibly fast.
• Mixing units. $c$, $\Delta x$, and $\Delta t$ must be in consistent units, or your CFL number lies to you. Common SI bug: $c$ in m/s, $\Delta x$ in cm.
• Reflections from open boundaries. Plain Dirichlet turns "outflow" into a wall. Use Mur or PML if you do not want echoes.
• Under-resolved wavelengths. Fewer than ~10 grid points per wavelength makes dispersion catastrophic. Refine the grid or use a higher-order scheme.

Summary

The wave equation, once reserved for clean geometries and clever substitutions, is now in your hands as a thirty-line solver. The ingredients:

1. Discretize on a regular grid; approximate $u_{tt}$ and $u_{xx}$ by central differences.
2. Solve for the future value to get the five-point stencil $u_i^{n+1} = 2u_i^{n} - u_i^{n-1} + r^{2}\Delta_x^{2}u_i^{n}$.
3. Use the half-step Taylor trick for the first iteration so the scheme stays second-order accurate.
4. Respect the CFL bound $r = c\,\Delta t / \Delta x \le 1$. Magic accuracy at $r = 1$; explosion just past it.
5. Pick a boundary condition that matches your physics: fixed for clamped, Neumann for free, Mur for open.
6. Vectorize in NumPy for clarity, port to PyTorch for speed and differentiability, and a Conv1d view drops you straight into modern ML research.

The next sections take this scheme into the wild: musical instruments, electromagnetic waves, and seismic propagation. With the FDM solver in hand, every one of those becomes a small modification, not a new chapter of mathematics.