Software 5.5.1 (Spherical waves). The residual of the radial wave equation \(u_{tt} = c^2(u_{rr} + \frac {n - 1}ru_r)\) for an outgoing wave in three and in two dimensions:

>>> import sympy as sp
>>> r, t, c = sp.symbols('r t c', positive=True)
>>> F = sp.Function('F')
>>> def radial(u, n):      # u_tt - c^2 (u_rr + (n - 1) u_r / r)
...     return u.diff(t, 2) - c**2 * (u.diff(r, 2) + (n - 1) / r * u.diff(r))
...
>>> sp.simplify(radial(F(r - c*t) / r, 3))         # spherical wave in 3D
0
>>> sp.simplify(radial(F(r - c*t) / sp.sqrt(r), 2))   # none in 2D
-c**2*F(-c*t + r)/(4*r**(5/2))

Reading the output. \(F(r - ct)/r\) is an exact solution in three dimensions for every profile \(F\). Its two-dimensional analogue \(F(r - ct)/\sqrt r\) leaves a residual of order \(r^{-5/2}\): in two dimensions an outgoing wave cannot keep its shape, and is a travelling profile only far from the source. This is the same fact as the wake.

MATLAB/Octave: with the Symbolic Math Toolbox, the same diff expressions.

Software 5.5.2 (Kirchhoff’s formula). Example 5.3.1 in three dimensions, with \(a = 1\), \(R = 5\): the spherical mean in Kirchhoff’s formula estimated from a million random directions, against the exact pulse:

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> w = rng.normal(size=(10**6, 3)); w /= np.linalg.norm(w, axis=1)[:, None]
>>> a, R = 1.0, 5.0              # g = 1 in the unit ball, f = 0, c = 1
>>> x = np.array([R, 0, 0])
>>> def kirchhoff(t):            # u = t * (mean of g over |y - x| = t)
...     return t * np.mean(np.linalg.norm(x + t * w, axis=1) < a)
...
>>> for t in (3.5, 4.5, 5.0, 5.5, 6.5):
...     exact = max(a**2 - (R - t)**2, 0) / (4 * R)
...     print(t, round(kirchhoff(t), 4), round(exact, 4))
...
3.5 0.0 0.0
4.5 0.0371 0.0375
5.0 0.0498 0.05
5.5 0.037 0.0375
6.5 0.0 0.0

Reading the output. The solution is \(0\) before \(t = 4\) and after \(t = 6\), and in between follows \((1 - (5 - t)^2)/20\), with peak \(0.05\) at \(t = 5\), to within the Monte Carlo error, which is \(t\sqrt {p(1 - p)/10^6} \approx 0.0005\) near the peak, where \(p = u/t \approx 0.01\) is the fraction of directions that hit the ball.

MATLAB/Octave: w = randn(1e6,3); w = w./vecnorm(w,2,2); and the same mean.

Software 5.5.3 (Poisson’s formula and the wake). Example 5.3.1 in two dimensions. Writing \(\rho = |y - x|\), Poisson’s formula becomes \(u = \frac 1{2\pi }\int \theta (\rho )\rho \,d\rho /\sqrt {t^2 - \rho ^2}\), with \(\theta (\rho )\) the angle of the circle of radius \(\rho \) about \(x\) that lies in the disc; the substitution \(\rho = t\sin \phi \) removes the square-root singularity:

>>> from scipy.integrate import quad
>>> def poisson(R, t, a=1.0):              # 2D, g = 1 in the unit disc, c = 1
...     def arc(rho):                # angle of circle |y - x| = rho in disc
...         if R == 0: return 2 * np.pi * (rho < a)
...         q = (R**2 + rho**2 - a**2) / (2 * R * rho)
...         return 2 * np.arccos(np.clip(q, -1, 1))
...     lo, hi = max(R - a, 0.0), min(R + a, t)
...     if hi <= lo: return 0.0
...     f = lambda phi: arc(t * np.sin(phi)) * t * np.sin(phi)  # rho=t sin
...     lims = np.arcsin(lo / t), np.arcsin(hi / t)
...     return quad(f, *lims, limit=200)[0] / (2 * np.pi)
...
>>> for t in (0.5, 2.0, 10.0, 100.0):    # centre: t - sqrt(t^2 - a^2)
...     exact = t - np.sqrt(max(t * t - 1, 0))
...     print(t, round(poisson(0.0, t), 6), round(exact, 6))
...
0.5 0.5 0.5
2.0 0.267949 0.267949
10.0 0.050126 0.050126
100.0 0.005 0.005
>>> [round(poisson(5.0, t), 4) for t in (3.5, 5.0, 6.5, 10.0, 30.0)]
[0.0, 0.1736, 0.1271, 0.058, 0.0169]

Reading the output. At the centre the quadrature reproduces \(t - \sqrt {t^2 - 1}\) to six digits, and the value at \(t = 100\), \(0.005\), is the asymptotic \(1/2t\). At distance \(5\) the solution is still \(0.058\) at \(t = 10\) and \(0.0169\) at \(t = 30\), close to \(1/2t = 0.0167\); in three dimensions it would be exactly \(0\) for all \(t > 6\).

MATLAB/Octave: integral with the same integrand.