Highorder discretization of a stable timedomain integral equation for 3D acoustic scattering^{†}^{†}thanks: \fundingSupported in part by NSF Grant DMS1418871 and the U.S. Department of Energy ASCR Applied Mathematics program. Any conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the NSF or DOE.
Abstract
We develop a highorder, explicit method for acoustic scattering in three space dimensions based on a combinedfield timedomain integral equation. The spatial discretization, of Nyström type, uses Gaussian quadrature on panels combined with a special treatment of the weakly singular kernels arising in nearneighbor interactions. In time, a new class of convolution splines is used in a predictorcorrector 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.
Timedomain integral equation for 3D acoustic scatteringA. H. Barnett, L. Greengard, and T. Hagstrom
coustic scattering, timedomain integral equations, highorder methods
65R20, 65M38
1 Introduction
Problems involving the scattering of waves by obstacles have countless applications in science and technology. For timeharmonic data, numerical methods based on integral equations are welldeveloped and widely used. Advantages of using an integral equation formulation, rather than the partial differential equation itself, include superior conditioning when secondkind 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 timedomain 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 spacetime history, making them somewhat unwieldy in the absence of fast algorithms (see Remark 1 below).
Surprisingly, much of the literature is focused on equations involving the single layer potential (see, for example, the review by HaDuong [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 quasiexplicit marching scheme is developed for the time domain magnetic field integral equation, and [38], in which highorder accurate Calderonpreconditioned 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 highorder 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 nearsingularities associated with eigenvalues of the interior Helmholtz problem [24] [10, p. 48–49]. For timedomain 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 illconditioned. 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 inducing by an incoming acoustic wave impinging on an obstacle in a uniform medium. The function satisfies the equation
(1) 
for where is the domain exterior to a compact obstacle with boundary . In this work, for definiteness, we set unit sound speed and consider the Dirichlet (i.e. soundsoft) problem,
where on the surface, and is a given incident wave. Our methods could also be applied to the Neumann (soundhard) problem.
For a target point and target time , the retarded single and double layer potentials applied to a density , , are defined by [21, Sec. 10.7] [33, Ch. 1]
(2)  
(3) 
Introducing surface weight functions , , for (analogous to the combined field parameter in the timeharmonic case [24]), we represent the scattered wave in the form
(4) 
Taking the limit of (4) as , using the jump relation for the doublelayer operator [33, Sec. 1.3–1.4], we obtain a linear integral equation for , ,
(5) 
where is the principal value part of the double layer potential (3) spatially restricted to , and is the single layer potential (2) spatially restricted to . When is smooth, both and have weakly singular kernels with singularities bounded by , as would be the case for the Laplace equation. Note that (5) is of Volterra type in time, that is, at time the and operators involve only the density history for . The first term suggests that it is also of the second kind. However, even for the case , it does not fall into the standard Fredholm theory [21, Sec. 10.7], since 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 and and their effect on the long time behavior of is a focus of [16]. There it is shown that:
 i.

For convex obstacles, 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 is sufficient to remove the damaging zerofrequency interior Neumann resonance. For sufficiently large, the long time exponential decay rate of the density 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 where is its radius.
Below, we examine the effect of simple choices for and for scattering by a nonconvex obstacle. In particular, numerical experiments in Section 4 compare results for various and when is the boundary of a torus or a perturbed torus.
Since the boundary is timeinvariant, 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 highorder Nyströmlike approximation in space combined with a straightforward time marching procedure. For this, we introduce a set of surface nodes, , and a regular sequence of time steps , and represent the density by interpolation from its discrete values
To evolve forward in time, the th timestep consists of applying an explicit, fixed linear rule which gives the vector in terms of the history and the current data vector . Because of the strong Huygens principle, the number of previous times needed is essentially the diameter of the obstacle divided by . The contributions of these past values come from approximating the retarded spatial integrals in (5). We account for this history dependence by constructing matrices , and such that for each target node ,
(6)  
where is a measure of the spatial grid spacing, and the convergence order is determined by the scheme. Here is a set of nodes, described in Section 3, comprising all nodes far from , plus a set of auxiliary nodes designed to handle the weakly singular contribution near to (see Fig. 4(c)). The point samples in the expressions above involve the density and its time derivative at retarded times, e.g.
(7) 
so that the retarded evaluation times do not correspond to points at previously computed time steps, . 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
(8) 
where is the current time. However, we use different basis functions, . The temporal interpolating functions and some properties of the timestepping schemes defined in [11, 12] are discussed in Section 2.
The experiments in Sections 3 and 4 verify the highorder convergence of the fully discrete algorithm. In addition, the stability properties of the time marching scheme, in its explicit predictorcorrector form, are investigated, revealing a rather surprising “inverse CFL” constraint of the form
(9) 
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 work to factor, with subsequent solves required at each time step requiring flops. Clearly, an explicit method (including predictorcorrector iterations) is cheaper, as no matrix factorizations are necessary and the cost per time step is . Finally, in Section 5 we summarize our results, and point to future enhancements and generalizations of our method.
Remark (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 (Fast algorithms)
A critical bottleneck in using retarded layer potentials for the solution of the wave equation is that they require memory and work per timestep. Fortunately, over the last two decades fast algorithms have been developed which reduce both of these costs to . 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 timestepping timedomain integral equations are Galerkin methods, discussed extensively in HaDuong’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 timemarching 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 to be smooth, compactlysupported splines (both standard Bsplines and more exotic bases), thus restoring the finite time history of the layer potential operators. However, their approximations (8) are quasiinterpolatory; that is, the basis functions do not satisfy
(10) 
As a result, their proposed methods are limited to second order accuracy, even if Bsplines 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 Dsplines for short) [2] are defined as follows.
 i.

Let be an integer and be a set of nodes, with , and for each . These nodes need not be the from the previous section. Let be the set of nearest nodes to be used in the difference stencil for . There are two cases:
 a.

(Interior node, ): .
 b.

(Boundary node, ): .
 ii.

Let be the degree Lagrange interpolating polynomial defined by data, , on the stencil . That is, is the unique polynomial of degree satisfying for all with .
 iii.

On the interval the Dspline interpolant of the data , , is defined as the degree Hermite interpolant of and . Precisely, is the unique polynomial of degree satisfying for :
(11) As the interpolation operators are linear, and depends only on the stencil data, there exist degree polynomials , , such that
(12)  iv.

Let be the index of the rightmost node not larger than , i.e. such that . Then the Dspline basis function associated with node is defined as the piecewise degree polynomial
(13)
By construction , and, for uniformly spaced nodes in time, they are translates of a single simple basis function away from . In addition, the Dspline interpolant reproduces polynomials of degree , which by standard results implies accuracy of order for function values and for derivatives.^{1}^{1}1Here we will identify the interpolants by the polynomial degrees they exactly reproduce; that is the degree Dspline is the piecewise degree function described above.
In the uniform grid case, we plot the interior functions for to , along with boundary functions for , in Figure 1. A Fortran90 implementation (with MATLAB interface) of the above Dsplines 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 spatiotemporal problem, we illustrate the use of the above temporal basis functions to create a collocation timestepping scheme for a simple 2ndkind Volterra integral equation with “tophat” kernel,
(14) 
We assume is smooth, and is known for and smooth for all . As in Section 1 let , be the time grid. Let be the nodes and be the corresponding weights of a node Gauss–Legendre quadrature rule on . Enforcing (14) at , and inserting the quadrature,
(15) 
holds to high accuracy because the integrand is smooth. Let the unknowns be , for . Let be the degree Dspline basis, as defined in the previous section, for the regular grid of spacing , which one may think of as stepping backwards in time from the current time . The resulting interpolant is
(16) 
where, abusing notation slightly, we take for the prehistory to be populated with the known solution . Substituting this interpolant into (15) defines a lowertriangular Toeplitz linear system
(17) 
with weights computed by
(18) 
Note that in (17) the upper limit for can be replaced by , which is sufficiently large to capture all history dependence in (14).
The system (17) is naturally best solved sequentially for unknowns , i.e. by timestepping, since each row can be written
(19) 
where the “righthand side” is explicitly given in terms of data and previous values. Of course, in this simple example, the “solve” for at each time step is trivial:
(20) 
Figure 2 shows that this scheme achieves an empirical order . In this experiment, is constructed numerically (via an accurate quadrature) such that the solution is a known function. Note that, although in this test case is not strictly zero for , it is of order or less, so that initialization by zero data for the prehistory of induces negligible error.
The figure also shows stability for arbitrarily small , even with fixed number of quadrature nodes. The mean mesh spacing is the closest analog to in the full spatiotemporal 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 tested, the quadrature (18) is completely unresolved for each basis function , yet the scheme is accurate—since itself is resolved—and stable.
Remark
Recall that in the context where there are spatial variables, the analog of the solution (20) to (19) requires the solution of an linear system at each time step , i.e. an “implicit” time step, which can incur a large cost. As mentioned above, for the full problem we will focus on predictorcorrector schemes, which replace this by a fixed number of explicit steps and require only effort.
2.2 Modal problem on the unit sphere
As a second example of a simple Volterralike scalar problem, but one which is more directly related to the full problem we aim to solve, we consider our combined field equation with restricted to the amplitudes, , of a spherical harmonic expansion of the full solution on the unit sphere. As shown in [16], the integral equation to be solved is
(21) 
where denotes the degree 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 inverseCFL 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 point Gauss rule and the Dspline interpolant of as well as its derivative is evaluated at the quadrature nodes to produce a scheme
(22) 
with weights
(23) 
Choosing to be, for every that we test,
we solve to , for varying between and in increments of and for all from to . The method is deemed unstable for a given time step if the maximum value of the solution for any grows beyond the expected value by a factor of roughly . If we take the average mesh width to be and plot versus the minimum stable CFL number, , for we obtain the results shown in Figure 3. We see for small the minimum CFL number is approximately for and for . We note that for (a thorder method) the stability region shrinks considerably since a maximum limit of about also appears; that is, the method is stable in this test only for . Thus we restrict ourselves to methods with in the range 2 to 6 (orders to ), 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 doublelayer potentials (2) and (3). Finally we present the full interpolation scheme and Volterra timestepping for (5).
3.1 Retarded layer potentials for offsurface targets
The case of an exterior evaluation target , far from , is simple: for densities that are smooth with respect to both and time , the integrands are also smooth, and a standard quadrature scheme using the density interpolatory nodes will be accurate. In this work we restrict ourselves to smooth toruslike surfaces. Such surfaces can be parameterized by an infinitely differentiable, doubly periodic function , where . We now describe a simple composite (i.e. panelbased) rule to generate the nodes and weights , which we use in later tests, such that for any smooth ,
(24) 
holds to highorder accuracy. Here and are the partials of . We cover the parameter space with a uniform by grid of rectangular patches (“quads”). Each patch is covered by a tensor product grid comprising a node Gauss–Legendre rule in each of the two directions. Let , , where , be the images of these parameter nodes under the map . Figs. 4(a) and (b) show two surfaces with some of their resulting nodes. The corresponding weights are found as follows. Let , , be Gauss–Legendre weights on the interval . For surface node , let and be its parameter values, and be its indices in the two directions within the appropriate panel. Then
(25) 
The expected convergence order for (24) is , where is the resolution (combining [13, (2.7.12)] with a theorem on composite rules [13, Sec. 2.4]).
3.2 Retarded layer potentials for onsurface targets
When the target is on , as needed in the integral equation (5), then the integrand has the following type of weak singularity. If one smoothly parametrizes via local polar coordinates centered at the target point , for each (i.e. radial line) the integrand is times a smooth function of . 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 highorder 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 neartarget contribution. To handle the latter they exploit the fact that the polar metric cancels the singularity in the integrand, so integrate using auxiliary quadrature nodes on a polar grid centered at the target. For general highorder 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 .
Building on the above, we present a simple highorder 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 be a target node, in panel . This panel and its eight neighbors form a block of “near” panels, containing node indices (the black and green nodes in Fig. 4(a)), leaving “far” panels. With respect to each of the latter panels, the singularity of the kernel is distant, so their native quadrature nodes and weights (25) may be accurately used, as in the previous section. This explains the first term in our approximation
(26) 
where ; the second term is explained below. Note that (26) is of the form (6), with . The expressions for the other two kernels are analogous.
The second term in (26) accounts for the contribution of the near panel block via a new targetspecific set of auxiliary nodes. For this it is convenient to switch to the standard parametrization for this near block, with the target panel preimage being ; see Fig. 4(c). Precisely, there is a simple affine map to the global parameters and , where and are the toroidal and poloidal indices (i.e. integer coordinates) of panel within the panel grid. Then denote by the (dependent) map from to , which is the above affine map composed with the map . 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 ; see Fig. 4(c).
Specifically, let the polar coordinates be centered at the target preimage in , and consider one of the four triangles, , that lies in the angle range and whose far edge is given by in polar coordinates. For the singlelayer (2), the integral of a retarded density as in (7) over the image of this triangle on is
(27) 
where is the Jacobian of the map to the surface, and we abuse notation slightly so that means , etc. Let and , , be respectively the nodes and weights of a point Gauss–Legendre rule on . Let and , , be the nodes and weights of a point Gauss–Legendre rule on . Then (27) is approximated by
(28) 
where , etc, indicates the value at the auxiliary node indexed by and . Highorder convergence is expected for (28) since, although has a conical singularity at the polar origin, along constant rays the integrand times is smooth. Summing the four triangles, there are auxiliary nodes for each target point. The weights in (26) may be read off by associating each with a term in (28), and taking all factors except the density sample . The weights and in (6) are found in an analogous way.
Remark (order of convergence)
The convergence of this onsurface scheme is subtle, due to its split into far and near source panels. Consider the error due to the node smooth rule on one of the nearest “far” panels (i.e. just outside of the block): as , 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 , the formal order is low (). Yet, this matters little in practice because its prefactor is expected to be exponentially small in . This follows by analogy with 1D node Gauss quadrature on , for which the error is of order , where is the size parameter of a Bernstein ellipse in which the integrand is analytic [36, Thm. 19.3]. Since the nearest singularity is at scaled to the standard panel , solving gives , hence for , , and for , . Combining with the accuracy of the rest of the far panels, and assuming , one can summarize the expected error as . We postpone a rigorous analysis for future work.
Remark
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 selfinteraction 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 spatiotemporal 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 diffeomorphic to a torus. Given , a doubly periodic smooth height modulation function, then the map is given in Cartesian coordinates by
(29) 
Thus, the major radius is 1. The plain torus of Fig. 4(a) is given by the constant , and is close to being the most benign surface of this topology. The “cruller” of Fig. 4(b) has , 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 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 satisfy (1) in the closure of the exterior domain , for all , and let and indicate respectively the value and normal derivative on , then
(30) 
(Note that, since is a solution for all time , there is no need for initial conditions as in [33, Sec. 1.4].) We use this to test both the offsurface (Sec. 3.1) and onsurface (Sec. 3.2) quadrature schemes. We use for an exterior solution to the wave equation generated by a generic point source inside, but far from, the surface , namely
(31) 
with a source signal . Fig. 5(a) shows a snapshot of the resulting retarded doublelayer density for an onsurface target . In this section we show results only for the cruller, omitting the more accurate and predictable results for the plain torus.
Offsurface test. Figs. 5(b) and (c) shows (with dots) the convergence of (30) for an exterior target , which is a generic point a distance 0.36 from (see subfigure (a)). We use the smooth panel scheme of Sec. 3.1. For both shapes, gives panels with low aspect ratios. Recalling that each quadrature panel has nodes, (b) shows (for ranging from 384 to 11616), while (c) shows (for from 1586 to 24567). The horizontal axis shows the mean spatial node spacing, or resolution, which we define as
(32) 
Although there is variation, convergence is asymptotically consistent with the expected order . The variation is absent and the errors much smaller for the plain torus (not shown). For , 12digit accuracy is reached at the highest .
Onsurface test. A surface target 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 , and thus involved worstcase auxiliary triangle aspect ratios. For panel aspect ratios up to around 2, we found for the auxiliary quadrature that was adequate, so we fixed this ratio. We explore three choices of for each choice of . For , Fig. 5(b) shows that, for (), errors due to the singular scheme are negligible relative to the overall error for . For , a higher is needed: causes negligible errors for (). In general, for all but the largest available, the errors of the singular scheme are negligible for the choice .
Remark (Choice of singular scheme order )
Based on such tests, we fix for the torus and for the cruller. Thus the auxiliary scheme is of much higher order than the underlying panels. This allows us to explore larger without loss of accuracy or stability. At smaller one could reduce and hence the effort for the auxiliary scheme.
With converged as above, the convergence is as expected from Remark 3.2: errors drop with order roughly (dashed lines), until they saturate to a loworder convergence at around (for ) or (for ).
3.4 Interpolation and explicit timestepping
We now describe how the above spatial quadrature for retarded potentials is combined with the timestepping of Sec. 2 to solve the timedependent BIE (5). Enforcing the BIE on the time grid (as in Sec. 2.1), and on the spatial nodes , gives
(33) 
Now, fixing constant and , we approximate the action of the retarded integral operators on the density interpolated from the spacetime data by a set of matrices , thus
(34) 
where, analogously to Sec. 2.1, is large enough to capture all history dependence via Huygens’ principle plus the support of the Dsplines.
Spatiotemporal interpolation. Each above matrix is filled as follows. For simplicity, consider only the plain singlelayer () contribution in (33). Applying spatial quadrature (26), then time interpolation (8), gives
(35)  
The first term ( far from ) is already in the form (34), so the contribution to in this case is . We now use spatial interpolation on each time slice to turn the auxiliary term also into a weighted sum over . Let , for , be a set of basis functions that interpolate over the preimage of the near patch. I.e., for any smooth function ,
(36) 
holds to high order, where are the values at the interpolatory nodes (preimages of ). Let be the standard parameters of the th auxiliary node for the target node . Spatial interpolation of then approximates the 2nd term of (35) by
which is now also of the form (34). Proceeding as above with the other two terms of (33) (and recalling that has two terms (3)), finally gives the formula
(37) 
For the basis we use the product Lagrange basis for whichever of the nine panels lies in, and zero elsewhere. Precisely, let be the panel in which node lies, let and be its two index coordinates within that panel, and let be the parameter offset of panel relative to the target panel ( and are either , , or ). Then,
where are the usual 1D Lagrange polynomials for the Legendre nodes on . This form aids bookkeeping since it decouples all interactions into independent panelpanel pairs, each of which is either near or far. The different by blocks of given by all nodes in a single source panel interacting with all nodes in a single target panel may be filled together. In practice the nonzero entries of are precomputed once and for all, then for each target the sum over in (37) is performed as an efficient matrixmatrix multiplication (GEMM). Note that the nearpanel interpolation error is expected to be .
Remark
One might be tempted to interpolate with respect to the retarded densities such as ; however, this would fail to be highorder accurate due to the conical singularity around the target . Instead one must interpolate in both space and time, as above, since as a function of space and time is smooth.
Fig. 6 shows the sparsity patterns of a selection of the resulting matrices. As expected by Huygens’ principle, is concentrated around the diagonal, but as the time delay increases, the shell of influence spreads across the panels, departing at the most distant (furthest offdiagonal) panels.
Predictorcorrector scheme. Using the notation for the vector of densities at time step , we rewrite (34) as a linear system to be solved at each time step,
(38) 
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 th time step. Firstly a “predictor” is generated via a fixed order extrapolation rule in time applied to the density vector,
(39) 
where we choose to match the Dspline order. Here the are simply the values of the Lagrange polynomials associated with the time nodes evaluated at . Then is evaluated according to the righthand side of (38); this is the most expensive task. Finally “corrector” steps are performed on , each of which is a Jacobi iteration with shift ,as follows. The system matrix is split into the diagonal matrix and matrix , defined by
(40) 
The linear system (38) to be solved is then