How to solve wave equations on unbounded domains

This post is a practical guide to numerically solving wave equations on unbounded domains. Using standard numerical methods, I demonstrate a technique called hyperboloidal compactification to solve unbounded-domain wave equations directly (no truncation, no artificial boundaries) by mapping the unbounded domain to a finite domain.

Code for the demonstrations uses public libraries and is available as a Google Colab notebook so that you can run everything from a browser.

Why waves?

Our most important senses, sight and sound, rely on waves. Electromagnetic and acoustic waves carry information from an external reality far from the immediate neighborhood of our body to our audiovisual system. Our current understanding of the fundamental building blocks of the Universe is based on waves. Quantum field theory models all elementary particles as wave-like excitations of quantum fields. Most remarkably, we have recently gained the ability to detect gravitational waves, or spacetime waves, from the most energetic phenomena in the Universe happening billions of light years away.

The ability to understand waves by solving the equations that describe them is essential also for practical applications, such as modeling of earthquakes (seismology), ocean waves and tsunamis (oceanography), sonar systems (acoustics), radar scattering off airplanes and antenna design (electromagnetism), among others.

But you probably didn’t need that motivation if you’re reading this.

The mathematical model

The simplest equation that describes a wave is aptly named the wave equation. We’ll consider the case with one space dimension $$ \left(-\partial_t^2 + \partial_x^2\right) u(x,t) = 0. \tag{1} \label{1}$$ The acceleration (or curvature in time), $\partial_t^2$, is balanced by the curvature in space, $\partial_x^2$. This delicate balance leads to solutions that propagate while maintaining their form. In one dimension, there are only two directions a wave can propagate, and we get the general solution $$ u(x,t) = f_+(x-t) + f_-(x+t). $$ The functions $f_\pm$ travel unchanged along the curves $x\pm t$. We call these curves characteristics. The characteristics extend infinitely in both directions.

We solved \eqref{1} explicitly, but more interesting wave equations cannot be solved like this. They include higher dimensions, varying coefficients, potential terms, source functions, nonlinearities, etc. For those cases, we use numerical methods.

But there’s a problem. The wave equation is posed on a spatially unbounded domain, $x\in(-\infty, \infty)$¹. Before we can solve the equation numerically, we need to decide what to do with the unbounded domain.

Why boundaries?

An argument that one finds in some textbooks and papers is that computational resources are finite, so we must truncate the solution domain artificially to some $x\in(-L, L)$ with $L$ suitably large. We then apply boundary conditions to the solution.

This artificial truncation creates two problems:

• boundary treatment, and
• radiation extraction.

The construction of boundary conditions for wave equations has been its own research direction since the 1970s, with thousands of papers written on the problem. The corresponding conditions have names such as absorbing, transparent, non-reflecting, or radiative². Due to the difficulty of the problem, people also constructed damping layers, which are not boundary conditions but layers of artificial dissipation that absorb the wave, such as the perfectly matched layer.

Even when you have a somewhat satisfactory treatment of the boundary, you still don’t have access to the full solution. The full solution, however, is relevant when you’re interested in the radiative fields that are defined at infinity. Whether you want to calculate radar cross-sections of airplanes or gravitational waves from black holes, you need to know the solution at infinity.

Then why don’t we solve the full problem on the full domain? The argument of finite computational resources is clearly wrong. We can map the unbounded domain to a bounded domain. Let’s see how that doesn’t work.

Spatial compactification

There are many options to map the infinite domain $x\in(-\infty,\infty)$ to a finite domain in some other coordinate, $\rho$. A well-known option is the tangens function: $x = \tan \rho$ with $\rho\in[-\pi/2,\pi/2]$³. You can also perform an algebraic compactification, $x=\rho/(1-\rho^2)$, or use the hyperbolic tangent, $x=\mathrm{arctanh}\ \rho$ with $\rho\in[-1,1]$. The basic idea is to map the infinite domain $x$ to a compact domain in $\rho$ and then solve the wave equation numerically on this finite, compact domain.

It sounds obvious and, indeed, has been tried in the 1970s⁴. The problem is that solutions oscillate out to infinity. The spatial mapping converts an unbounded domain problem to an unbounded oscillation problem, and we gain nothing.

Consider a sine-wave solution to the wave equation, $u(x,t) = \sin(x-t)$. If you compactify using the hyperbolic tangent, you get $u(\rho,t) = \sin(\mathrm{arctanh}\ \rho-t)$ which has infinite frequency at the boundaries. We can’t even plot the function, let alone solve for it numerically.

On the level of the wave equation, the propagating wave slows down as it approaches the boundary of the compactified domain, where its propagation speed vanishes. With $\partial_x = (1-\rho^2) \partial_\rho$, we get for Eq. \eqref{1} $$ \left[-\partial_t^2 + (1-\rho^2)^2 \partial_\rho^2 + 2 \rho (1-\rho^2) \partial_\rho \right] u(\rho,t) = 0. $$ The speed of the wave is now $c_\pm = \pm (1-\rho^2)$. Both speeds vanish at both boundaries, which implies that the outgoing characteristics get stuck. The characteristics are $x\pm t = \mathrm{arctanh} \ \rho \pm t$ plotted below.

You can see that the outgoing characteristics pile up near the boundary (blue on the left; green on the right). This slowing down of waves near the boundary makes intuitive sense. The original wave has a constant speed and cannot reach infinity at a finite time. Therefore, the wave gets slower on the mapped domain while its wavelength becomes shorter until we cannot resolve it anymore.

In summary, waves cannot propagate to infinity when we compactify an infinite mathematical domain to a finite numerical domain. They get stuck and blueshifted. We then conclude that spatial compactification is not useful.

Unless, maybe, we’re missing something.

Time transformation

The one-dimensional wave equation \eqref{1} was first written by d’Alembert in 1747, looong before we had any inkling that “space by itself, and time by itself, are doomed to fade away into mere shadows”⁵. Remarkably, the wave equation enforcing balance between curvature in time and curvature in space seems to know about spacetime.

What we miss when we only consider spatial compactification is that the wave equation is a spacetime equation. We must transform time as well.

Consider the outgoing sine wave again. The reason for the infinite frequency problem was that the characteristics blew up near the boundaries. Using $\mathrm{arctanh}\ \rho = \frac{1}{2} (\ln (1+\rho) - \ln(1-\rho))$, we can catch the problematic logarithms near $\rho=\pm 1$ by defining a new time coordinate $\tau$ such that⁶ $$ \tau = t + \frac{1}{2} \ln (1-\rho^2).$$ Now, the characteristics are given by $$ x\pm t = \mp \ln(1\mp \rho) \pm \tau. $$ This behavior is the one we need: only outgoing, no incoming characteristics at both boundaries.

For the wave equation, consider the transformation of the differential operators $$ \partial_t = \partial_\tau, \qquad \partial_x = - \rho \ \partial_\tau\ + (1-\rho^2) \ \partial_\rho. $$ The wave equation becomes

$$ \left(-\partial_\tau^2 - 2 \rho\ \partial_\tau \partial_\rho + (1-\rho^2) \partial_\rho^2 - \partial_\tau - 2 \rho \ \partial_\rho \right) u(\rho,\tau) = 0. \tag{2}\label{2} $$ It looks a bit more complicated than \eqref{1}, but now you can solve the equation numerically on the entire domain! Below is the dynamical solution for Gaussian initial data.

The waves leave the spacetime smoothly through the outer boundaries, which correspond to the two infinities at either side. The solution at these boundaries is the radiative field. No extrapolation is necessary! Combining the spatial compactification with a suitable time transformation solved both the boundary and the radiation problems.

The time transformation that we used is not unique. There are many other options. These time functions are called hyperboloidal. Below is the generalization.

Hyperboloidal compactification

Hyperboloidal compactification consists of a spatial compactification $$ \rho = g(x), $$ together with a suitable time transformation $$ \tau = t + h(x), $$ To see what suitable means here, let’s define $$ H:= \frac{dh}{dx}, \qquad G:= \frac{dg}{dx}. $$ As a compactification, $G$ vanishes at the boundaries. The wave equation becomes after a division by $G$ $$ \left(-\frac{1-H^2}{G} \partial_\tau^2 + 2 H \partial_\tau \partial_\rho + G \partial_\rho^2 + H’ \partial_\tau + G’\partial_\rho \right) u(\rho,\tau) = 0. $$ For this equation to be regular, we require $1-H^2 \sim G$ at the boundaries, which means $|H|\sim 1$ at the leading order. We also want $|H|<1$ away from the boundaries to keep the principal part intact⁷. This is the hyperboloidal compactification.

To solve the hyperboloidal wave equation \eqref{2}, it’s best to write it in symmetric hyperbolic form. We define the auxiliary variables $$ \psi := \partial_\rho u, \qquad \pi := \frac{1-H^2}{G} \partial_\tau u - H \partial_\rho u. $$ The full system reads $$\partial_\tau \pi = \partial_\rho \left[\frac{G}{1-H^2} ( \psi + H \pi ) \right], $$ $$\partial_\tau \psi = \partial_\rho \left[\frac{G}{1-H^2} (H\psi + \pi )\right],$$ $$ \partial_\tau u = \frac{G}{1-H^2} (H \psi + \pi).$$

The system is symmetric. Hyperbolicity follows from the properties of $H$ and $G$. For the case that we discussed above, we have $$ \rho = \tanh x, \qquad \tau = t - \ln(\cosh x) = t + \frac{1}{2} \ln (1-\rho^2). $$ The functions $G$ and $H$ are $$ G = 1-\rho^2, \quad H = -\rho \quad \Rightarrow \quad \frac{G}{1-H^2} = 1. $$ The equation becomes very simple $$\partial_\tau \pi = \partial_\rho \left( \psi - \rho \ \pi \right),$$ $$\partial_\tau \psi = \partial_\rho \left(- \rho\ \psi + \pi \right), \tag{3}\label{3}$$ $$\partial_\tau u = (-\rho\ \psi + \pi).$$

Next, we solve this equation numerically.

Numerical implementation

I picked a simple method that’s relatively accurate: 4th-order finite differences in space and 4th-order Runge-Kutta integration in time. The finite difference operators are one-sided at the boundaries.

def diff(u, h):
    kernel = np.array([-1, 8, 0, -8, 1])
    du = np.convolve(u, kernel, mode="same")
    # One-sided stencils at the boundaries
    du[0] = -25.0 * u[0] + 48.0 * u[1] - 36.0 * u[2] + 16.0 * u[3] - 3.0 * u[4]
    du[1] = -3.0 * u[0] - 10.0 * u[1] + 18.0 * u[2] - 6.0 * u[3] + u[4]
    du[-1] = 25.0 * u[-1] - 48.0 * u[-2] + 36.0 * u[-3] - 16.0 * u[-4] + 3.0 * u[-5]
    du[-2] = 3.0 * u[-1] + 10.0 * u[-2] - 18.0 * u[-3] + 6.0 * u[-4] - u[-5]
    return du / (12.0 * h)

The Runge-Kutta integrator is standard

def rk4(cur):
    nex = np.copy(cur)
    k1 = dt * rhs(nex)
    k2 = dt * rhs(nex + 0.5 * k1)
    k3 = dt * rhs(nex + 0.5 * k2)
    k4 = dt * rhs(nex + k3)
    nex += (k1 + 2.0 * (k2 + k3) + k4) / 6.0
    return nex

Below is the right-hand side calculation. I combine the variables $(u, \psi, \pi)$ into a single array u. The function rhs takes the array u and returns the right-hand side of the system for the symmetric hyperbolic formulation of the hyperboloidal wave equation \eqref{3}.

def rhs(u):
    rhs = np.zeros_like(u)
    hg = np.abs(x[1] - x[0])
    rhs[2, :] = u[0, :] - x * u[1, :]
    rhs[1, :] = diff(rhs[2, :], hg)
    rhs[0, :] = diff(u[1, :] - x * u[0, :], hg)
    return rhs

Then we integrate the equation iteratively in time.

def solve_wave_eq(initial):
    data = np.zeros(((nsteps + 1,) + initial.shape))
    data[0, :, :] = initial
    for i in range(1, nsteps + 1):
        data[i, :, :] = rk4(data[i - 1, :, :])
    print("Finished evolution at t = " + str(dt * nsteps))
    return data

That’s it. The rest is just setting up initial data and parameters. To play around and try different methods, use the Colab notebook as a starting point.

Epilogue

We discussed the simplest example of a wave equation. You may wonder about the wider applicability of this method. It has been applied in higher dimensions with varying coefficients, potential terms, source functions, and nonlinearities. I might write about black-hole perturbations in a future post.

While this technique is relatively well-known in relativity, it’s not used at all in other fields where wave propagation is relevant. It would be fascinating to see this relativistic method applied to a non-relativistic problem, such as the propagation of acoustic waves or the computation of radar cross-sections.