The 2D Wave Equation: Finite Differences, CFL and Standing Waves

One equation for every wave

Sound in air, ripples on a drum skin, light in a vacuum, a plucked string, a seismic wave through rock — all obey the same second-order partial differential equation. In two dimensions:

∂²u/∂t² = c² ( ∂²u/∂x² + ∂²u/∂y² )  =  c² ∇²u

u(x, y, t) = the displacement (or pressure) field
c          = the propagation speed of the medium

The physical reading is direct. ∇²u is the curvature of the field — the amount by which a point differs from the average of its surroundings. The equation says: a point that sits below the average of its neighbours is accelerated upward, and one above the average is pulled down. It is Hooke's law spread across a continuum, and everything wave-like follows from it.

Discretising it: two lines of finite differences

Put the field on a grid with spacing Δx and sample time at intervals Δt. Replace both second derivatives with the standard central difference, which is second-order accurate:

∂²u/∂t² ≈ ( u_i,j^(n+1) − 2u_i,j^n + u_i,j^(n−1) ) / Δt²

∇²u     ≈ ( u_(i+1),j + u_(i−1),j + u_i,(j+1) + u_i,(j−1)
            − 4·u_i,j ) / Δx²            // the 5-point stencil

Substitute, solve for the only unknown — the field at the next time step — and the entire solver falls out. It is fully explicit: no matrix, no linear system, no iteration. Note that it is a three-level scheme, because the equation is second order in time: you must keep the previous two frames.

const C2 = (c * dt / dx) ** 2;          // Courant number, squared

for (let j = 1; j < h - 1; j++)
  for (let i = 1; i < w - 1; i++) {
    const k = j * w + i;
    const lap = cur[k - 1] + cur[k + 1] + cur[k - w] + cur[k + w]
              - 4 * cur[k];
    next[k] = 2 * cur[k] - prev[k] + C2 * lap;
  }
[prev, cur, next] = [cur, next, prev];  // rotate the three buffers

That is a complete, working 2D wave solver in nine lines. This scheme is often called FDTD — finite-difference time-domain — the same family of methods used for Maxwell's equations in electromagnetics, applied here to the scalar wave equation.

The Courant condition, and the explosion that follows if you ignore it

An explicit scheme is only conditionally stable. Information in the numerical scheme travels exactly one cell per time step (the stencil only reaches its immediate neighbours); information in the physics travels c·Δt. If the physical wave outruns the stencil, the scheme cannot represent it, and the error does not merely grow — it grows without bound. The CFL condition (Courant–Friedrichs–Lewy, 1928) states the requirement precisely. With the Courant number C = c·Δt/Δx, a von Neumann stability analysis of this scheme in d dimensions gives:

C = c·Δt/Δx ≤ 1 / sqrt(d)

1D   C ≤ 1        (and C = 1 is exact — the "magic time step")
2D   C ≤ 1/√2  ≈ 0.7071
3D   C ≤ 1/√3  ≈ 0.5774

Violate it and the simulation does not degrade gracefully: within a few dozen steps the field saturates into a checkerboard of alternating ±∞ and every value becomes NaN. If your wave simulation explodes, the Courant number is the first thing to compute — it is almost always the cause. A useful safety margin is to run at C ≈ 0.5 in 2D rather than at the theoretical ceiling.

The 1D curiosity is worth knowing: at exactly C = 1 the discrete scheme reproduces the exact solution of the continuous equation, because the numerical and physical propagation speeds coincide precisely. That is the only case where a finite-difference scheme is exact, and it is also a warning — it does not generalise to 2D, where C = 1 is comfortably unstable.

Numerical dispersion: the wave that arrives at the wrong speed

Even when it is stable, the scheme is not perfect. Substituting a plane wave into the discrete equations gives a numerical dispersion relation that does not quite match the physical one: on the grid, the propagation speed depends on the wavelength and on the direction of travel. Short wavelengths — those resolved by only a few grid points — travel measurably slower than long ones, and they travel at slightly different speeds along the axes than along the diagonals.

The visible symptom is a sharp pulse that develops a trailing ripple as it propagates, spreading out even though the true equation is non-dispersive, and a circular wavefront that becomes subtly square. The cure is resolution: keep at least 10–20 grid points per wavelength of the shortest wave you care about, and the error becomes negligible. Higher-order stencils (a 9-point Laplacian) reduce the anisotropy for the same grid spacing.

Boundaries: reflect, free, or vanish

What you write into the edge cells defines the physics of the boundary, and there are three useful choices:

Dirichlet   u = 0 at the edge.  A CLAMPED boundary — a drum skin
            nailed to its rim, a string tied down. Waves reflect
            with an inversion (a crest returns as a trough).

Neumann     ∂u/∂n = 0 at the edge, implemented by copying the
            adjacent interior cell into a ghost cell. A FREE
            boundary — an open organ pipe. Waves reflect WITHOUT
            inversion.

Absorbing   The edge should behave as if the domain continued
(Mur, PML)  forever. A first-order Mur condition applies the 1D
            one-way wave equation at the edge; a Perfectly Matched
            Layer surrounds the domain with an artificially lossy
            region that swallows outgoing waves with almost no
            reflection, and is the standard in serious codes.

Doing nothing at all — leaving the border cells at their initial value — is a Dirichlet condition by accident, which is why an unfinished wave solver always shows a box full of reflections. And it is the reflections that produce the standing waves and normal modes of a bounded domain: on a rectangle, a discrete set of frequencies fit the boundary exactly, and the field settles into stationary patterns with fixed nodal lines. Sprinkle sand on a vibrating plate and it collects along those lines — the Chladni figures — which is the same phenomenon this solver reproduces on a grid.

Frequently asked questions

Try it live

Everything above runs in your browser — open Wave Equation and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.