Solving the 2D Wave Equation with Spectral Methods

The Equation That Models Everything That Vibrates

Sound propagating through air, light through space, seismic waves through the Earth’s crust, ripples on a pond — all of these obey variants of the same fundamental equation: the **wave equation**. It is one of the most important PDEs in physics and engineering, and it serves as a clean, transparent test case for numerical solvers.

In this post I set up a 2D wave equation on a periodic domain, solve it using **Dedalus** — a Python library for spectral PDEs — and animate the result as a matplotlib GIF. The interference patterns that emerge from four initial displacement bumps are both visually striking and mathematically instructive.

The PDE

The 2D wave equation governs the displacement $u(x, y, t)$ of a medium:

where $c > 0$ is the wave speed. The Laplacian $\nabla^2 u$ measures local curvature — regions that are higher than their surroundings are pulled down, and vice versa. This restoring force drives the oscillation.

To integrate equation (1) in time we rewrite it as a **first-order system** by introducing the velocity field $v = \partial_t u$:

This is the standard reduction used by virtually every PDE library. It turns one second-order equation into two coupled first-order equations that a standard time-stepper can handle directly.

Spectral Methods vs Finite Differences

Most introductory PDE courses cover **finite difference methods**: discretise the domain on a grid and approximate derivatives with neighbouring grid values. Spectral methods take a fundamentally different approach — they represent $u$ as a **sum of basis functions** (Fourier modes on a periodic domain, Chebyshev polynomials on a bounded one) and evolve the coefficients of those modes in time.

The payoff is spectral accuracy: for smooth solutions the error decays **exponentially** with the number of modes, rather than as a power of the grid spacing. For the wave equation on a periodic domain, Fourier modes are the natural basis — each mode evolves independently at a frequency dictated by the wavenumber.

**Dedalus** (version 3) is a Python library that implements this approach. You write the PDE in symbolic form, specify the basis and domain, and the library handles the spectral transformations, time-stepping, and parallelism automatically.

Simulation Setup and Code

The domain is $[0, 4\pi] \times [0, 4\pi]$ with periodic boundary conditions. The initial condition consists of four Gaussian displacement bumps placed at the four quadrant centres, with zero initial velocity. The solver uses the fourth-order Runge–Kutta scheme (RK443) with a fixed time step.

import os
os.environ["OMP_NUM_THREADS"] = "1"

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import dedalus.public as d3
import logging

logging.getLogger('dedalus').setLevel(logging.ERROR)

# Parameters
Lx, Ly     = 4 * np.pi, 4 * np.pi
Nx, Ny     = 128, 128
c          = 1.0
timestep   = 0.02
stop_time  = 14.0
n_frames   = 100
timestepper = d3.RK443

# Spectral basis (periodic in x and y)
coords = d3.CartesianCoordinates('x', 'y')
dist   = d3.Distributor(coords, dtype=np.float64)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx))
ybasis = d3.RealFourier(coords['y'], size=Ny, bounds=(0, Ly))

# Fields
u  = dist.Field(name='u',  bases=(xbasis, ybasis))  # displacement
ut = dist.Field(name='ut', bases=(xbasis, ybasis))  # velocity

# PDE: split into first-order system
problem = d3.IVP([u, ut], namespace=locals())
problem.add_equation("dt(u)  = ut")
problem.add_equation("dt(ut) = c**2 * lap(u)")

solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_time

# Initial conditions: four Gaussian bumps
x, y = dist.local_grids(xbasis, ybasis)
sig  = 0.6

def gauss(cx, cy):
    return np.exp(-((x - cx)**2 + (y - cy)**2) / (2 * sig**2))

u['g'] = (gauss(Lx * 0.25, Ly * 0.25) + gauss(Lx * 0.75, Ly * 0.25) +
          gauss(Lx * 0.25, Ly * 0.75) + gauss(Lx * 0.75, Ly * 0.75))
ut['g'] = 0.0

# Run solver and collect frames
frame_times = np.linspace(0, stop_time, n_frames)
frames = []
frame_idx = 0

while solver.proceed and frame_idx < n_frames:
    if solver.sim_time >= frame_times[frame_idx]:
        frames.append(u['g'].copy())
        frame_idx += 1
    solver.step(timestep)

# Animate and save as GIF
vmax = 0.6
fig, ax = plt.subplots(figsize=(7, 6), facecolor='#0d0d0d')
ax.set_facecolor('#0d0d0d')
im = ax.imshow(frames[0].T, origin='lower', cmap='seismic',
               vmin=-vmax, vmax=vmax, extent=[0, Lx, 0, Ly],
               interpolation='bilinear')
plt.colorbar(im, ax=ax).set_label('Displacement u(x,y,t)', color='white')
ax.set_title('2D Wave Equation — Dedalus Spectral Solver', color='white')
time_text = ax.text(0.02, 0.96, '', transform=ax.transAxes, color='white')

def update(i):
    im.set_data(frames[i].T)
    time_text.set_text(f't = {frame_times[i]:.2f}')
    return im, time_text

ani = animation.FuncAnimation(fig, update, frames=len(frames), interval=80, blit=True)
ani.save('wave2d_dedalus.gif', writer='pillow', fps=15, dpi=120)
plt.close()

Interference in Action

The animation below shows the simulation over $t \in [0, 14]$. Each bump launches a circular wave. As the waves expand and the domain is periodic, they wrap around and collide repeatedly — producing the characteristic **interference fringes** visible in the figure. Red regions indicate positive displacement (crest), blue regions negative (trough).

The periodic boundary condition means energy is conserved exactly — no absorption, no dissipation. The total $L^2$ norm of $u$ and $v$ is constant throughout the simulation, which is a useful sanity check for any wave solver.

Where This Equation Appears

The 2D wave equation and its variants arise across a surprisingly wide range of fields:

**Acoustics and room simulation.** Architectural acoustics software solves the wave equation to predict how sound behaves in concert halls, recording studios, and open spaces. Accurate solutions require handling complex geometries, reflections, and absorption boundaries.

**Seismology.** Ground motion from earthquakes is modelled as elastic wave propagation in a heterogeneous medium. Seismic inversion — working backwards from surface measurements to infer subsurface structure — is one of the most computationally intensive applications of PDE solvers in the physical sciences.

**Electromagnetic wave propagation.** Maxwell’s equations in source-free regions reduce to the wave equation for each field component. Antenna design, radar cross-section modelling, and photonic crystal simulation all rely on fast, accurate wave solvers.

**Financial mathematics — the Black-Scholes PDE.** The Black-Scholes equation is a parabolic PDE (diffusion-type), but in certain limiting regimes — particularly near expiry and with high convexity — wave-like behaviour dominates pricing dynamics. The numerical techniques developed for the wave equation inform finite difference schemes for derivative pricing.

**Quantum mechanics.** The time-dependent Schrödinger equation is formally similar to the wave equation. Spectral methods are widely used in computational quantum chemistry and condensed matter physics to evolve wavefunctions accurately over long time horizons.

**Medical imaging.** Ultrasound imaging works by transmitting acoustic pulses and recording reflections. Reconstructing an image from those reflections requires solving an inverse wave problem — a growing area of clinical and computational research.

Takeaway

The wave equation is simple to state but rich in behaviour. Spectral methods — as implemented in Dedalus — solve it to machine precision for smooth data on regular domains, with minimal code. The animated GIF is not just a visualisation: it is a precise numerical solution to a fundamental PDE, computed in seconds on a standard laptop.

Need a reliable numerical implementation for your problem? Get in touch.