FAST, HIGH-ORDER NUMERICAL EVALUATION OF VOLUME POTENTIALS VIA POLYNOMIAL DENSITY INTERPOLATION

. This article presents a high-order accurate numerical method for the evaluation of singular volume integral operators, with attention focused on operators associated with the Poisson and Helmholtz equations in two dimensions. Following the ideas of the density interpolation method for boundary integral operators, the proposed methodology leverages Green’s third identity and a local polynomial interpolant of the density function to recast the volume potential as a sum of single-and double-layer potentials and a volume integral with a regularized (bounded or smoother) integrand. The layer potentials can be accurately and efficiently evaluated everywhere in the plane by means of existing methods (e.g. the density interpolation method), while the regularized volume integral can be accurately evaluated by applying elementary quadrature rules. Compared to straightforwardly computing corrections for every singular and nearly-singular volume target, the method significantly reduces the amount of required specialized quadrature by pushing all singular and near-singular corrections to near-singular layer-potential evaluations at target points in a small neighborhood of the domain boundary. Error estimates for the regularization and quadrature approximations are provided. The method is compatible with well-established fast algorithms, being both efficient not only in the online phase but also to set-up. Numerical examples demonstrate the high-order accuracy and efficiency of the proposed methodology; applications to inhomogeneous scattering are presented.

1. Introduction.This paper considers the numerical evaluation of volume integral operators of the form over bounded domains Ω ⊂ R 2 with piecewise-smooth boundary Γ, where f ∈ C s (Ω), s ∈ N 0 := N ∪ {0}, is a given function (termed the source density), and where the operator kernel G k is the free-space Green's function (1) 0 (k |x − y|) , k ̸ = 0, for the Laplace (k = 0) or Helmholtz (k ̸ = 0) equation.In fact, the method we propose is considerably more general than the examples presented, extending to (i) volume operators corresponding to these partial differential operators that possess stronger singularities such as the gradient of V k as appears e.g. in [51], where ∇ y denotes a gradient with respect to the y variable and f ∈ [C s (Ω)] 2 , (ii) higher-dimensional analogues of these operators, and (iii) a variety of elliptic partial differential operators of interest in mathematical physics, e.g. in elasticity and fluids.Despite the broad utility of these operators, for example, to solve a nonlinear or an inhomogeneous linear PDE, until very recently the computation of volume potentials has been relatively neglected in the context of complex geometries; in what follows we outline some difficulties of this problem and present some alternatives to the volume potential.
There are several well-known challenges to be met for the efficient and accurate evaluation of operators such as (1.1) or (1.3) in complex geometries-some mirroring those of the more well-studied problem of layer potential evaluation and some unique to volume problems.As is well known, the free-space Green's function G k given in (1.2) features a logarithmic singularity (whose location is, of course, dependent on the target evaluation point x) that hinders the accuracy of standard quadrature rules for the numerical evaluation of the operator (1.1) or (1.3).Despite the requirement for some degree of local numerical treatment, the potential is a volumetric quantity depending globally on the source density, and as such it is especially important to couple its computation to global fast algorithms-of which a variety have been developed including H-matrix compression and their directional counterparts [14,15], the fast multipole method (FMM) [43], and, more recently, interpolated-factored Green function methods [12].Generally, it is also desirable to construct numerical discretizations of the operators which are efficient for repeated application (e.g. in schemes involving iterative linear solvers, Newton-type iterations for nonlinearities, or time-stepping) or with values that can be efficiently accessed (e.g. for the construction of direct solvers).
The prototypical, but by no means exclusive, use case for volume operators of the form (1.1) is to tackle the slightly more specialized problem of producing a particular solution to an inhomogeneous constantcoefficient elliptic PDE.Restricting attention to work in this vein where the resulting homogeneous problem is solved with integral equation techniques, there have been a variety of methods proposed for producing a valid particular solution, typically on uniform grids using finite difference, finite element, or Fourier methodssee [5] for more details.With the proviso that methods based on finite elements or finite differences are often of limited order of accuracy, they nonetheless do apply to piecewise-smooth domains.Fourier methods, in turn, relying on various concepts of Fourier extension to generate a highly-accurate Fourier series expansion of the density function (after which the production of a particular solution is straightforward) in particular have received significant attention as they allow the use of highly-efficient FFT algorithms.Such Fourier-based methods are either fundamentally limited to globally smooth domains [68,69] or have only been demonstrated on such domains [22,34].Recently a scheme has been proposed [33] that allows extension on piecewise-smooth domains and, by coupling to highly-efficient volume Laplace FMMs [30] achieves an efficient and high-order accurate algorithm for producing a particular solution to Poisson's equation.Another significant application of the volume operator (1.1) and (1.3) is in the numerical solution of Lippmann-Schwinger integral equations that arise in some formulations of inhomogeneous scattering problems, see e.g.[2,20,63,74]; for completeness we mention also the recent work [21] for the same problem that proceeds along a different path by performing a direct PDE discretization rather than using volume potentials.
In this work, the volume potential is computed at each target point by utilizing Green's third identity to relate it to a linear combination of regularized domain integrals (over a standard triangularization), boundary integrals, and a local solution to the PDE that can be explicitly computed in terms of monomial sources.The method relies on local interpolants of the density function and explicit solutions to the underlying inhomogeneous PDE associated with those interpolated density functions.
Although the proposed approach is in spirit a natural extension of the density interpolation method [31,38,[57][58][59], it shares several elements with previous methods.We first describe prior work in constructing local PDE solutions and then discuss the manner in which such solutions have historically been coupled to yield particular solutions valid on the entirety of the domain of interest.The first work we are aware of in which particular solutions were constructed using monomials is that of reference [6]; of course, the solutions are also polynomials.Subsequently, a variety of contributions studied combinations of polynomial sources and various elliptic partial differential operators; for monomials there exist solutions for Laplace [6], Helmholtz [37], and elasticity operators [52].A particular solution formula that is applicable to second-order constant coefficient partial differential operators featuring a non-vanishing zeroth-order term and a monomial right-hand side was presented in [26].Many of these methods can be considered to be closely related to the dual reciprocity method [54,56]-such schemes rely on the density function being globally well-approximated by a collection of simple functions (e.g.polynomials, sinusoidal, or radial functions), and implicitly rely on the accurate and stable determination of coefficients in such expansions.This latter task has not generally been successfully executed to even moderate accuracy levels for complex geometries and/or density functions.
Thus, for an arbitrary complex geometry, in order to use the simple basis function expansions of the density function envisioned in reference [6] for a stable, high-order solver, one is naturally led to the use of subdivisions of the domain, wherein global solutions are stitched together from solutions defined with local data.The case of a rectangular domain that can be covered entirely by boxes is considered in reference [41] to construct a high-order Poisson solver.That work, similarly to the present method, uses local polynomial solutions (corresponding to Chebyshev polynomial right-hand sides) and layer potential corrections to provide a global solution.As described there, the local solution defined with data supported on a given element undergoes jumps across the boundary of that element, and it follows from Green's identities that the particular solution, i.e. the volume potential, can be represented outside of the local element in terms of single-and double-layer potentials with polynomial densities defined explicitly in terms of the local element solution.This approach was also taken up using multi-wavelets on domains covered by a hierarchy of rectangular boxes in reference [8] as well as, in large part, in the three-dimensional works [46,53].Like these other works, the method proposed here uses a piecewise-polynomial approximation of the source.
We discuss next a few relevant distinctions between the present methodology and other more closelyrelated schemes.The scheme proposed in the present work differs from that introduced in reference [5] (see also [6]), where the volume potential is evaluated over irregular domain regions using a direct numerical quadrature of the singular and near-singular integrals.While the method presented in [5] requires only singular quadrature rules appropriate for the singular asymptotic behavior of the function G k , its applicability may not include some of the more severely singular kernels which the present methodology and the density interpolation method naturally treat [31,38]; additionally, the efficiency of that method suffers when gradients become large (e.g.|k| large) with a wider near-field volume correction area needed.We turn next to methods that like ours leverage a polynomial solution to the underlying PDE and transform volume potentials to layer potentials.Unlike the previously-mentioned works [8,30,41,46,53] that rely on the method of images and/or pre-computed solutions for singular and near-singular evaluation points-since the grid and therefore the target locations have known structure in those settings-the present work applies to unstructured meshes wherein evaluation points may lay in arbitrary locations relative to element boundaries.Our method shares some elements with the work [66] that has contemporaneously appeared: both use Green's third identity for volume potential evaluation over unstructured meshes, but, unlike that work, the class of near-singular layer potential evaluation schemes which we employ extends to three dimensions, treats operators with stronger singularities, and is, in principle, kernel-independent (see e.g.[31]).We stress, however, that with regard to choice of methods for layer potential evaluation the scheme is agnostic.Indeed, alternative techniques such as [1, 9-11, 13, 44, 47, 67, 80] can be used instead of the general-purpose density interpolation method [31].
An important novel feature of this method is that it results in a significant reduction in the number of singular corrections compared to previous volume potential schemes (even as the total number of volume corrections remains the same).This is made possible by the use of Green's identity over a large region Γ = ∂Ω (illustrated in Figure 1) as opposed to individual element-wise boundaries, with the effect that, on the one hand, near-singular layer potential quadrature (e.g. with the density interpolation method) is only required in a thin region abutting Γ while, on the other hand, volume quadrature for resulting regularized volume integrals can be effected to high-order accuracy with standard quadratures.
To support this last claim the paper presents a detailed error analysis culminating in the main Theorem 6.6 with high-order (h-refinement) estimates that are observed in the experiments performed for moderate interpolation degrees n ∈ {0, 1, 2, 3, 4} (certain numerical instabilities may be more significant for large interpolation degrees, an issue that is discussed in detail in Section 8 with a prescription outlined for managing the instability).For context on the significance of these results, the guarantees established here-most precisely, the main theorem's supporting Lemma 6.9 that establishes high order accuracy (an order of accuracy that is, in a sense made precise, optimal) as the mesh size h → 0 when a quadrature rule of fixed order m is used for all elements not containing the evaluation point-must be understood in contrast to the state of the theory for methods for nearly-singular integral evaluation in widespread use.Indeed, in the community it is routine, both for volume and layer potentials, to refer to the 'near-field' as a neighborhood of diameter proportional to the element / local patch / panel size h; the complement of this neighborhood is called the 'far-field'.Common practice is to use specialized corrective quadratures only over an O(h)-sized 'near field' in conjunction with a fixed-order (not spectral) quadrature rule for the 'far-field'.Unfortunately, such methods are often accompanied by unjustified claims of high-order accuracy because the 'far-field' encroaches closer to the singularity under h-refinement and errors in the 'far-field' integrals are not controlled to high-order.Simply put, as the mesh-size h → 0, either (i) the order of accuracy of far-field quadratures is expected to eventually degrade, or (ii) if, instead, true high-order accuracy is desired then Fig. 1: Comparison of the regions for which specialized quadrature is needed in two volume quadrature methods.Top: With the use of standard quadrature rules for V k [f ], nearby triangles for which the volume integral will be inaccurate for specific evaluation points collocated at a quadrature node and marked with a red star.Computation of the potential at all nodes displayed in the volume requires developing many (near-)singular quadrature rules for every triangle.Bottom: In the proposed methodology, all volume evaluation points which require near-singular layer potential evaluation and indeed specialized quadrature of any kind are marked in blue; no specialized quadrature is needed for evaluation points inside the domain.
. corrections are made on a growing number of elements or points (with the need for the inclusion of these corrections impacting on efficiency or even leading to superlinear complexity [19,79]).For example, the representative reference [67] even notes "crucially, this scaling of parameters with h also leads to the desired O(1) work per target in the local correction step".A certainly incomplete list of further recent examples in this 'locally-corrected' vein include [16-18, 42, 72, 80] and the volume potential schemes [5,66].Nevertheless, the practical effect of this issue may be concealed by the higher-order quadratures frequently used in the far-field so that significant impacts on accuracy are not necessarily seen in practice (reference [66, Fig. 9] does report seriously sub-optimal convergence rates in light of the theory in this paper).In connection with this discussion it is interesting to observe that in the boundary element literature [65, §5.3] similar problems are controlled through use of variable quadrature orders, that is, the use of quadrature orders that increase logarithmically with h-refinement.The remainder of this paper is structured as follows.An introduction to basic identities and background information on quadrature is presented in Section 2. Section 3 presents a general overview of the method including its compatibility with fast algorithms.The resulting complexity analysis is given in Section 4. Section 5 focuses on the construction of the interpolation polynomial f n (the conditioning of the linear system used to recover the interpolant is analyzed in Appendix A).A detailed error analysis presented in Section 6, whose outcome is Theorem 6.6, establishes the regularizing effect of the polynomials on the domain integrals and ultimately gives rise to a provably high-order accurate numerical scheme.Finally, numerical examples of the proposed approach are presented in Section 7, some limitations due to certain instabilities and methods to manage them are given in Section 8, and a few comments on future work are given in the conclusion.

Preliminaries.
In what follows we briefly outline the principles underlying the source approximation and singular integration strategy employed in this paper.The method relies on polynomial approximations of the density function but is somewhat agnostic to the interpolation strategy.We present a specific strategy wherein approximation is performed by Lagrange interpolating polynomials over the entirety of an element K in a triangulation T h .The method locally constructs, on every element K, a piece-wise (regularizing) polynomial approximation to the source, but, as will be explained in detail, is represented in a global basis on Ω.Without loss of generality, we assume the domain's centroid lies at the origin, which mitigates the scale of this (polynomial) basis over Ω.
Next, the method builds associated regularization polynomial solutions that are particular solutions to the PDE with the interpolant as a right-hand side, and the interpolants and solutions are then used to regularize the volume integral in (1.1).Precisely, we first explicitly construct a family of polynomial PDE solutions Φ n (•; K) : R 2 → C, the family being parametrized by K ∈ T h , that satisfies where f n (•; K) is the polynomial of total degree n that interpolates f on a discrete set of points within the (possibly curvilinear) triangle K which is closest to the target point x ∈ R 2 (typically, x ∈ K); details about the construction of f n and Φ n are given in Section 5 and [3] respectively.Employing Green's third identity and the solution Φ n we arrive at the relation where with γ(x) denoting a function proportional to the interior angle at x ∈ Γ that is equal to the familiar In addition, the regularization of the operator W k can be achieved similarly by interpolating the vector density f component-wise.Indeed, letting f n (•; K) denote the corresponding polynomial interpolant of f at/around the target point x ∈ R 2 , and letting Φ n (•; K) denote a solution of the vector PDE we obtain from the identities and Just like the case of V k [f ], the interpolation property of f n enables the evaluation of the regularized volume integral in (2.5) using a quadrature rule independent of the target point location.Using (2.5) in conjunction with existing methods for evaluating the layer potentials we obtain an accurate and efficient algorithm for evaluating W k [f ].We mention in passing that the divergence theorem can be useful in reducing computational effort when both V k and W k are needed.
Remark 2.1.In general, in the derivations and in the analysis that follow, we will assume smoothness of the density functions f and f over Ω, e.g. for a given integer n ∈ N 0 we assume f ∈ C s (Ω) and f ∈ [C s (Ω)] 2 , for a certain integer s > n (specific s values can be found in Theorem 6.6 but we do not claim the required value of s = n + 3 is in general tight; if all elements in the mesh are convex the required regularity for the theory is only s = n + 1 and should be sharp).In addition, the methodology applies in the event the source densities are only piecewise smooth over Ω, in which case Ω must be considered as a union of disjoint domains where f (resp.f ) is smooth.Thus, provided the location of possible discontinuities of a density function is known, such concerns can be addressed using the proposed technique.
2.1.Elementary numerical integration.We will work over a tessellation of the domain that consists of both straight and curved triangles.Let T h = {K ℓ } L ℓ=1 be a triangulation of Ω, i.e., Ω = ∪ L ℓ=1 K ℓ and Kℓ ∩ Kℓ ′ = ∅, ℓ ̸ = ℓ ′ , which we assume throughout consists of elements with maximum diameter (in the sense of the metric distance) h > 0 and with a minimum diameter of inscribed circles of ρ > 0; at times, local element diameters h K and inscribed circle sizes ρ K will be used.
2.1.1.Mappings for straight and curvilinear triangles.All elements are mapped from the closed reference triangle K with vertices and the transformation to physical coordinates is denoted by T : K → K. Here, and throughout this paper, superscript-ed hat notation, e.g.T , denotes quantities in reference space on the reference triangle K; undecorated symbols, e.g.T or K, denote maps into or quantities in physical space.When K is a straight triangle, T is an affine transformation that maps the vertices of K to the vertices of K.When K is a curved triangle with vertices v 1 = (v 11 , v 12 ), v 2 = (v 21 , v 22 ), and v 3 = (v 31 , v 32 ) that has exactly one curved edge (i.e. on the boundary), say along the edge connecting v 1 to v 2 , we employ, as recently suggested [5], the blending mappings [39,40] defined on the standard simplex ∆.To do this, let T ∆ denote the affine map T ∆ : K → ∆ and let the curved element map T : K → K be given by T = T γ • T ∆ , with .
Here, γ = (γ 1 , γ 2 ) : [0, 1] → Γ ∩ K is the parametrization of the coordinates of the curved edge connecting v 1 and v 2 that satisfy γ(0) = v 1 and γ(1) = v 2 .It can easily be seen that T is a C 1 -invertible map of K onto K for each element K; we denote {T ℓ } L ℓ=1 the set of all such curvilinear and affine mappings for elements in T h .2.1.2.Quadrature for boundary integrals and for smooth functions on mapped triangles.The regularized volume integral operators in (2.3) and (2.5) require a quadrature rule for smooth functions defined on the elements K ℓ ∈ T h , ℓ = 1, . . ., L. For sufficiently regular functions ϕ : Ω → C, a volume integral can be expressed via the sum, over all L elements of the triangulation T h , (2.6) where J ℓ is the Jacobian matrix of the transformation T ℓ prescribed above.Quadrature for each of the integrals over K in (2.6) is performed via the q n -numbered (see (3.3)) Vioreanu-Rokhlin quadrature nodes [75] I n = x j : x j ∈ K qn j=1 with associated weights {ω j } qn j=1 .A Vioreanu-Rokhlin quadrature rule is capable of exactly integrating polynomials of total degree m (values of m are tabulated [75]).Use of other quadratures is possible but we prefer Vioreanu-Rokhlin nodes because they doubly function for us as relatively well-conditioned interpolation nodes (see Appendix A and Figure 8).Letting N = Lq n , real-space nodes are obtained via the mappings {T ℓ } L ℓ=1 yielding the node-weight set of global quadrature nodes and weights, and the quadrature approximation can thus be written as (Note that the notation introduced in (2.8) allows us to refer to quadrature approximations Q E applied to integrals over subsets E ⊂ Ω.) The layer potentials in (2.3), on the other hand, can be accurately evaluated by means of existing techniques [31,38,57,58] (or the other layer potential schemes given in the introduction).For sufficiently regular trace functions φ, ψ : Γ → C, such techniques yield high-order approximations (2.9) where B Γ stands for (any) sufficiently accurate approximation of the boundary integrals.When x is wellseparated from Γ, of course, standard quadrature rules are sufficient.
3. Method description and algorithm sketch.In this section, we outline a numerical procedure that exploits (2.3) to develop an efficient, highly accurate method to numerically evaluate (1.1) everywhere in Ω, i.e., here we assume that x ∈ K ∈ T h ; the method directly extends (component-wise) to allow use of (2.5) in evaluating (1.3).The generalization of the procedure to evaluate volume potentials over narrow regions outside Ω (where the volume integrals remain nearly singular) is discussed in Remark 3.1.
To explain more concretely, we introduce the notation used in the remainder of the paper: we will make use of the standard multi-index notation where, for any α = (α 1 , α 2 ) ∈ N 2 0 , we set α! = α 1 !α 2 !, |α| = α 1 + α 2 , y α = y α1 1 y α2 2 when y = (y 1 , y 2 ) ∈ R 2 , and 0 .Our goal is to use a well-conditioned nodal basis set to produce a Lagrange interpolation polynomial f n (•; K) : Ω → C that regularizes the singular and near-singular integrals over and near to the element K ∋ x, K ∈ T h .(The same local polynomial density interpolant f n (•; K) will be used for all points x ∈ K.) A significant difficulty in using local interpolants with (2.3) is that the volume integrals (2.3) depend on the target point x ∈ K explicitly through the Green's function as well as implicitly through f n (•; K).Clearly, evaluating the volume integral over Ω at all target points x ∈ K for every K ∈ T h by means of the quadrature rule (2.8), would result in a highly inefficient algorithm even using fast methods.These problems are also present for the layer potentials in (2.3).To mitigate this cost, we express, using the binomial theorem, every local polynomial f n (•; K) in a fixed, global basis of normalized monomials, i.e., (3.1) f n (y; By doing so, V k [f − f n (•; K)](x) can be evaluated by means of the volume quadrature but via the intermediate application of the volume potential to each monomial p α .The application of Q Ω to G k (x, •)p α can be precomputed for a given fixed Ω and then re-used for each new source f , merely with different weights {c α [f ](K)} |α|≤n -thus leading to an algorithm with offline (f -independent) and online (f -dependent) components.The offline computations can be accelerated with fast algorithms, while the f -dependent weights themselves are efficiently computable in the online phase; details on computational costs are provided in Section 4.
The coefficients of the interpolant are obtained in a shifted and scaled coordinate system, because, as we explain in detail in Section 5, in Lagrange interpolation it is useful (for stability reasons) to solve the multivariate Vandermonde system in coordinates where the triangle fits inside a unit ball (see Figure 8 in Appendix A).The interpolation coefficients {c α [f ](K)} |α|≤n in (3.1) are then contained in the vector where the vector d K ∈ C qn is filled with samples of f at the q n quadrature nodes on K, and V K ∈ R qn×qn is the Vardermonde matrix in the shifted and scaled coordinate system.Further details on both the linear system V K c K = d K associated with (3.2) and the change of coordinates map represented by the matrix P K ∈ R qn×qn are provided in Section 5 (more precisely, P K is given by (5.7)).Unfortunately, the improvements to conditioning of Vandermonde systems that are obtained by the employed translation-and-scaling method do come with them a potential instability in the change of coordinates map P K .A viable strategy to manage this instability is discussed in Section 8 by using domain decomposition of the domain such that appropriate norms of P K are controlled.In most practical cases and for interpolation degrees n ∈ {0, 1, 2, 3, 4}, however, a simple re-centering of the domain suffices to manage this instability.
The square shape of the Vandermonde matrix V K in (3.2) displays an advantage of utilizing the Vioreanu-Rokhlin quadrature rule, which is the fact that the number of quadrature points q n coincides with the cardinality of the monomial basis {p α } |α|≤n , i.e., (3.3) q n := Once the coefficients {c α [f ](K)} |α|≤n in (3.1) are found, the polynomial in (2.1) is, by linearity, where the polynomials {P α } |α|≤n satisfy the inhomogeneous Laplace/Helmholtz equation: The construction of the polynomial PDE solutions P α corresponding to the monomial sources p α is based on the procedures presented in [3] and those solutions are used here.
With the polynomial f n (•; K) and polynomial PDE solution Φ n (•; K) in hand, we turn to their use in (2.3) and (2.5).The numerical approximation of the regularized volume integrals in (2.3) and (2.5) is performed in this work using general-purpose quadratures for smooth functions, both in the online and offline phases of the algorithm.We perform two approximations: first, the contribution from the element K containing the evaluation point is neglected and then smooth quadrature (i.e. a quadrature meant for smooth functions that is oblivious to any possible singularities present) on the remaining region is performed.In the description and analysis of each of these approximations that follows, it will therefore be useful to introduce the operators In the sequel we denote by the approximations of these operators on a triangulation T h using reference element quadratures that can exactly integrate polynomials of maximum total degree m ∈ N 0 , see e.g.[75,78].In fact, these quadratures are independent of x, implying they are compatible with fast algorithms, and the quadrature over the whole triangular mesh for V k can be written in the Section 2 notation for generic volume quadratures: that, by leveraging the expansion (3.1) of f n in the global basis {p α } |α|≤n , leads to the following numerical method for computing Here, the first approximation commits regularization error, the second approximation commits quadrature error, and the third equality identifies how the volume potential approximation which is to be analyzed can be efficiently computed, for a fixed geometry, in an online (f -dependent) / offline (f -independent) setting.
It may be useful at this stage to point out that the approximations in (3.8) can be expected to be of highquality in view of the fact that f − f n (•; K) is small (in a sense made precise in Lemma 6.8 and Theorem 6.6) over the region near to x where G k is singular or experiences large gradients.Based on the intuition that since f − f n (•; K) and its derivatives are small near x the integrand should be significantly smoother than the kernel G k (•, x) itself, we simply integrate with standard quadrature rules Q Ω (2.8) intended for smooth functions ϕ defined over Ω. Precise analysis of this intuition is provided in Section 6, with convergence rates for the approximations in (3.8) established.(Analogous statements to (3.8) for W k can be written, and the error analysis for these approximations is also given in Section 6.).
Combining (3.8) with (2.9) in (2.3) it then follows that V k [f ](x) for x ∈ Ω, can be approximated as As usual, the symbol ∂ ν in (3.9) denotes the derivative along the exterior unit normal to the boundary Γ.
Clearly, being that the set ) is independent of f , it can be precomputed and used for any suitable input density.
Remark 3.1.The regularization procedure outlined above for x ∈ Ω can be directly extended to arbitrary target point locations.Indeed, for target points x ∈ R 2 \Ω lying close to Γ, the same numerical regularization procedure can be followed, using the Lagrange interpolation polynomial f n (•; K) of the density f on the quadrature nodes contained in a triangle K ∋ x ⋆ where x * := argmin 4. Complexity analysis.The efficiency of the proposed methodology becomes evident from considering the evaluation of the approximation (3.9) to V k [f ] at all the N -numbered volume quadrature nodes {ξ j } N j=1 in the mesh.Table 1 provides an overview of the computational costs associated with these evaluations, assuming that the FMM is utilized to accelerate the volume and boundary integral computations.(Clearly, the layer and volume potential pre-computations inherent in (3.9), i.e. the computations listed in the f -independent rows in Table 1, can be computed in an embarrassingly-parallel fashion across the polynomial degree multi-index α.) The costs in Table 1 for evaluating ) are the costs of FMM-accelerated summation at the N volume target points; the FMM-accelerated general-purpose DIM [31] is utilized to evaluate the boundary integrals B Γ , with the boundary discretized using N b ∝ N 1/2 points.Turning to the evaluation of the coefficients 2) entails the LU-factorization of the Vandermonde matrix V K for each of the L elements K ∈ T h , which amounts to a cost of O(q Table 1: Computational costs for both the online (f -dependent) and offline (f -independent) components of the method; both the online and offline components are required.Here N = Lqn is the number of volume evaluation points in a mesh with L elements.The cardinality of the polynomial basis is denoted by qn per (3.3).
map P K for all of the L elements costs a total of O(q n N ) in view of the prescription (5.7) for this map.
Precomputing C K for every mesh element thus costs O(q 2 n N ), while evaluation of Overall, the operation count estimate shows that, for a given interpolation degree n ∈ N 0 , the VDIM methodology achieves quasi-optimal complexity which is confirmed by the numerical results in Section 7.
5. Interpolant construction.We construct a Lagrange interpolation polynomial f n (•; K) of total degree at most n ∈ N 0 that interpolates data at a prescribed nodal set (with q n defined in (3.3)): For detailed discussions concerning multivariate polynomial interpolation see [36,55,64]; we discuss the interpolation problem in more depth later but remark here only that given a nodal set I n the solvability of the Lagrange interpolation problem (in which case I n is called poised) in multiple dimensions is not assured.It suffices that the associated multivariate Vandermonde matrix is nonsingular, a requirement which is achieved, as we will prove, by the choice of the interpolation set I n as the mapped Vioreanu-Rokhlin nodes [75], i.e., I n = T ( I n ), T : K → K, with the reference space nodes I n = { x j } qn j=1 ⊂ K being precisely the Vioreanu-Rokhlin nodes.(The multivariate Vandermonde matrices associated with these nodes on the reference triangle have favorable condition numbers [75, Tab.5.1-2], at least on a certain basis, suggesting their use here.) As a matter of practicality, one can observe on the one hand that directly interpolating using the physical nodes leads to highly ill-conditioned matrices (in view of the vast differences in scale in the columns), while on the other hand interpolation on the reference simplex K, where conditioning can be considered ideal, is simultaneously not trivially compatible with the use of normalized monomials {p α } |α|≤n introduced in (3.1) throughout the domain (especially for curved triangles on the boundary) and can also involve a potentially inefficient change of basis calculations.We strike a middle ground by interpolating on a translated and scaled triangle , where c K and r K are K-dependent translation and scale factors, respectively, whose prescription is given in Figure 2. We utilize translation formulae for the purposes of using a global basis.An important consequence for the theory presented in Appendix A is that this translation and scaling operation results in the triangle K that possesses a minimum bounding circle of unit radius.A different translation and scaling procedure, and not with the use of translation formulae to retain a global basis, has been employed in the recent work [66].
The interpolating polynomial is thus obtained by solving the multivariate Vandermonde system (5.1) . ., f (x qn )] ⊤ ∈ C qn the vector of data samples at I n .Here V is the multivariate Vandermonde matrix corresponding to a normalized monomial basis set {p α } |α|≤n and the set { x j := ( x j,1 , x j,2 )} qn j=1 of translated-and-scaled points: (5.2) . . . .
Appendix A develops bounds on the condition number of V K (denoted there simply as V), with the dependence tightly linked to the mesh element quality.
The resulting translated-and-scaled interpolation polynomial by construction satisfies f n ( x j ; K) = f (x j ), j = 1, . . ., q n , and can be expressed as where the coefficients { c β [f ](K)} |β|≤n , sorted in (say) lexicographical order, are contained in the vector c K ∈ C qn solution of (5.1).Therefore, by setting x = r −1 K (x − c K ) above, we get that the sought Lagrange interpolation polynomial is given by (5.4) The needed separable expansions for f n (•; K) in (5.4) as well as for its associated polynomial particular solution Φ n (•; K), can be found by using the binomial formula that results in with the coefficients given by (5.7) The potential W k in (1.3), which features a vector-valued function f in the integrand, requires a vectorvalued interpolant f n (•; K) that can be constructed component-wise; explicit expressions for f n (•; K) and Φ n (•; K) are omitted for brevity.
6. Numerical quadrature accuracy for regularized volume integrals.This section examines the regularizing effect of the density interpolant on the volume potential evaluated at some point x ∈ K, K a (curvilinear) triangular element.Its relevance is that part of the method's efficiency owes to the avoidance of the use of (near-)singular volumetric quadrature in the entirety of the domain, as explained above.Indeed, we develop the arguments that allow us to restrict the use of nearly-singular quadrature (for layer potentials) to a bounded region close to the boundary and still control the error introduced by using generic quadratures meant for smooth functions on regularized volume integrals.The outcome of this analysis is error estimates, given thereafter in Theorem 6.6.These error estimates are optimal in the sense that the order of convergence established is the same as that for the regularization error, for example that is, the order of accuracy obtained by the use of exact evaluation of the nearly-singular volume integrals.Moreover, they are closely matched by the results from numerical experiments given in Section 7.
6.1.Main convergence result.We start by requiring our family of meshes T to satisfy the usual quality constraints [65].Definition 6.1.A family of triangulations of Ω denoted by T = {T h } h>0 with h > 0 denoting the maximum element size in T h is called shape-regular if the shape-regularity constant is finite and quasi-uniform if the quasi-uniformity constant For such mesh families, Theorem 6.6 applies for evaluation on subsets of Ω that are interior to an element with respect to other elements of the mesh, an idea made precise in the next definition.Definition 6.2.Consider a family of triangulations T = {T h } h>0 of Ω, and let S h = K∈T h ∂K.A family of evaluation point sets E = {E h } h>0 , E h ⊂ Ω for all h > 0, possibly intersecting the boundary Γ = ∂Ω, is called well-separated with respect to T if A typical example is evaluation at the point set consisting of (interior) Vioreanu-Rokhlin interpolation/quadrature nodes: the theorem thus estimates the error in the volume potential evaluated at all the interpolation/quadrature points ξ j N j=1 , N = Lq n .Indeed, any set of interior points on K leads to a wellseparated family of evaluation points over the triangulation, as the following proposition shows.Another use-case of the definition would be the evaluation at point-sets laying on the boundary Γ = ∂Ω, e.g. for the solution of boundary value problems.
Now, using the fact that ρ −1 K ≤ κ T /h K and that h −1 K ≤ q T /h, we obtain ρ −1 K ≤ κ T q T /h, which by the shaperegularity and quasi-uniform assumptions on T is a finite constant.It hence follows from the inequalities above that Therefore, taking infimum over all the points x ∈ E h and y ∈ S h , and then infimum over all the triangulations T h , we arrive at where d T ,E is defined in (6.3).The proof is now complete.
Remark 6.4.Similar statements can be easily made for triangulations including curved triangles as well as for families of evaluation point sets E that lay on Γ (and which are yet separated from S h \ Γ).
A final ingredient for the statement of the theorem is an assumption on the existence and boundedness of Lagrange polynomials on mesh elements in the family.Assumption 6.5.For a given family T = {T h } h>0 of triangulations of Ω, the Lagrange interpolation polynomials λ i (y) = λ i (y; K ℓ ) (1 ≤ i ≤ q n ) exist on each element K ℓ ∈ T h ∈ T and, further, are uniformly bounded; that is, (6.4) sup with H(K ℓ ) being the convex hull of the vertices and the set of interpolation nodes I n .
This is essentially the same concept as the Λ-poisedness introduced in Appendix A for Theorem A.2, for the purposes of analyzing the conditioning of the polynomial interpolation problem.Indeed, Assumption 6.5 is satisfied for shape-regular families T of straight triangles, a fact that can be seen using elements of the proof of Theorem A.2: the estimate (A.6) together with shape-regularity providing an upper bound on the quantity R = h K /ρ K there.We do not pursue here guarantees of the assumption over curved triangles.
We are now in a position to state the theorem.
Theorem 6.6.Let the respective (positive weight) quadrature and interpolation degrees m, n ∈ N 0 , m > n, and the wavenumber k ∈ C be given and let T = {T h } h>0 denote a family of triangulations of the connected domain Ω that is shape-regular and quasi-uniform and satisfies Assumption 6.5.Let E = {E h } h>0 be a family of well-separated evaluation point sets for T , and denote by ), s = max{n + 3, m + 1}, with interpolation in some K ∈ T h enforced on the associated interpolation node-set where C V and C (2) V are positive constants independent of h but dependent on Ω, κ T , q T , f , k, m, n, d T ,E and Λ(n, T ); the constants depend linearly on Λ(n, T ).Likewise, for the W k operator, for all x ∈ K ∩ E h it holds that where C W and C W are positive constants independent of h but dependent on Ω, κ T , q T , f , k, m, n, d T ,E and Λ(n, T ); the constants depend linearly on Λ(n, T ).
Fig. 3: Depiction of the meshed domains employed in the proof of Theorem 6.6.
Remark 6.7.The Vioreanu-Rokhlin quadrature rules serve also as interpolation schemes, allowing more concrete error expressions where m and n are related.For interpolation degrees n ≥ 2, the tabulated Vioreanu-Rokhlin rules satisfy m ≥ n + 2 [75, Tab.5.1], so the error estimate (6.5) in such cases simplifies to for some constant C V > 0. Similarly, for n ≥ 2 the the error estimate (6.6) simplifies to for some constant C W > 0. Estimates (6.5) and (6.6) yield concrete error estimates when n < 2 via the relations [75, Tab.5.1] m = 1 when n = 0 and m = 2 when n = 1; note that the case n = 1 happens to be, comparing to Figures 4 and 5, the case where super-convergence is observed in numerical experiments for the V k operator.As a consequence, the analysis that follows shows, for V k and n ≥ 2 or for W k and n ≥ 1, that these estimates for the numerical Vioreanu-Rokhlin quadrature evaluation of the regularized volume integral yield precisely (modulo a | log h| loss for W k with n = 1, m = n + 1) the order of accuracy in (6.14) and (6.15), respectively, of the regularization error alone-that is, the approximation error introduced by neglecting the regularized weakly-singular integral over K ∋ x.
In preparation for proving Theorem 6.6, we note that we must account for contributions from not only singular and near-singular integration but also more distant regions, and for these purposes it is helpful to handle integration regions separately.We thus define (6.9) 3).Here and in the sequel, the analysis that follows (1) establishes that the contribution of the first integral for the potential is small (the regularization error), (2) shows that employing standard quadratures for the second integral leads to errors that depend on the interpolation degree and that match the order of accuracy of the regularization error, and (3) develops estimates for the accuracy of standard quadrature rules when applied to the last integral.Based on the estimates that we develop here we could, for evaluation at a given x, neglect all elements in the region N h (x) and retain an optimal order of convergence; in practice, for both simplicity and accuracy, only the contribution from K will be discarded, consistent with our definitions (3.7) of the punctured Green's function and (3.6) of the punctured volume potentials.Before proving Theorem 6.6, we recall the latter definitions, and for convenience write them in the form (Also, recall from Section 3 that k denote numerical quadrature approximations to these operators over T h with a quadrature rule capable of integrating polynomials of total degree at most m.) Proof of Theorem 6.6.The proof, with exposition primarily restricted to the operator V k and relevant differences for the W k case mentioned as needed, is divided into in two main parts.The first part proceeds by showing that the error incurred by neglecting the weakly-singular integral over K in ( V k is an O(h n+3 |log h|) quantity as h → 0, and then next that the near-singular volume integrals over each of the elements within the neighborhood N h of x are each O(h n+3 |log h|) quantities so that the integral over N h \ K approximates the integral over N h with errors that behave as O(h n+3 |log h|) as h → 0.Then, the second part considers the error from numerical quadrature that arises both from integration on the elements in the near-and far-fields Ω \ N h .Throughout the proof, C j , j = 0, . . ., 6, will denote positive constants independent of h.
Part 1: contributions from N h .Let x ∈ K ∩ E h denote a fixed evaluation point.In what follows the Lagrange interpolation error estimate [71, Thm.1] (that holds since the interpolant here is polynomial and thus easily satisfies the required uniformity condition needed there) will prove useful; it provides (cf.also [23,Thm. 2]) (6.11) |f (y) − f n (y; K)| ≤ C 0 h n+1 for all y ∈ N h (x), and, for vector-valued functions the corresponding estimate (6.12) where C 0 and C 1 denote constants independent of h but possibly dependent on the quotient of the diameter of N h and h.Note that these estimates hold not merely on K but as well over an O(h) neighborhood; indeed, the result of [23, Thm.2] can be applied over the entirety of N h with interpolation conditions enforced at The first step of the proof provides bounds on components of the regularized volume potential.Considering each of the integrals over the elements of N h , we have from a change to polar variables centered at x ∈ K ∩ E h ⊂ N h and the bound (6.11), the estimate (6.13) In a similar vein, it is easy to see that where we used the fact Similarly, since N h \ K ⊂ N h we get from the latter inequality of (6.13) the bound and a similar argument for W k yields (6.17) Turning to numerical quadratures and using the rule (2.8), for the elements comprising N h \ K, using (6.16) and (6.11) in conjunction with the triangle inequality, we find where the inequalities above follow because, firstly, the quadrature rule has positive weights that satisfy (ξ j ,ωj )∈χ:ξ j ∈N h \K ω j ≲ h 2 and secondly, since x ∈ E h lies at a distance from N h \ K that scales linearly with h, sup (ξ j ,ωj )∈χ:ξ j ∈N h \K G k (x, ξ j ) ≲ | log h|, with an implied constant dependent on d T ,E (see Definition 6.2).
In a similar vein for the ( W k kernel, using the triangle inequality and (6.17) together with the estimate (6.12) yields the bound where the last inequality follows because, recalling x ∈ E h , sup (ξ j ,ωj )∈χ:ξ j ∈N h \K ∇ y G k (x, ξ j ) ≲ h −1 with an implied constant dependent on d T ,E .
Part 2: contributions from Ω \ N h .This part relies on a result concerning the accuracy of ordinary quadratures over Ω \ N h , given thereafter in Lemma 6.9.That result in turn requires sharp estimates on the derivatives of the regularized volume integrands (6.20) x, y ∈ Ω, and (6.21) x, y ∈ Ω.
Such estimates for the derivatives of ϕ l (•, y), l = 1, 2, are given in the following lemma whose proof (deferred to Section 6.2) is based on error estimates for multivariate interpolation.
Lemma 6.8.Take as given the assumptions and setting of Theorem 6.6.In particular, n ∈ N 0 is the interpolation degree and m ∈ N, m > n, is the degree of exactness of the quadrature rule.Letting α ∈ N 2 0 be a multi-index satisfying |α| = m + 1 and for ϕ 1 given in (6.20), the following estimate holds for all δ ∈ (0, min{1, (diam Ω) 1/(n−m) }] and any x ∈ K: where C T = C T (m, n, k, f, Ω) and C L = C L (m, n, k, κ T , q T , f, Ω) denote positive constants both independent of x, δ, K, T , h and Λ (Λ defined in Assumption 6.5).
The estimates given in Lemma 6.8 for the derivatives of the integrands ) are on themselves insufficient to obtain convergence rates of higher order than O(h n+1 ) (resp.O(h n )) in connection with classical error estimates [45, Sec.7.4]: while the integrand is globally smooth in Ω \ N h , the growth of derivatives D α y ϕ j (x, •) (j = 1, 2) in the Lemma 6.8 estimates, for derivative orders |α| = m + 1 > n + 1, exactly offsets any possible gain in order of accuracy from use of a higher-order quadrature rule.The following lemma, whose proof is given in Section 6.3, addresses this issue, recovering near-optimal convergence rates (optimal up to a possible | log h| multiplicative factor in some limited cases).It does so by developing sharp bounds for the quadrature error over Ω \ N h .Lemma 6.9.Take as given the assumptions and setting of Theorem 6.6.In particular, n ∈ N 0 is the interpolation degree and m ∈ N, m > n, is the degree of exactness of the reference quadrature rule leading to the global composite quadrature rule Q Ω in (2.7)-(2.8).For all x ∈ E h and for sufficiently small h > 0, it holds that (6.24) with positive constants C j = C j (m, n, k, κ T , q T , f, Ω) (j = 1, 2) independent of h and Λ (Λ defined in Assumption 6.5).
Proof completion.Collecting (6.14), (6.18), and (6.24) from Lemma 6.9, we conclude that for all x ∈ K ∩ E h it holds that Similarly, for ( W k by collecting (6.15), (6.19), and (6.25) from Lemma 6.9 we find that for all x ∈ K ∩ E h it holds that The proof is complete.
Remark 6.10.Weaker versions of Theorem 6.6 can be easily established using a corresponding version of Lemma 6.8 with quadrature order m ≤ n (and, indeed, with no need for Lemma 6.9).However, since the resulting quadrature error of O(h m+1 ) would not saturate the order of accuracy of the regularization error that is established in Theorem 6.6, the use of such low-powered quadratures is disfavored and so not studied in this paper.Remark 6.11.It may be useful for intuition to grasp the basic implication of Lemma 6.9 in a simplified context, ignoring all regularization techniques employed in this paper.The proof of Lemma 6.9 (see Section 6.3) in essence establishes that composite quadrature of log |y| over B 1 (0), when ignoring the integral over the mesh element that contains the origin and when using a rule of exactness degree m ≥ 2, results in errors of O(h 2 | log h|); the (pessimistic) classical error estimate [45,Sec. 7.4] does not predict convergence.The techniques involve use of bounds on concentric rings (see Figure 3) surrounding the singularity over which classical estimates are repeatedly used, implementing the intuition that the integrand is 'well-behaved' on 'most' mesh elements.Such composite quadrature error estimates over unstructured meshes in the presence of an isolated singularity appear to be novel.We find also low-order accurate estimates on quadrature over uniform structured grids in the presence of weak singularities [74] (some of the proof techniques appear similar to those that we independently developed and employed); all higher-order estimates provided there (which are, additionally, restricted to parallelepipeds) make essential use of composite quadratures on graded meshes.We also find interesting connection to a different body of work in the broader context of non-composite quadrature in the presence of isolated singularities, much of it involving Rabinowitz [27,49,60,61,77], for Gauss or Clenshaw-Curtis type quadratures.
6.2.Proof of Lemma 6.8.Our goal is to find suitable bounds for D α y ϕ 1 (x, y), which, applying the product rule, can be expressed as (6.26) We start by noticing that the partial derivatives of G k in (1.2) satisfy the bounds (6.27) that follow from the properties of the Hankel/logarithm function and induction [28, §10.6, 10.7].
Our proof strategy to bound the other factors in (6.26) is based on expressing (6.28) x,n [f ](y), where, denoting by T x,n [f ](y) the nth total degree Taylor polynomial of f centered at x, we define the two remainder functions (6.29) The term R x,n [f ] can be simply expressed in terms of the Taylor remainder, while R x,n [f ] will be rewritten using a certain 'multipoint' Taylor formula for Lagrange polynomials [23].
As expected, for η ∈ N 2 0 , the Taylor remainder R x,n [f ](y) satisfies (6.30) This estimate can be proven for |η| < n + 1 by first leveraging the identity D η R (1) x,n−|η| [D η f ] which follows directly from the fact that D η T x,n [f ] is the (n − |η|)th-degree Taylor polynomial of D η f .From the integral form of the Taylor remainder R (For the sake of simplicity of presentation, we have assumed here that x and y can be joined by a straight line, an assumption that is fulfilled for convex domains Ω.However, a simple generalization of the remainder formula [29, Thm.A1], which is valid for path-connected domains, can be used instead.)Then, from the bounds we readily arrive at (6.33) which gives an expression for the constant C R (1) in the case |η| ≤ n + 1.For n + 1 < |η| ≤ m + 1, on the other hand, the bound follows directly from the fact that D η T x,n [f ] = 0, which allows us to select Remembering that α ∈ N 2 0 is a multi-index satisfying |α| = m + 1 > n + 1 and utilizing the bounds (6.30) with η = β in conjunction with (6.27) with η = α − β , we have To bound the term with the logarithm, we first note that log δ T. G. ANDERSON, M. BONNET, L. FARIA, AND C. P ÉREZ-ARANCIBIA for a suitable implied constant.
To bound the other two terms in (6.34), we note that since |x − y| ≥ δ, for m > n it holds that |x − y| n−m ≤ δ n−m .Then, for |α| = m + 1 and for suitable implied constants we have, using the relation and where in the latter the fact that δ ≤ 1 was also used.Therefore, using the above relations we conclude that (6.34) simplifies to (6.35) where C T = C T (m, n, k, f, Ω).This completes the estimate concerning R (1) x,n [f ] we will rewrite it using [23,Thm. 1]; we note that the hypotheses of that theorem are (i) That the interpolation set I n is poised (termed 'k-unisolvent' there) and (ii) That H(K) (defined in Assumption 6.5) is star-shaped with respect to every element in I n (termed 'K is a Σ-admissible set' there).Both are assured by the present hypotheses of Assumption 6.5; note that by construction H(K) is always star-shaped with respect to each interpolation node x i ∈ I n .Thus, proceeding first with estimates for R (2) x,n [f ], the theorem provides for γ ∈ N 2 0 that the derivatives of the Lagrange interpolation error can be expressed as (6.36) with η i (y) = θ i y + (1 − θ i )x i for some 0 < θ i < 1 and x i ∈ I n .Here D (n+1) f (η i (y); x i − y) denotes the (n + 1)th-order total differential of f at η i (y) [50,Ch.7];i.e., (6.37) (In the event that H(K) ̸ ⊂ Ω -that is, if K is a curved triangle with concave boundary and so f is not defined for all possible η i (y) -we identify f here and in what follows with its Sobolev extension to R 2 assured by the continuous linear extension operator E s : H s (Ω) → H s (R 2 ), s ∈ N 0 , provided by [70, §6, Thm.5] for Lipschitz Ω.This, together with the Sobolev lemma [32] accounts for the assumed additional technical need of f ∈ C n+3 (Ω) so that the extension is also an element of C n+1 (H(K)) with norm bounded by ∥E n+3 ∥ times the norm of f ∈ C n+3 (Ω).)Now, since both T x,n [f ](y) and f n (y; K) are polynomials of degree n in y, we can write where the coefficients {c γ } |γ|≤n that measure error in approximating the Taylor coefficients are given by Note that we have used here fact that T x,n [f ] is the nth-degree Taylor polynomial of f at x so that We then obtain from (6.36) that (6.39) We next provide bounds on the coefficients {c γ } |γ|≤n .In view of (6.37) we first observe that (6.40) . In order to bound the derivatives of the Lagrange polynomials λ i , we use the following Markov-like inequality: for the convex hull H(K), K ∈ T h , and a polynomial p of total degree at most n satisfying sup x∈H(K) |p(x)| ≤ Λ, it holds that (6.41) sup (where κ T and q T are the shape-regularity and quasi-uniformity constants of the mesh family).This inequality stems from iterating |γ| times the inequality chain [76, eq. 3.2] sup in which the width ω H(K) of a compact H(K) ⊂ R 2 , defined [76] as the smallest distance between two parallel lines touching ∂H(K) without intersecting its interior, here verifies h/(κ T q T ) ≤ ω H(K) ≤ h on account of the fact that two parallel lines touching ∂H(K) are separated by at least the diameter of the largest inscribed circle in H(K) and are separated by at most the diameter h of the triangle (see also Definition 6.1, from which it follows that h/q T provides a lower bound on the mesh element size in the preceding argument).
Thus, using (6.49) together with (6.50), the estimate (6.47) becomes To obtain a quadrature error estimate over J j=1 A j one simply sums over the rings A j , j = 1, . . ., J: the last sum being estimated above, for m > n + 2 by a constant (as , for m = n + 2 by the harmonic numbers H J that grow as log J ≲ | log h|, and for m = n + 1 by J < δ 0 h −1 . To estimate the quadrature error associated with the integral over A J+1 , we observe that δ 0 < (J +1)h ≤ |y − x| for y ∈ A J+1 .Therefore, using δ = δ 0 in (6.22) and from the assumption h < δ 0 (see beginning of proof) we obtain It hence follows from (6.47) that In view of (6.46), from (6.52) and (6.54) we have that the desired bound in (6.24) is established with The proof for ϕ 2 proceeds in identical fashion.To bound M (2) m+1 (j) for each j = 1, . . ., J, we use the fact that jh ≤ |y − x| for all y ∈ A j , so from (6.23) with δ = jh < δ 0 we obtain Using (6.55) together with (6.50), the estimate (6.47) The quadrature error estimate over J j=1 A j results from summing over the rings A j , j = 1, . . ., J: where we have used once again the boundedness of the series as well as the relation J < δ 0 h −1 together with the logarithmic growth of the harmonic series for m = n + 1.Finally, to estimate the quadrature error over A J+1 we use the fact that |x − y| > (J + 1)h > δ 0 for all y ∈ A J+1 which allows us to use the bound (6.23) with δ = δ 0 .Therefore, from that bound and the assumption h < δ 0 we obtain so from (6.47) in the case l = 2 it follows that This establishes (6.25) with C 2 = C 7 + C 8 , and the proof is thus complete.

Validation examples.
In order to validate the proposed methodology, we manufacture potential evaluations with scalar and vectorial volume densities, f : Ω → C and f := f d : Ω → C 2 (d ∈ R 2 being a constant vector), respectively, that do not involve the evaluation of volume integrals and so are more easily computable as a reference solution for assessing error levels.To do this, we start off by selecting a known smooth function u : Ω → C and define the scalar source density as f := (∆ + k 2 )u : Ω → C. By Green's third identity it readily follows that V k [f ] = v ref in Ω, where the reference potential is given by in terms of the known function, u, and the single-and double-layer Laplace/Helmholtz potentials applied to its Neumann and Dirichlet traces, respectively.Similarly, we have Since both v ref and w ref only entail the evaluation of layer potentials and their gradients, a highly accurate numerical approximation of the boundary integrals yields a reliable expression for the volume potentials that can be used as a reference to measure the numerical errors.In all the numerical examples considered in this section, we employ where |θ| = 1 and z j ∈ Γ, j = 1, . . ., 4, which gives rise to the density function and make the selections k = 2π, k 0 = 3π, d = (1, 1), s = 10, and θ = (cos π 3 , sin π 3 ).The numerical results are demonstrated on six specific domains.We first consider three domains with smooth boundary parametrizations, namely, the unit disk, the kite [24], and the jellyfish [38], whose boundaries are parametrized, respectively, by the smooth curves γ 0 (t) = (cos t, sin t), (7.4a) γ 1 (t) = (cos t + 0.65 cos 2t − 0.65, 1.5 sin t), and (7.4b) These domains encompass curved triangles as described in Section 2.1.1.Finally, our last three domains, namely, the windmill, the Nazca bird, and the snowflake, are polygonal.The latter encompasses a large number of small-scale features, its meshes are quite heterogeneous (max K∈T h h K / min K∈T h h K > 10), and it includes several poor-quality triangles in some cases.
In the numerical examples that follow, the relative numerical errors are measured by means of the formulae which provide approximate relative errors in the natural C 0 (Ω)-norm, where {ξ j } N j=1 are the volumetric quadrature nodes and where v approx and w approx denote respectively the approximate potentials V k [f ] and W k [f ].The evaluation of each operator is performed in the experiments that follow using Vioreanu-Rokhlin quadrature nodes for integration and interpolation.The results are presented in Figure 4 and Figure 5 for the smooth and polygonal domains, respectively, showcasing the obtained errors, Error V and Error W , for meshes of varying sizes h and interpolation degrees n ∈ {0, 1, 2, 3, 4}.We limit our focus to these interpolation degrees for two primary reasons.On one hand, the method's efficiency decreases as the interpolation degree increases, and on the other, to manage the error-amplifying effect of the translation formula (5.7), discussed in detail below in Section 8. Notably, we observe clear convergence orders for interpolation degrees n ∈ {0, 2, 3, 4}, that match the error estimates established in Theorem 6.6.However, interestingly, for the V k operator, we observe evidence of super-convergence in the case n = 1 even as the expected rate is observed for the W k operator (see Remark 6.7).

Timing results.
The method has been implemented in Julia, in the open-source Inti.jl [4] package also developed by the authors.Some single-core double-precision performance statistics of the method for computing V k are given in Tables 2 and 3 for a volume potential computed over a unit circle; all times reported are in seconds.The code was executed in Julia v1.10 on a Linux machine with an AMD Ryzen 7950X3D CPU.The column marked T B Γ,close lists the precomputation cost of setting up near-singular operator evaluation B Γ [φ, ψ].The column T QΩ,BΓ,Pα shows the cost of volume forward maps, boundary forward maps and particular solution at all points in the domain; T C denotes the cost of (i) Factorization of the Vandermonde-like matrix V, of (ii) Computing translation and scaling quantities c K , r K , and PV −1 , and finally of (iii) Construction of the sparse correction matrix.The column T pre displays the cost of all precomputations.The reported times do not make use of vectorized (or batched) use of the FMM, wherein the same FMM call works on multiple densities in a single pass.Turning to online costs, T corr denotes the time required for application of the sparse correction matrix while T app denotes the full cost of applying V k during the online phase.The online throughput is N/T app ; as evidenced by the %fmm column the method captures the full speed of the FMM at all orders.interacting with a scatterer Ω ⊂ R 2 with the examples thus showing how the high-order accurate operator evaluation previously demonstrated can lead to a high-order volume integral equation solver.We consider a scatterer that has a smooth refractive index function η > 0 defined within its boundaries.Outside of Ω, we assume a constant refractive index of η = 1.The scattering problem seeks a volumetric total field u :  1.3e+1 6.3e+1 5.5e−1 7.7e+1 6.2e+1 1.5e−3 3.6e+0 99.9 5.4e+4 1.25e−2 700545 2.7e+1 2.6e+2 1.9e+0 2.9e+2 2.5e+2 7.0e−3 1.3e+1 99.9 5.3e+4 6.25e−3 2796030 5.9e+1 9.4e+2 6.5e+0 1.0e+3 9.8e+2 3.5e−2 5.2e+1 99.9 5.4e+4 3.13e−3 11160315 1.2e+2 3.7e+3 2.5e+1 3.9e+3 3.9e+3 1.4e−1 2.0e+2 99.9 5.5e+4 with the scattered field u sc := u−u inc satisfying Sommerfeld's radiation condition.It can also be reformulated as the Lippmann-Schwinger volume integral equation for u inside the scatterer: with the solution outside the scatterer given by the representation formula The Lippmann-Schwinger equation (7.7) can be effectively solved using Nyström methods based on the high-order discretization approaches presented in this work using iterative methods, wherein the resulting N × N linear system for the approximate values of the total field at the quadrature nodes {ξ j } N j=1 is solved via GMRES [62] and the volume integral operator V k is used as a forward map.This methodology, which only necessitates repeated applications of the operator V k over the fixed domain Ω, allows us to exploit the algorithm's potential fully, as all source-independent steps can be precomputed and reused at each operator application (see Table 1).To assess the accuracy of such a Nyström method for the Lippmann-Schwinger equation, we first focus on a specific scenario where η maintains a constant value of 4 within Ω.In this case, the total field can be determined by solving the transmission problem (7.9) ∆u + ηk 2 u = 0 in Ω, that can itself be recast into a second-kind boundary integral equation [48].By solving the boundary integral equation (and employing an exponentially convergent Nyström method, based on the Martensen-Kussmaul quadrature rule [24,Sec. 3.5] for kernels with logarithmic singularities to do so), we obtain a highly accurate reference solution that serves as the benchmark to quantify the error in the solution of the volume integral equation for this piecewise constant refractive index medium.Figure 6a displays the real part of the total field u in this example as well as the pointwise errors within Ω in the volume integral equation solutions obtained utilizing the proposed method for interpolation degrees n = 1, 2, 3 and 4 for a fixed mesh with h ≈ 0.09.The resulting linear system was solved via GMRES with a tolerance of 10 −7 , which was achieved in approximately the same number of iterations (∼ 340) in all the examples presented in that figure.We consider, finally, a challenging scattering problem involving a piecewise-smooth refractive index η generated from an x-ray image of a human skull, seen in the left panel of Figure 6b.The right panel of that figure shows the real part of the total field resulting from the scattering of a Gaussian beam that impinges on the skull obstacle from the top-left corner.The Lippmann-Schwinger integral equation solution displayed in that figure, was obtained using the proposed method with n = 2.The scatterer diameter is about 9λ, where the wavelength is λ = 2π/k = 2/3 in this example.GMRES convergence was achieved after 753 iterations for a relative error tolerance of 10 −7 .In order to generate a reference solution to roughly estimate the error, the same problem was solved over the same mesh (h = 0.11) but with n = 4.The relative error, measured on a rectangular curve enclosing the obstacle, is approximately 10 −4 .
8. Limitations and Future Work.The scheme developed in this paper is not meant to be directly applied as-is on arbitrarily-large domains or for sources with arbitrarily-fine features leading to highly anisotropic meshes; this section explores the causes of these limitations and proposes means of managing them.Unfortunately, the ill-conditioning issues that arise when directly solving the Vandermonde system on the physical element K, which were in principle avoided by using the translation and scaling strategy, manifest in a different form as numerical instabilities when evaluating (5.6) and the translation formula (5.7).The formulas in (5.6) may encounter classical cancellation errors typical in polynomial evaluation, and, additionally, significant numerical errors can arise in computation of the coefficients {c α [f ](K)} |α|≤n arising from the translation formula in (5.7).The latter may indeed amplify the errors in the coefficients { c β [f ](K)} |β|≤n to a considerable extent and amounts to the main limitation of the proposed methodology.
To estimate the extent of this undesired amplification effect we can bound the norm of the linear transformation that maps { c β [f ](K)} |β|≤n to {c α [f ](K)} |α|≤n .In view of (5.7), this mapping is represented by the matrix P K ∈ R qn×qn with coefficients In order to estimate the matrix norm ∥P K ∥ ∞ , we note from the definition of p β−α and the inequality Translating the summation indices by α = (α 1 , α 2 ) and restricting our attention to the relevant case It therefore follows from the bound above that (8.1) This analysis provides a rough indication of the impact of the error amplifying effect of the translation formula on the overall error in the proposed methodology and it reveals a critical quantity on which the error for a given domain discretization depends, namely, max K∈T h |c K |/r K .To examine this in more detail we consider a numerical experiment where Error V is computed for integral operators V k integrating over a sequence of translated square domains Ω ℓ = (−1, 0) × (−1, 0) + t ℓ where t 0 = (0, 0) and t ℓ = 2 ℓ−1 (1, 1) for ℓ = 1, . . ., 5, which are applied to densities f (x − t ℓ ), where f is given by for the parameters s = 100, k = 2π, k 0 = 3π, and θ = (cos π 2 , sin π 3 ).To effectively handle the peaked Gaussian source component in f , which is centered at the square's top-right corner, the mesh T h of Ω 0 is locally refined around (0, 0) resulting in a ratio max K∈T h h K / min K∈T h h K ≈ (4.12 × 10 −2 )/(5.65 × 10 −4 ) = 72.92(the meshed domain is displayed on the right panel in Figure 7).The remaining meshes of the domains Ω ℓ , ℓ = 1, . . ., 5, are constructed by a direct translation of T h , i.e., as T h + t ℓ = {K + t ℓ : K ∈ T h }.The plot in Figure 7 displays Error V for each one of the domains Ω ℓ versus the parameter max K∈t ℓ +T h |c K |/r K that, in view of (8.1), appears to be relevant to assess the impact of translation formula on the error.These results show that for sufficiently large max K∈t ℓ +T h |c K |/r K values, the error in evaluating the interpolation polynomial coefficients is amplified to the extent that it becomes dominant, leading to Error V ∝ max K∈t ℓ +T h (|c K |/r K ) n as expected from the estimate (8.1).This verifies that the translation formulae limits the achievable accuracy of the proposed methodology, especially for interpolation degrees n > 4. Indeed, these results show that even in the most benign case, Ω 0 , where max K∈T h |c K | ≤ 1, for interpolation degrees n > 4 the overall error is dominated by the amplified errors in the coefficients of the interpolation polynomial in the monomial basis, a phenomenon also revealed by the estimate (8.1).
The analysis and experiments presented above suggest a strategy to control the accuracy deterioration due to the use of the translation formula, at least to some extent.Such strategy consists of decomposing the mesh T h as D j=1 T (j) h , and find t j that minimizes max K∈T (j) h |c K − t j |/r K for j = 1, . . ., D. This min-max problem, however, can be difficult to solve in practice, but one may instead choose to use a heuristic for the location of the translation point t j to approximate it as a weighted center of mass of T (j) h that accounts for the 1/r K factor, i.e., tj = K∈T Provided D is sufficiently large (i.e.there are enough subdomains), one can manage the size of max 1≤j≤D max K∈T (j) h |c K − t j |/r K , thus alleviating the error-amplifying effect of the translation formula stemming from the resulting (|c K −t j |/r K ) n factor in (8.1).In the limit case when T (j) h = {K j }, j = 1, . . ., D = L, this domain decomposition procedure reduces to the methodology put forth in [66] which is known not to suffer from the aforementioned instabilities.It is worth noticing, though, that the experiments in Section 7 show that, in practice, for interpolation degrees n ≤ 4, one can achieve Error V ≈ 10 −11 without any rescaling or domain decomposition, for domains exhibiting both max K∈T h |c K |/r K and max K∈T h r −1 K well above 5 • 10 3 .In fact, the error analysis of regularized volume integrals developed here justifies a slightly modified alternative approach where the regularization is done more locally, with the boundary of a regularization region laying on a boundary of the union of elements at a distance sufficiently far from the evaluation point so that those boundary integrals can be treated with regular quadrature, while the volume integrals in the complementary region are sufficiently smooth so as to not require regularization.As a byproduct of the locality, no fast algorithms would be needed in the off-line phase of such a variant.We are currently pursuing this idea, and anticipate that this will allow the positive features of the method to seamlessly hold in settings with extreme mesh anisotropy or in large domains with no additional cost.9. Conclusions.We have presented a provably high-order scheme for the numerical evaluation of volume integral operators that is amenable to fast algorithms.Future work may consider the macro forms of domain subdivision mentioned in Section 8 and their effect on computational efficiency.Applications to other PDEs will be explored in future work and is in some cases quite natural given the generality afforded by [3,31].The reduction in the dimensionality of the region of singular quadrature should prove especially effective in three dimensions.We look forward to more finely tuned and parallelized schemes built on these concepts and coupling them to more adaptive methods.
Appendix A. Invertibility and conditioning of the translated-and-scaled Vandermonde matrix.This appendix addresses the invertibility and conditioning of the multivariate Vandermonde matrix V := V K in (5.1) and (5.2) corresponding to a set of distinct points I n .It is well-known [36,55,64] that there is a subtle (and increasingly so at high order) interplay between the geometry of the nodal set and the solvability of the interpolation problem (nevertheless, this is often still the case [64]).The poisedness for node sets transformed from known-poised sets is, in our case, a consequence of reference [23,Lem. 1], by an argument which is partially used in the proof of Theorem A.2 below.However, until the effect of the transformation on the node-sets is understood, the conditioning could in principle be arbitrarily poor.
To address these matters, drawing in part on results in [25], in Theorem A.2 below we show that the linear system in (5.1) for an arbitrary element K generated by an affine transformation T , is not only necessarily invertible ( I n is poised) but has a condition number that can be bounded above by a function that depends explicitly on a standard measure of triangle quality and the degree n of the polynomial basis.For simplicity, we restrict the theorem statement to straight elements; a similar statement could be developed for curved elements generated by (or approximated by) polynomial transformations.It will be useful to consider the geometrical parameters (see Figure 8) (A.1) h K := diam K and ρ K := sup a>0 2a : B a (x) ⊂ K, x ∈ R 2 , that denote the diameter (length of longest edge) of the triangle and the maximum diameter of inscribed disks B a (x) := {y ∈ R 2 : |x − y| < a} in the triangle, respectively.They, together through their ratio, provide a common measure of triangle quality.Note that by construction we have √ 3 ≤ h K ≤ 2; the parameters h K = √ 3 and ρ K = 1 denote analogous quantities for the reference triangle K.In order to state and prove the result, it will help to recall a few facts from the work [25].First, a given nodal interpolation set Y = {ζ 1 , ζ 2 , . . ., ζ qn } ⊂ B 1 (0) ⊂ R 2 is said to be Λ-poised if and only if the Euclidean norm ∥λ(•)∥ 2 of the vector function λ : R 2 → R qn with each element a member of the Lagrange polynomial basis for that nodal set Y, satisfies sup ζ∈B1(0) ∥λ(ζ)∥ 2 ≤ Λ.The quantity sup ζ∈B1(0) ∥λ(ζ)∥ 2 can be bounded above and below by constants times the Lebesgue constant of common parlance, and so in effect Λ functions as such.Then, using the specific factorial-normalized monomial basis set {p α } |α|≤n , p α (y) = y α /α!, the work [25] develops a bidirectional connection between the condition number of the multivariate Vandermonde matrix V (for that basis and the nodal set Y) and the size of Lagrange polynomials associated with Y, i.e. the quantity Λ for which Y is Λ-poised.
We repeatedly use the following lemma, which is a direct adaptation of [25, Thm.1] to our notation; essentially, it represents one direction of the just mentioned bidirectional connection.
Lemma A.1.If V is nonsingular and V −1 ≤ M , then the set Y is Λ-poised with Λ = √ q n M in the unit disk B 1 (0).Conversely, if the set Y is Λ-poised in the unit disk B 1 (0), then V is nonsingular and where θ n > 0 is dependent on n but independent of Y and Λ. Fig. 9: Dependence of the condition number κ( V) on the interpolation degree n and the triangle's shape characterized by the "aspect ratio" hK /ρK .Each of the data points marked by a colored plus sign (+) corresponds to a triangle K with randomly generated vertices contained in [−10, 10] 2 .A total of 1000 triangles were used in this example, each tested for every n in the study.The dashed lines mark the asymptotic behavior of the bound for κ( V)/κ( V) established in Theorem A.2.A remarkable agreement between the empirical data and theoretical bound is observed.
Theorem A.2.For a given n ∈ N, let the nodes I n be mapped from the Vioreanu-Rokhlin nodes I n of interpolation degree n via a non-degenerate affine transformation T .We have that the interpolation sets I n and I n are poised in the normalized monomial basis {p α } |α|≤n .Further, the associated multivariate Vandermonde matrix V in (5.2) is invertible and its condition number satisfies the bound where χ n is computable and depends only on n, and where V is the multivariate Vandermonde matrix in the normalized basis on I n .
To assess the sharpness of the estimate (A.2), we present Figure 9 which displays the quotient κ( V)/κ( V) plotted against 1 + h K /ρ K for matrices V corresponding to 1000 triangles K with random vertices uniformly distributed within [−10, 10] 2 and for various interpolation degrees.These results suggest that our bound (A.2) effectively captures the dependence of the condition number κ( V) on the triangle's shape characterized here by the "aspect ratio" h K /ρ K .
Remark A.3.Explicit bounds on the "constant" θ n in Lemma A.1 can be computed; reference [25, Eq. ( 13)] provides the necessary ingredient for n = 2, showing a path for bounding θ n for larger n as well.Concretely, while the quantity κ( V) can be computed once and for all in the basis {p α } |α|≤n and is of very small size, the theorem represents a somewhat modest result because χ n in Theorem A.2 depends linearly on θ n in Lemma A.1, and the latter quantity admittedly, as reference [25] notes, likely grows as much as exponentially with n.
Proof.As a first step, we estimate the Λ-poisedness of I n .It follows from Lemma A.1, the identity ∥ V −1 ∥ = κ( V)∥ V∥ −1 , and the fact that ∥ V∥ ≥ 1, that I n is Λ-poised with (A.3) Λ := √ q n κ( V), where κ( V) denotes the condition number of the multivariate Vandermonde matrix for the Vioreanu-Rokhlin nodes in the monomial basis.

Fig. 2 :
Fig. 2: Left: Selection of translation and scale parameters for triangles.Here v1, v2, and v3 denote the vertices of the triangle K which are numbered in such a way that the side lengths |v2 −v1|, |v3 −v2| and |v3 −v1| are in decreasing order.Right: The resulting positions of triangles in the unit ball after translationand-scaling.

Proposition 6 . 3 .
Let T = {T h } h>0 be a quasi-uniform and shape-regular family of straight triangulations of Ω.The family E = {E h } h>0 with E h being a set of interior quadrature nodes over a triangulation T h ⊂ T of Ω, is a well-separated family of evaluation point sets with respect to T .Proof.Let I denote the quadrature nodes on K and set δ = inf y∈∂ K, x∈ I | x − y| > 0, which is clearly independent of the families T and E. Consider a triangle K ∈ T h and let T : K → K be the affine transformation T ( x) = A x + b.Then I = T ( I) equals the set of quadrature nodes contained in K. From the bijectivity of T and [23, Lem.2] we have

Fig. 4 :
Fig. 4: Numerical accuracy in the evaluation of volume potentials V k and W k .Top Row: Plot of the real part of the density functions f , defined in (7.3), that are used in the numerical experiments for each of the three smooth domains Ω. Relative errors in the numerical evaluation of V k (middle row) and W k (bottom row) for various discretization sizes.

7. 3 .Fig. 5 :
Fig. 5: Numerical accuracy in the evaluation of volume potentials V k and W k using the proposed polynomial density interpolation method.Top Row: Plot of the real part of the density functions f , defined in (7.3), that are used in the numerical experiments for each of the three polygonal domains Ω. Middle and Bottom Rows: Relative errors in the numerical evaluation of V k (middle row) and W k (bottom row) for various discretization sizes.
< l a t e x i t s h a 1 _ b a s e 6 4 = " z 3 Y n Y u k w 4 O C 2 M H 0 Q J S r M u v w b d I g = " > A A A B 7 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F h P B K t y l i D Z C Q A v L C F 4 S S I 6 w t 9 l L l u z u H b t 7 Q j j y G 2 w s F L H 1 B 9 n 5 b 9 w k V 2 j i g 4 H H e z P M z A s T z r R x 3 W + n s L G 5 t b 1 T 3 C 3 t 7 R 8 c H p W P T 9 o 6 T h W h P o l 5 r L o h 1 p Q z S X 3 D D K f d R F E s Q k 4 7 4 e R 2 7 n e e q N I s l o 9 m m t B A 4 J F k E S P Y W M m v y h u y 0 n o r Q d g p s x n r V m 4 v / e b 3 U R N d B x m S S G i r J c l G U c m R i N P 8 c D Z m i x P C p J Z g o Z m 9 F Z I w V J s b m U 7 I h e K s v r 5 N 2 v e Y 1 a o 2 H e q V 5 l 8 d R h D M 4 h 0 v w 4 A q a c A 8 t 8 I E A g 2 d 4 h T d H O i / O u / O x b C 0 4 + c w p / I H z + Q O E t I 3 b < / l a t e x i t > n = 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " M d X I I m b L j 9 G l r f x 0 N 8 J 1 I q 1 a E T E = " > A A A B 7 H i c b V A 9 S w N B E J 2 N X z F + R S 1 t F h P B K t y l i D Z C Q A v L C F 4 S S I 6 w t 9 l L l u z t H b t 7 Q j j y G 2 w s F L H 1 B 9 n 5 b 9 w k V 2 j i g 4 H H e z P M z A s S w b V x n G 9 U 2 N j c 2 t 4 p 7 p b 2 9 g 8 O j 8 r H J 2 0 d p 4 o y j 8 Y i V t 2 A a C a 4 Z J 7 h R r B u o h i J A s E 6 w e R 2 7 n e e m N I 8 l o 9 m m j A / Ii P J Q 0 6 J s Z J X l T f 1 6 q B c c W r O A n i d u D m p Q I 7 W o P z V H 8 Y 0 j Z g 0 V B C t e 6 6 T G D 8 j y n A q 2 K z U T z V L C J 2 Q E e t Z K k n E t J 8 t j p 3 h C 6 s M c R g r W 9 L g h f p 7 I i O R 1 t M o s J 0 R M W O 9 6 s 3 F / 7 x e a s J r P + M y S Q 2 T d L k o T A U 2 M Z 5 / j o d c M W r E 1 B J C F b e 3 Y j o m i l B j 8 y n Z E N z V l 9 d J u 1 5 z G 7 X G Q 7 3 S v M v j K M I Z n M M l u H A F T b i H F n h A g c M z v M I b k u g F v a O P Z W s B 5 T O n 8 A f o 8 w e G O Y 3 c < / l a t e x i t > n = 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 T H T h w Z n H v A T P + Q g V Z z r r Y Y B P j E = " > A A A B 7 H i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C b a C p 7 J b o X o R C n r w W M G 1 h X Y p 2 T T b h i b Z J c k K Z e l v 8 O J B E a / + I G / + G 9 N 2 D 9 r 6 Y O D x 3 g w z 8 8 K E M 2 1 c 9 9 s p r K 1 v b G 4 V t 0 s 7 u 3 v 7 B + X D o 0 c d p 4 p Q n 8 Q 8 V p 0 Q a 8 q Z p L 5 h h t N O o i g W I a f t c H w z 8 9 t P V G k W y w c z S W g g 8 F C y i B F s r O R X 5 f V F t V + u u D V 3 D r R K v J x U I E e r X / 7 q D W K S C i o N 4 V j rr u c m J s i w M o x w O i 3 1 U k 0 T T M Z 4 S L u W S i y o D r L 5 s V N 0 Z p U B i m J l S x o 0 V 3 9 P Z F h o P R G h 7 R T Y j P S y N x P / 8 7 q p i a 6 C j M k k N V S S x a I o 5 c j E a P Y 5 G j B F i e E T S z B R z N 6 K y A g r T I z N p 2 R D 8 J Z f X i W P 9 Z r X q D X u 6 5 X m b R 5 H E U 7 g F M 7 B g 0 t o w h 2 0 w A c C D J 7 h F d 4 c 6 b w 4 7 8 7 H o r X g 5 D P H 8 A f O 5 w + H v o 3 d < / l a t e x i t > n = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " U 9 u QA d x A d Y Y X W f O W 6 V K t 8 t c W w r g = " > A A A B 7 H i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C b a C p 7 J b p H o R C n r w W M G 1 h X Y p 2 T T b h i b Z J c k K Z e l v 8 O J B E a / + I G / + G 9 N 2 D 9 r 6 Y O D x 3 g w z 8 8 K E M 2 1 c 9 9 s p r K 1 v b G 4 V t 0 s 7 u 3 v 7 B + X D o 0 c d p 4 p Q n 8 Q 8 V p 0 Q a 8 q Z p L 5 h h t N O o i g W I a f t c H w z 8 9 t P V G k W y w c z S W g g 8 F C y i B F s r O R X 5 f V F t V + u u D V 3 D r R K v J x U I E e r X / 7 q D W K S C i o N 4 V j r r u c m J s i w M o x w O i 3 1 U k 0 T T M Z 4 S L u W S i y o D r L 5 s V N 0 Z p U B i m J l S x o 0 V 3 9 P Z F h o P R G h 7 R T Y j P S y N x P / 8 7 q p i a 6 C j M k k N V S S x a I o 5 c j E a P Y 5 G j B F i e E T S z B R z N 6 K y A g r T I z N p 2 R D 8 J Z f X i W P 9 Z r X q D X u6 5 X m b R 5 H E U 7 g F M 7 B g 0 t o w h 2 0 w A c C D J 7 h F d 4 c 6 b w 4 7 8 7 H o r X g 5 D P H 8 A f O 5 w + J Q 4 3 e < / l a t e x i t > n = 4 (a) Planewave scattering by a obstacle, for k = 2π and η = 4 in Ω. Left: Real part of the total field u solution of the transmission problem (7.9); nearly-singular volume integration arising in near-field evaluation of the representation formula (7.8) is handled seamlessly by VDIM per Remark 3.1 Right: Logarithm in base ten of the absolute pointwise error in the Lippmann-Schwinger equation solution for interpolation degrees n = 1, 2, 3 and 4 over a fixed mesh of Ω with h ≈ 0.09.The maximum poinwise errors in those figures are around 3.0 • 10 −3 , 2.3 • 10 −5 , 8.8 • 10 −6 , and 9.2 • 10 −7 , respectively.(b) Gaussian beam scattering by an inhomogeneous obstacle, for k = 3π and a variable refractive index.Left: Piecewise smooth refractive index η generated from a skull x-ray image.Right: Real part of the total field, resulting from the scattering of a Gaussian beam by the skull-like inhomogeneous obstacle, obtained from an approximate Lippmann-Schwinger equation solution.

Fig. 6 :
Fig. 6: Application of the proposed polynomial density interpolation method to the Lippmann-Schwinger equation.Convergence experiments demonstrate the same high convergence orders for solution of the volume integral equation that are illustrated for the operator evaluation in Figure 4, and to our knowledge represent the first reported provably high-order accurate solutions of such integral equations with piecewise-smooth material parameters.

Fig. 7 :
Fig. 7: Effect of the translation formula (5.7) on the errors in the evaluation of the operator V k .Left: Relative errors in the evaluation of V k integrating over translations of the square domain (−1, 0)×(−1, 0).Right: Locally refined mesh of the square domain used to resolve the peaked density function (8.2) utilized in the corresponding numerical experiments.

Fig. 8 :
Fig.8: Coordinate systems and transformations.The condition number of the Vandermonde-like system for interpolation in the unit ball (lower) can be controlled.

Table 2 :
Performance table for Laplace using FMM (ϵ fmm = 10 −13 ); all experiments are performed on a single core, and h values are chosen to result in a range of N roughly from 10 thousand to 10 million.The offline throughput (N/Tpre) at large N is steady around 23, 000 for n = 2 and around 9200 targets/core/second for n = 4. Memory usage of the volume correction matrix is displayed in MB.

Table 3 :
Performance table for Helmholtz (k = 10π) using FMM (ϵ fmm = 10 −7 ); all experiments are performed on a single core, and h values are chosen to result in a range of N roughly from 10 thousand to 10 million.The offline throughput (N/Tpre) at large N is approximately 9300 targets/core/second for n = 2 and approximately 3000 targets/core/second for n = 4. Memory usage of the volume correction matrix is displayed in MB.