High-order discretization of a stable time-domain integral equation for 3D acoustic scattering

TL;DR

This paper introduces a high-order explicit time-domain integral equation method for 3D acoustic scattering, combining advanced spatial discretization and convolution splines to achieve high accuracy and stability.

Contribution

It presents a novel high-order discretization scheme using Nyström methods and convolution splines for stable, accurate 3D acoustic scattering simulations.

Findings

Achieved 5-9 digit accuracy in simulations

Demonstrated stability up to 8th order

Validated parameters for combined field formulation

Abstract

We develop a high-order, explicit method for acoustic scattering in three space dimensions based on a combined-field time-domain integral equation. The spatial discretization, of Nystr\"om type, uses Gaussian quadrature on panels combined with a special treatment of the weakly singular kernels arising in near-neighbor interactions. In time, a new class of convolution splines is used in a predictor-corrector algorithm. Experiments on a torus and a perturbed torus are used to explore the stability and accuracy of the proposed scheme. This involved around one thousand solver runs, at up to 8th order and up to around 20,000 spatial unknowns, demonstrating 5-9 digits of accuracy. In addition we show that parameters in the combined field formulation, chosen on the basis of analysis for the sphere and other convex scatterers, work well in these cases.

Equations

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u,u(x,0)=\frac{\partial u}{\partial t}(x,0)=0,$$

$$u(x,t)=g(x,t),x\in \Gamma ,t\in [0,T],$$

$$\mathcal{S}\mu (x,t)$$

$$\mathcal{D}\mu (x,t)$$

$$u(x,t)=\mathcal{D}\mu (x,t)+\mathcal{S}\left(a\frac{\partial \mu }{\partial t}+b\mu \right)(x,t).$$

$$\frac{\mu (x,t)}{2}+D\mu (x,t)+S\left(a\frac{\partial \mu }{\partial t}+b\mu \right)(x,t)=g(x,t),x\in \Gamma ,t\in [0,T],$$

$${\mu }_j^k\approx \mu (x_j,t_k).$$

$${\int }_{\Gamma }\frac{f(y)}{4\pi |x_i-y|}d{S}_y$$

$${\int }_{\Gamma }\frac{n_y\cdot (x_i-y)}{4\pi |x_i-y{|}^3}f(y)d{S}_y$$

$${\int }_{\Gamma }\frac{n_y\cdot (x_i-y)}{4\pi |x_i-y{|}^2}f(y)d{S}_y$$

$$f(y_{\ell })=\mu (y_{\ell },t_k-|x_i-y_{\ell }|),\frac{\partial \mu }{\partial t}(y_{\ell },t_k-|x_i-y_{\ell }|),$$

$${\mu }_j(t_k-\tau )=\sum _r{\omega }_r(\tau ){\mu }_j^{k-r},$$

$$\Delta t\ge c_{\text{mbox}\mathrm{iCFL}}h.$$

$${\omega }_r(t_k)={\delta }_{rk}.$$

$$\frac{d^{\ell }{D}_{k+1/2}}{d{t}^{\ell }}({\tau }_k)=\frac{d^{\ell }{P}_k}{d{t}^{\ell }}({\tau }_k),\frac{d^{\ell }{D}_{k+1/2}}{d{t}^{\ell }}({\tau }_{k+1})=\frac{d^{\ell }{P}_{k+1}}{d{t}^{\ell }}({\tau }_{k+1}).$$

$$D_{k+1/2}(\tau )=\sum _{r:{\tau }_r\in S_k\cup S_{k+1}}{y}_r{\omega }_{k+1/2,r}(\tau ),\tau \in ({\tau }_k,{\tau }_{k+1}).$$

$$\omega_{r}(\tau)=\left\{\begin{array}[]{ll}\omega_{k(\tau)+1/2,r}(\tau),&\mbox{if $r$ is such that }\tau_{r}\in S_{k(\tau)}\cup S_{k(\tau)+1},\\ 0,&\mbox{otherwise}.\end{array}\right.$$

$$u(t)+{\int }_0^{1}u(t-s)ds=f(t),t\in [0,T].$$

$$u(t_k)+\sum _{\ell =1}^p{w}_{\ell }u(t_k-{\tau }_{\ell })\approx f(t_k),k=1,\dots ,T/\Delta t,$$

$$u(t_k-\tau )\approx \sum _{r=0}^{\infty }{\omega }_r(\tau )u^{k-r},$$

$$u^k+\sum _{r=0}^{\infty }{W}^r{u}^{k-r}=f(t_k),k=1,\dots ,T/\Delta t,$$

$$W^r=\sum _{\ell =1}^p{w}_{\ell }{\omega }_r({\tau }_{\ell }),r=0,1,\dots$$

$$(1+W^0)u^k=f(t_k)-\sum _{r=1}^{\infty }{W}^r{u}^{k-r}=:{\tilde{f}}^k,k=1,\dots ,T/\Delta t,$$

$$u^k=(1+W^0{)}^{-1}{\tilde{f}}^k.$$

$$\frac{{\mu }_n(t)}{2}+\frac{1}{4}{\int }_0^2{P}_n(1-s^2/2)\left({\mu }_n(t-s)+(2-s){\mu }_n^{\prime }(t-s)\right)ds=g_n(t),$$

$$\mu_{n}^{k}=(1/2+V^{0})^{-1}\biggl{(}g_{n}(t_{k})-\sum_{r=1}^{\infty}V^{r}\mu_{n}^{k-r}\biggr{)}~{},$$

$$V^{r}=\frac{1}{4}\sum_{\ell=1}^{p}w_{\ell}P_{n}(1-\tau_{\ell}^{2}/2)\bigl{(}\omega_{r}(\tau_{\ell})+(2-\tau_{\ell})\omega_{r}^{\prime}(\tau_{\ell})\bigr{)}~{}.$$

$$g_n(t)=10e^{-10(t-2{)}^2},$$

$${\int }_{\Gamma }f(y)d{S}_y={\int }_0^{2\pi }{\int }_0^{2\pi }f(\varphi ,\theta )\Vert z_{\varphi }\times z_{\theta }\Vert d\varphi d\theta \approx \sum _{j=1}^N{w}_{j}f(x_j)$$

$$w_j=\frac{(2\pi {)}^2}{n_{\varphi }{n}_{\theta }}{\eta }_{m_j}{\eta }_{m_j^{\prime }}\Vert z_{\varphi }({\varphi }_j,{\theta }_j)\times z_{\theta }({\varphi }_j,{\theta }_j)\Vert .$$

$${\int }_{\Gamma }\frac{f(y)}{4\pi |x_i-y|}d{S}_y\approx \sum _{j\notin J_{\text{mbox}\mathrm{near},i}}\frac{w_j}{4\pi r_{ij}}f(x_j)+\sum _{\ell =1}^{N_{\text{mbox}\mathrm{aux}}}{\tilde{S}}_{i\ell }f(y_{i\ell }),$$

$${\iint }_T\frac{J(u,v)f(u,v)}{4\pi |x_i-\tilde{z}(u,v)|}dudv={\int }_{{\phi }_0}^{{\phi }_1}{\int }_0^{r(\phi )}\frac{J(r,\phi )f(r,\phi )}{4\pi |x_i-\tilde{z}(r,\phi )|}rdrd\phi$$

Full text

{AMS}

65R20, 65M38

1 Introduction

Problems involving the scattering of waves by obstacles have countless applications in science and technology. For time-harmonic data, numerical methods based on integral equations are well-developed and widely used. Advantages of using an integral equation formulation, rather than the partial differential equation itself, include superior conditioning when second-kind formulations are used, reduced dimensionality, the lack of a need for mesh generation in the volume, and the availability of fast solvers [8, 9, 31]. They also impose outgoing radiation conditions exactly, avoiding the need for artificial boundary conditions on a truncated computational domain when considering exterior problems.

By contrast, numerical methods for time-domain integral equations are not so widely used, and the supporting theory is not nearly as complete. This is due, in large part, to the dependence of the relevant layer potentials on their space-time history, making them somewhat unwieldy in the absence of fast algorithms (see Remark 2 below).

Surprisingly, much of the literature is focused on equations involving the single layer potential (see, for example, the review by Ha-Duong [22] and the recent monograph by Sayas [33]), which lead to first kind equations in either the time or frequency domains. Some exceptions that make use of second kind formulations include [37], in which a quasi-explicit marching scheme is developed for the time domain magnetic field integral equation, and [38], in which high-order accurate Calderon-preconditioned versions of the electric field integral equation are developed.

Here, we focus on the scalar wave equation; our main contributions are (1) the development of a high-order Nyström discretizations for the combined field integral equation proposed in [16] and (2) a numerical study of the stability of explicit marching schemes for the resulting system. In the frequency domain, combined field equations are used to avoid singularities or near-singularities associated with eigenvalues of the interior Helmholtz problem [24] [10, p. 48–49]. For time-domain calculations, it is argued in [16] that these interior resonances will always dominate the long time behavior of the densities, rendering the recovery of the solution in the volume increasingly ill-conditioned. Similar conclusions for electromagnetic problems are reached by Shanker et al in [34], and potential benefits of direct formulations to avoid numerical dispersion appear in the work of Banjai [1].

More precisely, we consider the computation of the scattered wave $u(x,t)$ inducing by an incoming acoustic wave impinging on an obstacle in a uniform medium. The function $u$ satisfies the equation

[TABLE]

for $(x,t)\in\Omega\times[0,T]$ where $\Omega\subset\mathbb{R}^{3}$ is the domain exterior to a compact obstacle with boundary $\Gamma$. In this work, for definiteness, we set unit sound speed $c=1$ and consider the Dirichlet (i.e. sound-soft) problem,

[TABLE]

where $g=-u_{\mbox{\tiny\rm inc}}$ on the surface, and $u_{\mbox{\tiny\rm inc}}$ is a given incident wave. Our methods could also be applied to the Neumann (sound-hard) problem.

For a target point $x\in\Omega$ and target time $t$, the retarded single and double layer potentials applied to a density $\mu(y,t)$, $y\in\Gamma$, are defined by [21, Sec. 10.7] [33, Ch. 1]

[TABLE]

Introducing surface weight functions $a(y)$, $b(y)$, for $y\in\Gamma$ (analogous to the combined field parameter $\eta$ in the time-harmonic case [24]), we represent the scattered wave $u$ in the form

[TABLE]

Taking the limit of (4) as $x\rightarrow\Gamma$, using the jump relation for the double-layer operator [33, Sec. 1.3–1.4], we obtain a linear integral equation for $\mu(x,t)$, $x\in\Gamma$,

[TABLE]

where $D$ is the principal value part of the double layer potential (3) spatially restricted to $\Gamma\times\Gamma$, and $S$ is the single layer potential (2) spatially restricted to $\Gamma\times\Gamma$. When $\Gamma$ is smooth, both $D$ and $S$ have weakly singular kernels with singularities bounded by $1/|x-y|$, as would be the case for the Laplace equation. Note that (5) is of Volterra type in time, that is, at time $t$ the $D$ and $S$ operators involve only the density history $\mu(\cdot,t^{\prime})$ for $t^{\prime}<t$. The first term suggests that it is also of the second kind. However, even for the case $a\equiv b\equiv 0$, it does not fall into the standard Fredholm theory [21, Sec. 10.7], since $D$ is not compact. Rigorous proofs of existence and uniqueness, along with some regularity estimates, for both a single and double layer formulation, may be based on the combination of Laplace transformation in time and subsequent analysis of the parametrized spatial integral equations [22, 33]. Direct proofs of convergence can also be carried out (using the double layer alone) by fixed point iteration [21].

The choice of the parameters $a$ and $b$ and their effect on the long time behavior of $\mu$ is a focus of [16]. There it is shown that:

i.

For convex obstacles, $a(y)\equiv 1$ is a natural choice, since an asymptotic analysis at high frequency indicates that it leads to an optimal cancellation of the leading part of the delay term; see additional discussions in Section 2;

ii.

Any positive $b$ is sufficient to remove the damaging zero-frequency interior Neumann resonance. For $b$ sufficiently large, the long time exponential decay rate of the density $\mu$ will match that of the dominant physical scattering resonances, though for cases other than the sphere an optimal choice of this parameter is unclear. For the sphere one should take $b\equiv R^{-1}$ where $R$ is its radius.

Below, we examine the effect of simple choices for $a$ and $b$ for scattering by a nonconvex obstacle. In particular, numerical experiments in Section 4 compare results for various $a$ and $b$ when $\Gamma$ is the boundary of a torus or a perturbed torus.

Since the boundary is time-invariant, the integral equation (5) is separable in space and time, so it is natural to separate the spatial and temporal discretizations. Here we propose a high-order Nyström-like approximation in space combined with a straightforward time marching procedure. For this, we introduce a set of $N$ surface nodes, $x_{j}\in\Gamma$, $j=1,\ldots,N$ and a regular sequence of time steps $t_{k}=k\Delta t\in[0,T]$, and represent the density $\mu$ by interpolation from its discrete values

[TABLE]

To evolve forward in time, the $k$th time-step consists of applying an explicit, fixed linear rule which gives the vector $\{\mu_{j}^{k}\}_{j=1,\ldots,N}$ in terms of the history $\{\mu_{j}^{k-r}\}_{j=1,\ldots,N}^{r=1,\ldots,n}$ and the current data vector $\{g(x_{j},t_{k})\}_{j=1,\ldots,N}$. Because of the strong Huygens principle, the number of previous times $n$ needed is essentially the diameter of the obstacle divided by $\Delta t$. The contributions of these past values come from approximating the retarded spatial integrals in (5). We account for this history dependence by constructing matrices $S^{h}$, $D^{h}$ and $W^{h}$ such that for each target node $x_{i}$,

[TABLE]

where $h={\mathcal{O}}(N^{-1/2})$ is a measure of the spatial grid spacing, and the convergence order $\gamma$ is determined by the scheme. Here $\{y_{\ell}\}_{\ell=1,\ldots,N^{\prime}}$ is a set of $N^{\prime}>N$ nodes, described in Section 3, comprising all nodes $x_{j}$ far from $x_{i}$, plus a set of auxiliary nodes designed to handle the weakly singular contribution near to $x_{i}$ (see Fig. 4(c)). The point samples $f(y_{\ell})$ in the expressions above involve the density and its time derivative at retarded times, e.g.

[TABLE]

so that the retarded evaluation times do not correspond to points at previously computed time steps, $t_{k-r}$. To handle this, and to approximate the time derivatives, we introduce temporal interpolants at each spatial grid point. As also suggested by Davies and Duncan [11, 12], our interpolants take the form

[TABLE]

where $t_{k}$ is the current time. However, we use different basis functions, $\omega_{r}$. The temporal interpolating functions and some properties of the time-stepping schemes defined in [11, 12] are discussed in Section 2.

The experiments in Sections 3 and 4 verify the high-order convergence of the fully discrete algorithm. In addition, the stability properties of the time marching scheme, in its explicit predictor-corrector form, are investigated, revealing a rather surprising “inverse CFL” constraint of the form

[TABLE]

Such conditions have been noted before, but typically in the context of finite volume methods [4, 7, 25].

We note that most authors use fully implicit time marching schemes, although there are exceptions [37]. While the linear systems to be solved involve only local interactions, a direct solver on a surface, even with optimal ordering, is likely to require ${\mathcal{O}}(N^{3/2})$ work to factor, with subsequent solves required at each time step requiring ${\mathcal{O}}(N\ln{N})$ flops. Clearly, an explicit method (including predictor-corrector iterations) is cheaper, as no matrix factorizations are necessary and the cost per time step is ${\mathcal{O}}(N)$. Finally, in Section 5 we summarize our results, and point to future enhancements and generalizations of our method.

Remark 1** (Software).**

*An open source MATLAB (with some Fortran90) implementation of the methods of this paper is freely available from the repository https://github.com/ahbarnett/BIE3D, with code for most of the figures from this paper in the timedomainwaveeqn/paper directory.

Remark 2** (Fast algorithms).**

*A critical bottleneck in using retarded layer potentials for the solution of the wave equation is that they require ${\mathcal{O}}(N^{2})$ memory and ${\mathcal{O}}(N^{2})$ work per timestep. Fortunately, over the last two decades fast algorithms have been developed which reduce both of these costs to ${\mathcal{O}}(N\log N)$. These are described, for example, in [3, 17, 9, 28, 35, 41] and can be used to accelerate most Galerkin or Nyström discretization schemes, including the method developed here. In the present work, we make use of direct summation methods for the sake of simplicity.

2 Temporal discretization

Two commonly used approaches to time-stepping time-domain integral equations are Galerkin methods, discussed extensively in Ha-Duong’s review [22], and convolution quadrature [27], discussed extensively in Sayas’ monograph [33]. Although provably stable with exact integration, Galerkin methods have been found to be sensitive to quadrature errors arising from the need to compute integrals in cut elements when discontinuous polynomial bases in time are employed. As a result, various smooth nonstandard basis functions have been proposed [39, 32]. Convolution quadrature methods, on the other hand, have been found to be more robust. Their construction, based on combining standard time-marching schemes for ordinary differential equations with representations of convolution in the Laplace domain, leads to methods which do not respect the strong Huygens principle obeyed by the layer potentials (2)–(3). That is, the solution updates involve the entire time history of the potentials. As such, special methods must be introduced to alleviate the storage and computation costs for long time computations; see, e.g., [1].

Most relevant to our approach is the convolution spline method of Davies and Duncan [11, 12]. They present a simple reinterpretation of convolution quadrature as a temporal approximation (8), and suggest taking $\omega_{r}(t)$ to be smooth, compactly-supported splines (both standard B-splines and more exotic bases), thus restoring the finite time history of the layer potential operators. However, their approximations (8) are quasi-interpolatory; that is, the basis functions $\omega_{r}$ do not satisfy

[TABLE]

As a result, their proposed methods are limited to second order accuracy, even if B-splines of high degree are used. Higher order accuracy can be achieved using “marching on in time” (MOT) schemes, which are more closely related to the approach presented here (see [38] and references therein). We propose the use of smooth piecewise polynomial bases which are designed to satisfy (10). The functions, which we term “difference splines” (or D-splines for short) [2] are defined as follows.

i.

Let $q$ be an integer and $\{\tau_{k}\}_{k=0,1,\dots}$ be a set of nodes, with $\tau_{0}=0$, and $\tau_{k+1}>\tau_{k}$ for each $k$. These nodes need not be the $t_{k}$ from the previous section. Let $S_{k}$ be the set of $2q+1$ nearest nodes to be used in the difference stencil for $\tau_{k}$. There are two cases:

a.

(Interior node, $k\geq q$): $S_{k}:=\{\tau_{k-q},\ldots,\tau_{k+q}\}$.

b.

(Boundary node, $k<q$): $S_{k}:=\{\tau_{0},\ldots,\tau_{2q}\}$.

ii.

Let $P_{k}(\tau)$ be the degree $2q$ Lagrange interpolating polynomial defined by data, $\{y_{r}\}$, on the stencil $S_{k}$. That is, $P_{k}(\tau)$ is the unique polynomial of degree $2q$ satisfying $P_{k}(\tau_{r})=y_{r}$ for all $r$ with $\tau_{r}\in S_{k}$.

iii.

On the interval $(\tau_{k},\tau_{k+1})$ the D-spline interpolant of the data $\{y_{r}\}$, $D_{k+1/2}(\tau)$, is defined as the degree $4q+1$ Hermite interpolant of $P_{k}(\tau)$ and $P_{k+1}(\tau)$. Precisely, $D_{k+1/2}(\tau)$ is the unique polynomial of degree $4q+1$ satisfying for $\ell=0,\ldots,2q$:

[TABLE]

As the interpolation operators are linear, and $D_{k+1/2}(\tau)$ depends only on the stencil data, there exist degree-$(4q+1)$ polynomials $\omega_{k+1/2,r}(\tau)$, $\tau_{r}\in S_{k}\cup S_{k+1}$, such that

[TABLE]

iv.

Let $k(\tau)$ be the index of the rightmost node not larger than $\tau$, i.e. such that $\tau_{k(\tau)}\leq\tau<\tau_{k(\tau)+1}$. Then the D-spline basis function $\omega_{r}(\tau)$ associated with node $r$ is defined as the piecewise degree-$(4q+1)$ polynomial

[TABLE]

By construction $\omega_{r}(\tau)\in C^{2q}$, and, for uniformly spaced nodes in time, they are translates of a single simple basis function away from $\tau=0$. In addition, the D-spline interpolant reproduces polynomials of degree $2q$, which by standard results implies accuracy of order $2q+1$ for function values and $2q$ for derivatives.111Here we will identify the interpolants by the polynomial degrees they exactly reproduce; that is the degree-$2q$ D-spline is the piecewise degree $4q+1$ function described above.

In the uniform grid case, we plot the interior functions $\omega_{r}$ for $2q=2$ to $6$, along with boundary functions for $2q=4$, in Figure 1. A Fortran90 implementation (with MATLAB interface) of the above D-splines on regular grids is available in the timedomainwaveeqn/timeinterp directory; see Remark 1.

2.1 Application to a simple Volterra equation

Before attacking the full spatio-temporal problem, we illustrate the use of the above temporal basis functions to create a collocation time-stepping scheme for a simple 2nd-kind Volterra integral equation with “top-hat” kernel,

[TABLE]

We assume $f$ is smooth, and $u(t)$ is known for $t\leq 0$ and smooth for all $t\leq T$. As in Section 1 let $t_{k}=k\Delta t$, $k\in\mathbb{Z}$ be the time grid. Let $\{\tau_{\ell}\}_{\ell=1}^{p}$ be the nodes and $\{w_{\ell}\}_{\ell=1}^{p}$ be the corresponding weights of a $p$-node Gauss–Legendre quadrature rule on $[0,1]$. Enforcing (14) at $t=t_{k}$, and inserting the quadrature,

[TABLE]

holds to high accuracy because the integrand is smooth. Let the unknowns be $u^{k}:=u(t_{k})$, for $k=1,\ldots,T/{\Delta t}$. Let $\omega_{r}(\tau)$ be the degree-$2q$ D-spline basis, as defined in the previous section, for the regular grid of spacing ${\Delta t}$, which one may think of as stepping backwards in time from the current time $t_{k}$. The resulting interpolant is

[TABLE]

where, abusing notation slightly, we take $u^{k}$ for the pre-history $k\leq 0$ to be populated with the known solution $u(t_{k})$. Substituting this interpolant into (15) defines a lower-triangular Toeplitz linear system

[TABLE]

with weights computed by

[TABLE]

Note that in (17) the upper limit for $r$ can be replaced by $n\geq 1/{\Delta t}+q$, which is sufficiently large to capture all history dependence in (14).

The system (17) is naturally best solved sequentially for unknowns $k=1,2,\dots$, i.e. by time-stepping, since each row can be written

[TABLE]

where the “right-hand side” $\tilde{f}^{k}$ is explicitly given in terms of data and previous values. Of course, in this simple example, the “solve” for $u^{k}$ at each time step is trivial:

[TABLE]

Figure 2 shows that this scheme achieves an empirical order $2q+2$. In this experiment, $f$ is constructed numerically (via an accurate quadrature) such that the solution $u$ is a known function. Note that, although in this test case $u(t)$ is not strictly zero for $t\leq 0$, it is of order $10^{-16}$ or less, so that initialization by zero data for the pre-history of $u$ induces negligible error.

The figure also shows stability for arbitrarily small ${\Delta t}$, even with fixed number $p$ of quadrature nodes. The mean mesh spacing $1/p$ is the closest analog to $h$ in the full spatio-temporal scheme of Section 3. Thus, in contrast to the full scheme, and also to the modal problem for the sphere in Section 2.2, we observe no inverse CFL condition. In fact, for most of the ${\Delta t}$ tested, the quadrature (18) is completely unresolved for each basis function $\omega_{r}$, yet the scheme is accurate—since $u$ itself is resolved—and stable.

Remark 3**.**

*Recall that in the context where there are $N$ spatial variables, the analog of the solution (20) to (19) requires the solution of an $N\times N$ linear system at each time step $k$, i.e. an “implicit” time step, which can incur a large cost. As mentioned above, for the full problem we will focus on predictor-corrector schemes, which replace this by a fixed number of explicit steps and require only ${\mathcal{O}}(N)$ effort. *

2.2 Modal problem on the unit sphere

As a second example of a simple Volterra-like scalar problem, but one which is more directly related to the full problem we aim to solve, we consider our combined field equation with $a=b=1$ restricted to the amplitudes, $\mu_{n}(t)$, of a spherical harmonic expansion of the full solution on the unit sphere. As shown in [16], the integral equation to be solved is

[TABLE]

where $P_{n}$ denotes the degree-$n$ Legendre polynomial. The key difference between this example and the one considered above is the presence of the time derivative of the density. We believe this term leads to the inverse-CFL constraint observed for the full model, and want to verify that it appears in this simpler case where no spatial integrations are required. We note that in [16] the time derivative term is removed via integration by parts. Such a procedure would remove the difficulty, but it is unclear how it could be accomplished for the full problem with a general scatterer.

Our goal will simply be to determine the stability limits, if any, on the choice of time step. We discretize exactly as above: the integral is approximated by a $p$-point Gauss rule and the D-spline interpolant of $\mu_{n}$ as well as its derivative is evaluated at the quadrature nodes to produce a scheme

[TABLE]

with weights

[TABLE]

Choosing $g_{n}$ to be, for every $n$ that we test,

[TABLE]

we solve to $T=100$, for $p$ varying between $64$ and $1920$ in increments of $64$ and for all $n$ from [math] to $p/2$. The method is deemed unstable for a given time step if the maximum value of the solution for any $n$ grows beyond the expected value by a factor of roughly $1.1$. If we take the average mesh width to be $h=2/p$ and plot $p$ versus the minimum stable CFL number, ${\Delta t}_{{\mbox{\tiny\rm min}}}/h$, for $2q=2,4,6$ we obtain the results shown in Figure 3. We see for $h$ small the minimum CFL number is approximately $2.6$ for $2q=4,6$ and $2.5$ for $2q=2$. We note that for $2q=8$ (a $10$th-order method) the stability region shrinks considerably since a maximum limit of about $3$ also appears; that is, the method is stable in this test only for $2.5h<{\Delta t}<3h$. Thus we restrict ourselves to methods with $2q$ in the range 2 to 6 (orders $4$ to $8$), where the time step, once above its minimum stable value, appears to be controlled by accuracy and not stability, although we have no formal proof of such a stability result.

3 Spatial discretization and full scheme

First we present then test a quadrature scheme to apply the retarded single- and double-layer potentials (2) and (3). Finally we present the full interpolation scheme and Volterra time-stepping for (5).

3.1 Retarded layer potentials for off-surface targets

The case of an exterior evaluation target $x\in\Omega$, far from $\Gamma$, is simple: for densities $\mu(y,t)$ that are smooth with respect to both $x\in\Gamma$ and time $t$, the integrands are also smooth, and a standard quadrature scheme using the density interpolatory nodes $x_{j}\in\Gamma$ will be accurate. In this work we restrict ourselves to smooth torus-like surfaces. Such surfaces can be parameterized by an infinitely differentiable, doubly $2\pi$-periodic function $x=z(\phi,\theta)$, where $z:[0,2\pi)^{2}\to\mathbb{R}^{3}$. We now describe a simple composite (i.e. panel-based) rule to generate the nodes $x_{j}$ and weights $w_{j}$, which we use in later tests, such that for any smooth $f:\Gamma\to\mathbb{R}$,

[TABLE]

holds to high-order accuracy. Here $z_{\phi}$ and $z_{\theta}$ are the partials of $z$. We cover the parameter space $[0,2\pi)^{2}$ with a uniform ${n_{\phi}}$-by-${n_{\theta}}$ grid of rectangular patches (“quads”). Each patch is covered by a tensor product grid comprising a $p$-node Gauss–Legendre rule in each of the two directions. Let $x_{j}$, $j=1,\dots,N$, where $N={n_{\phi}}{n_{\theta}}p^{2}$, be the images of these parameter nodes under the map $z$. Figs. 4(a) and (b) show two surfaces with some of their resulting nodes. The corresponding weights $w_{j}$ are found as follows. Let $\eta_{m}$, $m=1,\dots,p$, be Gauss–Legendre weights on the interval $[0,1]$. For surface node $j$, let $\phi_{j}$ and $\theta_{j}$ be its parameter values, and $m_{j},m^{\prime}_{j}\in\{1,p\}$ be its indices in the two directions within the appropriate panel. Then

[TABLE]

The expected convergence order for (24) is ${\mathcal{O}}(h^{2p})$, where $h={\mathcal{O}}(N^{-1/2})$ is the resolution (combining [13, (2.7.12)] with a theorem on composite rules [13, Sec. 2.4]).

3.2 Retarded layer potentials for on-surface targets

When the target $x$ is on $\Gamma$, as needed in the integral equation (5), then the integrand has the following type of weak singularity. If one smoothly parametrizes $\Gamma$ via local polar coordinates $(r,\varphi)$ centered at the target point $x$, for each $\varphi$ (i.e. radial line) the integrand is $1/r$ times a smooth function of $r$. Surprisingly, this is the same form as for the 3D elliptic BVP case: the kernels in (2) and (3) are identical to (or in the 2nd term of (3), less singular than) the Laplace kernels, and although the retardation introduces a conical singularity into the density, this does not change the singularity of their product.

For the elliptic case there exist many high-order Nyström quadrature schemes in 3D. For surfaces diffeomorphic to the sphere, a global spherical harmonic basis [40] [10, Sec. 3.6], or spherical grid rotation [18, 19] achieves spectral accuracy. For more general smooth surfaces, Bruno–Kunyansky [6, 42] use a smooth partition of unity to isolate the singular near-target contribution. To handle the latter they exploit the fact that the polar metric $rdrd\varphi$ cancels the $1/r$ singularity in the integrand, so integrate using auxiliary quadrature nodes on a polar grid centered at the target. For general high-order triangulations, including those with high aspect ratio, Bremer–Gimbutas [5] developed generalized Gaussian quadratures that again use auxiliary nodes. Note that in these elliptic schemes, the integral kernel is directly evaluated at each of the auxiliary nodes, but the density must be spatially interpolated from its values at the nodes $x_{j}$.

Building on the above, we present a simple high-order accurate scheme to evaluate the retarded layer potentials (2) and (3) on surfaces discretized by a structured rectangular grid of quad panels of the type described above in Sec. 3.1. Let $x_{i}$ be a target node, in panel $k$. This panel and its eight neighbors form a $3\times 3$ block of “near” panels, containing node indices $j\in J_{{\mbox{\tiny\rm near}},i}$ (the black and green nodes in Fig. 4(a)), leaving $({n_{\phi}}{n_{\theta}}-9)$ “far” panels. With respect to each of the latter panels, the singularity of the kernel is distant, so their native quadrature nodes $x_{j}$ and weights (25) may be accurately used, as in the previous section. This explains the first term in our approximation

[TABLE]

where $r_{ij}:=|x_{i}-x_{j}|$; the second term is explained below. Note that (26) is of the form (6), with $N^{\prime}=({n_{\phi}}{n_{\theta}}-9)p^{2}+N_{\mbox{\tiny\rm aux}}$. The expressions for the other two kernels are analogous.

The second term in (26) accounts for the contribution of the near $3\times 3$ panel block via a new target-specific set of $N_{\mbox{\tiny\rm aux}}$ auxiliary nodes. For this it is convenient to switch to the standard parametrization $(u,v)\in[-3,3]^{2}$ for this near block, with the target panel $k$ preimage being $[-1,1]^{2}$; see Fig. 4(c). Precisely, there is a simple affine map to the global parameters $\phi=\pi(2i_{k}+u+1)/{n_{\phi}}$ and $\theta=\pi(2i^{\prime}_{k}+v+1)/{n_{\theta}}$, where $i_{k}$ and $i^{\prime}_{k}$ are the toroidal and poloidal indices (i.e. integer coordinates) of panel $k$ within the panel grid. Then denote by $\tilde{z}$ the ($k$-dependent) map from $(u,v)$ to $\mathbb{R}^{3}$, which is the above affine map composed with the map $z$. The auxiliary nodes comprise four grids, each of which integrates over one of the four triangles connecting the block walls to the preimage of the target. Together the four triangles cover $[-3,3]^{2}$; see Fig. 4(c).

Specifically, let the polar coordinates $(r,\varphi)$ be centered at the target preimage in $(u,v)$, and consider one of the four triangles, $T$, that lies in the angle range $[\varphi_{0},\varphi_{1}]$ and whose far edge is given by $r(\varphi)$ in polar coordinates. For the single-layer (2), the integral of a retarded density $f$ as in (7) over the image of this triangle on $\Gamma$ is

[TABLE]

where $J(u,v)=\|\tilde{z}_{u}(u,v)\times\tilde{z}_{v}(u,v)\|$ is the Jacobian of the map to the surface, and we abuse notation slightly so that $J(r,\varphi)$ means $J(u(r,\varphi),v(r,\varphi))$, etc. Let $\rho_{m}$ and $\eta_{m}$, $m=1,\dots,n_{r}$, be respectively the nodes and weights of a $n_{r}$-point Gauss–Legendre rule on $[0,1]$. Let $\varphi_{n}$ and $\xi_{n}$, $n=1,\dots,n_{\varphi}$, be the nodes and weights of a $n_{\varphi}$-point Gauss–Legendre rule on $[\varphi_{0},\varphi_{1}]$. Then (27) is approximated by

[TABLE]

where $J_{mn}:=J(\rho_{m}r(\varphi_{n}),\varphi_{n})$, etc, indicates the value at the auxiliary node indexed by $n$ and $m$. High-order convergence is expected for (28) since, although $f$ has a conical singularity at the polar origin, along constant-$\varphi$ rays the integrand times $r$ is smooth. Summing the four triangles, there are $N_{\mbox{\tiny\rm aux}}=4n_{r}n_{\varphi}$ auxiliary nodes for each target point. The weights $\tilde{S}_{i\ell}$ in (26) may be read off by associating each $\ell$ with a term $nm$ in (28), and taking all factors except the density sample $f_{nm}=f(y_{i\ell})$. The weights $D_{i\ell}^{h}$ and $W_{i\ell}^{h}$ in (6) are found in an analogous way.

Remark 4** (order of convergence).**

*The convergence of this on-surface scheme is subtle, due to its split into far and near source panels. Consider the error due to the $p\times p$-node smooth rule on one of the nearest “far” panels (i.e. just outside of the $3\times 3$ block): as $h\to 0$, its integrand has a singularity that remains at a fixed distance relative to the panel size, so its error is expected to drop only in proportion to the panel area times the typical integrand. Thus, for any fixed $p$, the formal order $\gamma$ is low ($\gamma<2$). Yet, this matters little in practice because its prefactor is expected to be exponentially small in $p$. This follows by analogy with 1D $p$-node Gauss quadrature on $[-1,1]$, for which the error is of order $\varepsilon=\rho^{-2p}$, where $\rho$ is the size parameter of a Bernstein ellipse in which the integrand is analytic [36, Thm. 19.3]. Since the nearest singularity is at $\pm 3$ scaled to the standard panel $[-1,1]$, solving $3=(\rho+\rho^{-1})/2$ gives $\rho\approx 5.8$, hence for $p=4$, $\varepsilon<10^{-6}$, and for $p=8$, $\varepsilon<10^{-12}$. Combining with the accuracy of the rest of the far panels, and assuming $n_{r}>p$, one can summarize the expected error as ${\mathcal{O}}(\varepsilon h^{q}+h^{2p})$. We postpone a rigorous analysis for future work.

Remark 5**.**

*Although similar to that of [5], our method is simpler and somewhat more efficient, since we exploit the structured nature of the panel grid to cover both self-interaction and neighboring panel interactions with a single auxiliary node rule. This is only possible because, for our class of surfaces, the charts from each panel extend to their neighbors in a known, smooth fashion. We also do not attempt to handle panels of aspect ratios much larger than 2. *

3.3 Validation of retarded layer potential evaluation

Before presenting the full spatio-temporal scheme, we pause to describe some numerical tests of the above spatial quadratures for retarded potentials.

Surfaces used for tests. In this work we use a simple class of smooth surfaces $\Gamma$ diffeomorphic to a torus. Given $H(\phi,\theta)$, a doubly $2\pi$-periodic smooth height modulation function, then the map $x=z(\phi,\theta)$ is given in Cartesian coordinates by

[TABLE]

Thus, the major radius is 1. The plain torus of Fig. 4(a) is given by the constant $H\equiv 0.5$, and is close to being the most benign surface of this topology. The “cruller” of Fig. 4(b) has $H(\phi,\theta)=0.5+0.1\cos(5\phi+3\theta)$, and is more challenging due to its higher curvature, and ridges which do not align with the parameter directions. We will not list the analytic partial derivatives of $z$ here, although of course they are needed for Jacobian computations. Note that both objects have a maximum diameter of approximately 3.

Green’s representation formula. We test the convergence of the above quadrature schemes by numerically checking Green’s representation formula (Kirchhoff’s formula) for the wave equation [33, (1.17)], which for convenience we now state. Let $u(x,t)$ satisfy (1) in the closure of the exterior domain $\Omega$, for all $t$, and let $u^{+}$ and $u_{n}^{+}$ indicate respectively the value and normal derivative on $\Gamma$, then

[TABLE]

(Note that, since $u$ is a solution for all time $t$, there is no need for initial conditions as in [33, Sec. 1.4].) We use this to test both the off-surface (Sec. 3.1) and on-surface (Sec. 3.2) quadrature schemes. We use for $u$ an exterior solution to the wave equation generated by a generic point source $x_{0}=(0.9,-0.2,0.1)$ inside, but far from, the surface $\Gamma$, namely

[TABLE]

with a source signal $T(t)=\cos 5t$. Fig. 5(a) shows a snapshot of the resulting retarded double-layer density for an on-surface target $x$. In this section we show results only for the cruller, omitting the more accurate and predictable results for the plain torus.

Off-surface test. Figs. 5(b) and (c) shows (with dots) the convergence of (30) for an exterior target $x=(1.3,0.1,0.8)$, which is a generic point a distance 0.36 from $\Gamma$ (see subfigure (a)). We use the smooth panel scheme of Sec. 3.1. For both shapes, ${n_{\theta}}=3{n_{\phi}}/2$ gives panels with low aspect ratios. Recalling that each quadrature panel has $p\times p$ nodes, (b) shows $p=4$ (for $N$ ranging from 384 to 11616), while (c) shows $p=8$ (for $N$ from 1586 to 24567). The horizontal axis shows the mean spatial node spacing, or resolution, which we define as

[TABLE]

Although there is variation, convergence is asymptotically consistent with the expected order $2p$. The variation is absent and the errors much smaller for the plain torus (not shown). For $p=8$, 12-digit accuracy is reached at the highest $N$.

On-surface test. A surface target $x\in\Gamma$ is shown in Fig. 5(a), and the singular auxiliary node scheme of Sec. 3.2 used. This target was in the corner of a panel, close to parameters $\phi=\theta=0$, and thus involved worst-case auxiliary triangle aspect ratios. For panel aspect ratios up to around 2, we found for the auxiliary quadrature that $n_{\varphi}=2n_{r}$ was adequate, so we fixed this ratio. We explore three choices of $n_{r}$ for each choice of $p$. For $p=4$, Fig. 5(b) shows that, for $h\leq 0.12$ (${n_{\phi}}\geq 12$), errors due to the singular scheme are negligible relative to the overall error for $n_{r}=12$. For $p=8$, a higher $n_{r}$ is needed: $n_{r}=20$ causes negligible errors for $h\leq 0.1$ (${n_{\phi}}\geq 9$). In general, for all but the largest $h$ available, the errors of the singular scheme are negligible for the choice $n_{r}=2p+4$.

Remark 6** (Choice of singular scheme order $n_{r}$).**

*Based on such tests, we fix $n_{r}=2p$ for the torus and $n_{r}=2p+4$ for the cruller. Thus the auxiliary scheme is of much higher order than the underlying panels. This allows us to explore larger $h$ without loss of accuracy or stability. At smaller $h$ one could reduce $n_{r}$ and hence the effort for the auxiliary scheme.

With $n_{r}$ converged as above, the $h$-convergence is as expected from Remark 4: errors drop with order roughly $2p$ (dashed lines), until they saturate to a low-order convergence at around $10^{-7}$ (for $p=4$) or $10^{-10}$ (for $p=8$).

3.4 Interpolation and explicit time-stepping

We now describe how the above spatial quadrature for retarded potentials is combined with the time-stepping of Sec. 2 to solve the time-dependent BIE (5). Enforcing the BIE on the time grid $t=t_{k}$ (as in Sec. 2.1), and on the spatial nodes $x_{i}$, gives

[TABLE]

Now, fixing constant $a$ and $b$, we approximate the action of the retarded integral operators on the density interpolated from the space-time data $\mu_{j}^{k}$ by a set of $N\times N$ matrices $A^{r}$, thus

[TABLE]

where, analogously to Sec. 2.1, $n\geq\mbox{diam}(\Gamma)/{\Delta t}+q$ is large enough to capture all history dependence via Huygens’ principle plus the support of the D-splines.

Spatio-temporal interpolation. Each above matrix $A_{r}$ is filled as follows. For simplicity, consider only the plain single-layer ($S\mu$) contribution in (33). Applying spatial quadrature (26), then time interpolation (8), gives

[TABLE]

The first term ($j$ far from $i$) is already in the form (34), so the contribution to $A^{r}_{ij}$ in this case is $\frac{w_{j}}{4\pi r_{ij}}\omega_{r}(r_{ij})$. We now use spatial interpolation on each time slice $t_{k-r}$ to turn the auxiliary term also into a weighted sum over $\mu_{j}^{k-r}$. Let $L_{j}(u,v)$, for $j\in J_{{\mbox{\tiny\rm near}},i}$, be a set of $9p^{2}$ basis functions that interpolate over the preimage of the near $3\times 3$ patch. I.e., for any smooth function $f$,

[TABLE]

holds to high order, where $f_{j}$ are the values at the interpolatory nodes (preimages of $x_{j}$). Let $(u_{i\ell},v_{i\ell})$ be the standard parameters of the $\ell$th auxiliary node for the target node $i$. Spatial interpolation of $\mu(\cdot,t_{k-r})$ then approximates the 2nd term of (35) by

[TABLE]

which is now also of the form (34). Proceeding as above with the other two terms of (33) (and recalling that $D$ has two terms (3)), finally gives the formula

[TABLE]

For the basis $L_{j}$ we use the $p\times p$ product Lagrange basis for whichever of the nine panels $j$ lies in, and zero elsewhere. Precisely, let $k_{j}\in\{1,\dots,9\}$ be the panel in which node $j$ lies, let $i_{j}$ and $i^{\prime}_{j}$ be its two index coordinates within that panel, and let $(u_{0}^{k},v_{0}^{k})$ be the parameter offset of panel $k$ relative to the target panel ($u_{0}^{k}$ and $v_{0}^{k}$ are either $-2$, [math], or $2$). Then,

[TABLE]

where $l_{i}(x)$ are the usual 1D Lagrange polynomials for the $p$ Legendre nodes on $[-1,1]$. This form aids bookkeeping since it decouples all interactions into independent panel-panel pairs, each of which is either near or far. The $r$ different $p^{2}$-by-$p^{2}$ blocks of $A^{r}_{ij}$ given by all nodes $j$ in a single source panel interacting with all nodes $i$ in a single target panel may be filled together. In practice the $N_{\mbox{\tiny\rm aux}}\times p^{2}\times p^{2}$ nonzero entries of $L_{j}(u_{i\ell},v_{i\ell})$ are precomputed once and for all, then for each target $j$ the sum over $\ell$ in (37) is performed as an efficient matrix-matrix multiplication (GEMM). Note that the near-panel interpolation error is expected to be ${\mathcal{O}}(h^{p})$.

Remark 7**.**

*One might be tempted to interpolate with respect to $y$ the retarded densities such as $\mu(y,t-|x_{i}-y|)$; however, this would fail to be high-order accurate due to the conical singularity around the target $y=x_{i}$. Instead one must interpolate in both space and time, as above, since as a function of space and time $\mu(y,t)$ is smooth. *

Fig. 6 shows the sparsity patterns of a selection of the resulting $A^{r}$ matrices. As expected by Huygens’ principle, $A^{0}$ is concentrated around the diagonal, but as the time delay $r$ increases, the shell of influence spreads across the panels, departing at the most distant (furthest off-diagonal) panels.

Predictor-corrector scheme. Using the notation ${\bm{\mu}}^{k}:=\{\mu_{j}^{k}\}_{j=1}^{N}$ for the vector of densities at time step $k$, we rewrite (34) as a linear system to be solved at each time step,

[TABLE]

note the similarity to (19). As discussed in the introduction, rather than an implicit solve of (38) for each time step, we prefer the following explicit scheme which achieves the same order.

Consider the $k$th time step. Firstly a “predictor” ${\bm{\mu}}^{(0)}$ is generated via a fixed order-$m$ extrapolation rule in time applied to the density vector,

[TABLE]

where we choose $m=2q$ to match the D-spline order. Here the $c_{r}$ are simply the values of the Lagrange polynomials associated with the time nodes $t_{k-r}$ evaluated at $t_{k}$. Then $\tilde{{\mathbf{g}}}^{k}$ is evaluated according to the right-hand side of (38); this is the most expensive task. Finally $n_{c}$ “corrector” steps are performed on ${\bm{\mu}}$, each of which is a Jacobi iteration with shift $d$ ,as follows. The system matrix is split into the diagonal matrix $B$ and matrix $\tilde{A}^{0}$, defined by

[TABLE]

The linear system (38) to be solved is then $(\tilde{A}^{0}+B){\bm{\mu}}^{k}=\tilde{{\mathbf{g}}}^{k}$. Writing the $\alpha$th iterate for this solution as ${\bm{\mu}}^{(\alpha)}=\{\mu^{(\alpha)}_{j}\}_{j=1}^{N}$, and initializing with (39), the Jacobi iteration is

[TABLE]

Once may interpret each iteration as decrementing $\mu_{j}$ by the $j$th component of its residual divided by $B_{jj}$. This iteration, if it converges, converges to the exact (implicit) solution. However, to make an explicit scheme we fix $n_{c}$, independent of $N$, so that the approximate solution to the time step is ${\bm{\mu}}^{k}={\bm{\mu}}^{(n_{c})}$.

This completes the description of the entire scheme for evolving the density. The cost per time-step is ${\mathcal{O}}(N^{2})$, and thus the total ${\mathcal{O}}(N^{2}\cdot T/{\Delta t})={\mathcal{O}}(h^{-4}{\Delta t}^{-1})$. The numerical wave equation solution is then evaluated with cost ${\mathcal{O}}(N)$ at any desired time and (not close) exterior target point using (4) with the quadrature of Sec. 3.1.

Remark 8**.**

*We chose a good diagonal shift empirically (by examining extremal eigenvalues) as $d=-0.25$: this vastly increases the corrector convergence rate in the case of ${\Delta t}\gg h$, yet does no harm in other situations. With this shift, $n_{c}=8$ was found to be sufficiently large to give stability and errors similar to that of a full implicit solution. Since $n_{c}\ll n$, the total corrector cost remains negligible compared to that of evaluating $\tilde{{\mathbf{g}}}^{k}$.

4 Numerical experiments

In this section we test the convergence of the full time-dependent BIE scheme for exterior wave equation BVPs for the torus and cruller surfaces. We set $n_{r}$ according to Remark 6, $n_{\varphi}=2n_{r}$, and $n_{c}$ and $d$ according to Remark 8.

4.1 Effect of coupling parameters $a$ and $b$

Choosing the torus surface, we pick an intermediate spatial order $p=6$, with a low-resolution discretization of ${n_{\phi}}=9$ by ${n_{\theta}}=6$ panels, thus $N=1944$, or $h=0.108$. For this experiment we fix a timestep ${\Delta t}=0.1$, with D-spline order $2q=p-2$ to match the spatial order. We solve the exterior BVP with data $g$ deriving from the unit-magnitude known solution given in the caption of Fig. 7, which, being a Gaussian pulse, dies away rapidly for $t>10$. The four panel pairs of Fig. 7 contrast the resulting behavior of the norm of the density, and of the pointwise error in $u(x,t)$, for the four combinations of constant coupling parameters $a\in\{0,1\}$ and $b\in\{0,2\}$. The latter value $b=2$ is chosen as the maximum principal curvature of the torus. These weights $a$ and $b$ are motivated in the introduction, summarizing [16]. Each combination results in a different behavior:

• (a)

$a=b=0$: Once the pulse has passed, the density $\mu$ grows asymptotically linearly in time. This secular growth is associated with a zero-frequency Neumann resonance [16]. Accurate evaluation of $u$ for long times at any target point is thus impossible, due to growing catastrophic cancellation (indeed, a growing error in $u$ is visible).

• (b)

$a=1$, $b=0$: The situation is better than (a), with the size of $\mu$ peaking then tending to a positive constant. The spatial function $\mu(x,t)$ (not shown) tends to a constant on $\Gamma$; this is consistent with the decay of $u$ in (4) since ${\mathcal{D}}$ acting on a constant vanishes in $\Omega$.

• (c)

$a=0$, $b=2$: The situation is similar to (b), except that the density continues to oscillate with smaller constant amplitude around 0.1. At late times a weak instability (exponential growth in error) is seen.

• (d)

$a=1$, $b=2$: In this case only, density decays exponentially after the peak of the pulse. The decay rate transitions to a slower exponential rate (for $t>19$) once $\mu$ has dropped to the typical size of the error.

Thus the story for (a), (b), and (d) is exactly the same as found previously for the sphere [16] (case (c) was not tested in that work). Only the last case leads to a density which dies exponentially, thus the possibility of potential evaluation with high relative accuracy. The last case also leads to the smallest $\mu$ values and lowest errors in $u$: the maximum $u$ error is $5\times 10^{-6}$ (compare (b) for which it is $4\times 10^{-5}$).

On a laptop with i7-7700HQ CPU, running a MATLAB implementation (using Parallel Toolbox with 8 threads, calling single-threaded Fortran90 for D-spline evaluation as in Sec. 2), the above calculation takes 16 seconds for the parallel assembly of the sparse quadrature matrices (requiring around 10 GB), then 7 s for each run (i.e. 0.025 s per time step), which is dominated by the single-threaded sparse matrix-vector products which evaluate $\tilde{{\mathbf{g}}}^{k}$ in the right-hand side of (38).

We remark that the exterior solution in this case decays much more rapidly than a typical solution of a scattering problem, which would be no faster than exponential, as determined by the scattering resonances. Therefore here one cannot demand a density decay rate matching that of the solution. In the more realistic scattering problem considered below, the decay rates are better matched.

4.2 Stability and convergence

Here we report on tests of accuracy and stability covering the $(h,{\Delta t})$ plane, for both torus and cruller domains. The same Gaussian pulse BVP as in Sec. 4.1 was used, and the same target point. The coupling parameters were fixed at $a=1$ and $b=2$ as motivated above. The temporal order $2q+2$ was set to match the spatial panel order $p$, with orders 4, 6, and 8 tested.

Fig. 8 shows the resulting numbers of accurate digits. Note that the peak in $u(x,\cdot)$ at the target $x$ is of order 1; thus panel (c) shows over 9-digit relative accuracy for $p=8$ at the minimum $h$. The main results are:

• •

In every case there is an unstable region (lower right, shaded in blue) in the approximate form of an inverse-CFL condition (9) for some ${\mathcal{O}}(1)$ constant $c_{\mbox{\tiny\rm iCFL}}$ that appears to be at most weakly dependent upon the order and upon the domain. Dashed red lines show a boundary of the form (9), to guide the eye; the contours and shading show some minor variations from this form.

• •

For the cruller at high orders, there is an additional stability constraint that ${\Delta t}$ not exceed a certain ($h$-independent) value ${\Delta t}_{\mbox{\tiny\rm max}}$. For $p=8$ it appears that ${\Delta t}_{\mbox{\tiny\rm max}}\approx 0.15$. For the torus, the scheme is stable for ${\Delta t}$ as large as 0.8 (a quarter of the domain diameter, and much bigger than the typical panel size).

• •

For any fixed ${\Delta t}<{\Delta t}_{\mbox{\tiny\rm max}}$, the scheme appears stable as $h\to 0$.

• •

When stable, convergence with the expected high orders is seen, with respect to both ${\Delta t}$ and $h$. The contours of constant error bend sharply from horizontal to vertical in each $(h,{\Delta t})$ plane, indicating the transition from ${\Delta t}$-dominated error (upper left region) to $h$-dominated (lower right). Fig. 9 quantifies the convergence rates: panels (a) and (b) show ${\Delta t}$-convergence at the fixed smallest $h$, matching the expected temporal order $2q+2$ in all stable regions. Because of the inverse-CFL condition, one cannot fix a small ${\Delta t}$ to test the $h$-convergence; instead we varied ${\Delta t}$ with $h$ along “$h$-dominated” curves, shown (grey dotted) in Fig. 8, which remains stable and in a region where contours are vertical. Panels (c) and (d) show that the resulting $h$-convergence is, barring some variation, consistent with ${\mathcal{O}}(h^{p})$, as expected for the interpolation order (Sec. 3.4). The ${\mathcal{O}}(h^{2p})$ spatial quadrature order is no longer seen.

Remark 9** (Other assessments of stability).**

*Stability was also assessed by extracting the growth/decay rate of $\|\mu(\cdot,t)\|_{\infty}$ near the end of the time interval, $T=50$. The results were consistent with the blue shading of Fig. 8 extracted simply from the error, so are not shown. More sophisticated attempts to assess stability did not prove useful. For instance, setting ${\mathbf{g}}^{k}=0$ in (38), the entire predictor-corrector scheme, followed by a backwards shift by 1 in $k$, can be viewed as a huge sparse square matrix acting on the density history vector $\{{\bm{\mu}}^{k},{\bm{\mu}}^{k-1},\dots,{\bm{\mu}}^{k-n}\}\in\mathbb{R}^{nN}$. We estimated extremal eigenvalues of this matrix via ARPACK [26] (eigs in MATLAB), but convergence was extremely slow, often taking much more time than the $T/{\Delta t}$ time-steps needed (we speculate that this is due to its block-companion matrix structure). *

Collecting the data in Fig. 8 required 966 solution runs, with between $63$ and $2000$ time-steps per run, and $N$ varying from $384$ to $13824$. This required 11 CPU-days on a server with two 14-core Intel Xeon 2.4GHz 2680v4 CPUs and 512 GB of RAM. Most of the time is taken by time-stepping (single-threaded sparse matrix-vector products). Because of our direct ${\mathcal{O}}(h^{-4})$ implementation, run time is dominated by the one or two smallest $h$ values. The RAM cost of ${\mathcal{O}}(h^{-4})$, dominated by forming and storing the set $\{A^{r}\}_{r=0}^{n}$, limited the minimum $h$ that could be tested; at our minimum $h$, 260 GB of RAM was used.

4.3 Plane wave pulse acoustic scattering example

Finally we demonstrate the solver for the Dirichlet (sound-soft) scattering application, using the cruller domain of Fig. 4(b). The boundary data is $g=-u_{\mbox{\tiny\rm inc}}$ on $\Gamma$, for the plane wave

[TABLE]

where ${\mathbf{d}}$ defines its direction and $T(t)$ its signal function. $u_{\mbox{\tiny\rm inc}}$ is a solution to the wave equation in $\mathbb{R}^{3}\times\mathbb{R}$. Fig. 10(a–b) shows two snapshots of the resulting physical (full) wave $u_{\mbox{\tiny\rm tot}}=u_{\mbox{\tiny\rm inc}}+u$, which vanishes on $\Gamma$ for all time, with ${\mathbf{d}}$ and $T(t)$ as stated in the figure caption. Note that the pulse $T(t)$ is about 7 times narrower than the pulse used in previous tests. Its full-width half-maximum is about 0.35, i.e. around 1/9 of the domain diameter.

Our experiment is done at spatial order $p=6$, with $2q=4$ to match in temporal order, and $n_{r}=2p+4$, $n_{\varphi}=2n_{r}$ as in Remark 6. Fig. 10(d) shows the pointwise convergence with respect to $h$, for ${\Delta t}=1.0h$. The median error is shown estimated by comparison against the smallest $h$ case (${n_{\phi}}=30$). It is quite consistent with 6th order convergence. The smallest-$h$ point shown, with median estimated error $2.5\times 10^{-6}$, used ${n_{\phi}}=27$ by ${n_{\theta}}=18$ panels, thus $N=17496$ spatial nodes. This took 44 minutes on the Xeon server mentioned above, using largely one thread. Of this, 23 minutes was sparse matrix filling (needing around two billion nonzero elements in $\{A^{r}\}$ matrices), and 21 minutes for the 279 time steps. 290 GB of RAM was used.

We note that if only 3-digit accuracy is needed, $n_{\phi}=9$ ($h\approx 0.11$) is sufficient, requiring only 1 minute, and 11 GB of RAM, on the laptop mentioned in Sec. 4.1.

An overall exponential decay of the solution is indicated by Fig. 10(c); notice that the signal also oscillates as it decays. Although domains without trapped rays have long been known to have exponential decay of the wave equation solution [29], the cruller has trapped rays—for example it has one in the plane $z=0$ reflecting off the five symmetric bumps encircling the hole. Assuming the hard-walled case is similar to the case of smooth potentials, recent analysis predicts exponential decay, in the high-frequency limit, at a rate determined by the smallest distance of scattering resonances from the real axis, which is controlled by the weakest trapped ray Lyapunov exponent [30] [14, Ch. 6]. The figure also plots the decay of the density 1-norm: it appears to be exponential at the same rate as the solution, which is thus the optimal rate. This indicates an absence of catastrophic cancellation in the solution representation, a major advantage of our formulation.

5 Conclusions and future work

We have implemented a recently-proposed combined-field time-domain integral equation formulation [16] to solve to high-order accuracy the Dirichlet BVP for the scalar wave equation in a general smooth exterior domain in $\mathbb{R}^{3}$. This showcases a method previously only studied for the sphere (exploiting separation of variables [16]). We provide evidence that in relatively general domains, as for the sphere, both of the coupling weights $a$ and $b$ must be positive to prevent long-lived resonances that cause catastrophic cancellation with conventional formulations. The retarded potentials are discretized in space using high-order quadratures for weakly-singular kernels, using new difference-spline temporal interpolants. The timestepping is explicit, via a predictor-corrector scheme. We verified the expected high order accuracy, explored stability in the $(h,{\Delta t})$ plane—which indicates an inverse CFL condition as found in the modal sphere case—and benchmarked a plane wave scattering example in which the exponential density decay rate appears to match that of the solution. In contrast to direct discretization methods, time-domain integral formulations need not have any CFL condition (upper bound on ${\Delta t}/h$), and we do not observe one.

The implementation presented is direct, since it needs ${\mathcal{O}}(N^{2})$ memory and effort per timestep, and assumes a structured grid of quad patches that is appropriate only for torus-like surfaces. In this setting, because of our high orders up to 8th, we achieve typically 5–9 digits of accuracy. Thus, for low-frequency-content data, needing $N\lesssim 10^{4}$, the scheme is quite useful for high-accuracy solutions on smooth surfaces.

The work suggests several directions for improvement.

• •

Handling more general shapes would require unstructured, adaptive triangles and/or quad patches. One route to handle the kernel singularity is then to borrow from on-surface quadrature schemes for harmonic layer potentials of this same singularity class [6, 42, 5].

• •

To get high accuracy evaluation near to the surface, a special quadrature scheme would be needed, building on a variety of existing schemes for the Laplace equation (see, for example, [23]).

• •

In order to effectively address large-scale problems, our marching scheme should be coupled with the fast algorithms mentioned in the introduction [3, 17, 9, 28, 35, 41].

• •

Extensions of our formulation and discretization to electromagnetic scattering problems using Debye sources, as discussed in [15, 20], seem possible.

• •

One may be able to optimize stability, or density decay rate, by a better choice of positive coupling parameters, or functions, $a$ and $b$.

• •

In our studies we tied the temporal order $2q+2$ to match the spatial order $p$, and observed a somewhat restrictive ${\Delta t}_{\mbox{\tiny\rm max}}$ for the cruller domain with $p=8$. Independently varying $2q$ and $p$ would expose whether this constraint is tied to the temporal or spatial high order.

• •

There are many opportunities for rigorous analysis of the discretization scheme.

Finally, we note that the inverse-CFL condition may pose a problem for spatially-adaptive quadratures in the low-frequency regime, since the large panels would place a lower bound on ${\Delta t}$ that would be inaccurate for the small features. It is possible that a modified temporal interpolant could remove this condition, or that locally-adaptive timestepping could circumvent it. In particular, since the temporal interpolants are constructed independently for each point on the spatial grid, it is possible to use different time grids wherever needed. However the stability of such a scheme requires further study.

Acknowledgments

We are grateful for useful discussions with Charlie Epstein, Maciej Zworski, and Andrew Hassell.