Modified Representations for the Close Evaluation Problem

Abstract

When using boundary integral equation methods, we represent solutions of a linear partial differential equation as layer potentials. It is well-known that the approximation of layer potentials using quadrature rules suffer from poor resolution when evaluated closed to (but not on) the boundary. To address this challenge, we provide modified representations of the problem’s solution. Similar to Gauss’s law used to modify Laplace’s double-layer potential, we use modified representations of Laplace’s single-layer potential and Helmholtz layer potentials that avoid the close evaluation problem. Some techniques have been developed in the context of the representation formula or using interpolation techniques. We provide alternative modified representations of the layer potentials directly (or when only one density is at stake). Several numerical examples illustrate the efficiency of the technique in two and three dimensions.

1. Introduction

One can represent the solution of partial differential boundary-value problems using boundary integral equation methods, which involves integral operators defined on the domain’s boundary called layer potentials. Using layer potentials, the solution can be evaluated anywhere in the domain without restriction to a particular mesh. For that reason boundary integral equations have found broad applications, including in fluid mechanics, electromagnetics, and plasmonics [1,2,3,4,5,6,7,8].

The close evaluation problem refers to the nonuniform error produced by high-order quadrature rules used to discretize layer potentials. This phenomenon arises when computing the solution close to the boundary (i.e., at close evaluation points). It is well understood that this growth in error is due to the fact that the integrands of the layer potentials become increasingly peaked as the evaluation point approaches the boundary (nearly singular behavior), leading in limited cases to an $O\left(1\right)$ error [9].

There exists a plethora of manners to address the close evaluation problem: using extraction methods based on Taylor series expansions [10], regularizing the nearly singular behavior of the integrand and adding corrections [11,12], compensating quadrature rules via interpolation [13], using Quadrature By Expansion related techniques (QBX) [9,14,15,16,17,18,19], using adaptive methods [20], using singularity subtraction techniques and interpolation [21,22,23], or using asymptotic approximations [24,25,26], to name a few. Most techniques rely on either providing corrections to the kernel (related to the fundamental solution of the PDE at stake), or to the density (solution of the boundary integral equation).

In the latter category, it is well-known that Laplace’s double-layer potential can be straightforwardly modified via a density subtraction technique based on Gauss’ law (e.g., [27]). This modification alleviates the close evaluation problem, and provides a better approximation for any given numerical method. However this identity technique is specific to Laplace’s double-potential. Other identities have been derived for other problems, such as for the elastostatic problem [28].

In this paper, we provide modified representations of layer potentials, and we give guidance to address the close evaluation problem in two and three dimensions. In particular, we modify Laplace’s single-layer potential (representing the solution of the exterior Neumann Laplace problem) and Helmholtz layer potentials (in the context of a sound-soft scattering problem). With some given quadrature rule, the resulted modified representations allow us to obtain better approximations compared to standard representations. The proposed modifications are based on subtracting specific solutions (or auxiliary functions) of the PDE at stake. The use of auxiliary functions have been developed in the context of Boundary Regularized Integral Equation Formulation (BRIEF) [29,30,31] to regularize the representation formula on the boundary, or in the context of density interpolation techniques [21,23,32] to regularize layer potentials (generalization of density subtractions). Those techniques commonly consider multiple auxiliary functions, and may require to solve additional problems to find such functions. The proposed work concentrates on regularizing nearly singular integrals using explicitly one analytic auxiliary function, and when representing the solution with layer potentials involving only one density (no representation formula). We provide several examples of auxiliary functions (and compare them), and provide guidelines to find them. The proposed modified representations are simple and easy to implement, and allow one to straightforwardly gain accuracy in evaluating the solution, especially when computational resources are limited. This work provides valuable insights into Laplace and Helmholtz layer potentials. Additionally this can also be applied to modify boundary integral equations to avoid weakly singular integrals.

The paper is organized as follows: Section 2 presents some context and motivation for the proposed modified representations. Section 3 establishes the modified representations and general guidelines to find appropriate auxiliary functions. Section 4 and Section 5 illustrate the efficiency of the modified representations for Laplace and Helmholtz in two and three dimensions, off and on boundary. Finally, Section 6 presents our concluding remarks, Appendices Appendix A and Appendix B provide a brief summary of the Nyström methods used in two and three dimensions, and Appendix C details some proofs for Section 3.

2. Motivation for Modified Representations

Consider a domain $D\subset {\mathbb{R}}^d$, $d=2,3$, that is a bounded simply connected open set with smooth boundary (of class ${\mathcal{C}}^2$), and a linear elliptic partial differential equation of the form $\mathcal{L}u=0$. It is common to represent the solution $\mathtt{v}$ of that PDE using the so-called representation formula (e.g., Theorem 6.5 in [33], Theorem 3.1 in [34]). In particular for $\mathtt{v}$ satisfying $\mathcal{L}\mathtt{v}=0$ in D, we have the following identities:

$$\begin{array}{cc}\hfill & {\displaystyle {\int }_{\partial D}{\partial }_{n_y}G(x,y)\mathtt{v}\left(y\right)d{\sigma }_y-{\int }_{\partial D}G(x,y){\partial }_{n_y}\mathtt{v}\left(y\right)d{\sigma }_y=\left\{\begin{array}{cc}-\mathtt{v}\left(x\right)\hfill & x\in D,\hfill \\ {\displaystyle -\frac{1}{2}\mathtt{v}\left(x\right)}\hfill & x\in \partial D,\hfill \\ 0\hfill & x\in E:={\mathbb{R}}^d\backslash \overline{D},\hfill \end{array}\right.}\hfill \end{array}$$

(1)

where G denotes the fundamental solution of considered PDE, $n_y$ is the unit outward normal of D at y, and $d{\sigma }_y$ is the integration surface element. For instance, (1) holds true for $\mathcal{L}:=\Delta$ and $\mathcal{L}:=\Delta +k^2$, the Laplace and the Helmholtz equation, respectively. The goal of this paper is to use (1) with well-chosen $\mathtt{v}$ to modify the representation of the solution of boundary value problems associated with $\mathcal{L}$. Let us illustrate the strategy with, for example, the Exterior Neumann Laplace problem:

$$\left|\begin{array}{c}\mathrm{Find} u\in {\mathcal{C}}^2\left(E\right)\cap {\mathcal{C}}^1(\overline{E}:={\mathbb{R}}^d\backslash D)\mathrm{such} \mathrm{that}:\hfill \\ \Delta u=0\mathrm{in}E,{\partial }_{n}u=g\mathrm{on}\partial D,\underset{\left|x\right|\to \infty }{lim}u\left(x\right)=o\left(1\right),\hfill \end{array}\right.$$

(2)

with some smooth data g (with null average). The solution of Problem (2) can be represented using Green’s formula [34,35]:

$$\begin{array}{cc}\hfill u\left(x\right)& {\displaystyle ={\int }_{\partial D}{\partial }_{n_y}G(x,y)u\left(y\right)d{\sigma }_y-{\int }_{\partial D}G(x,y){\partial }_{n_y}u\left(y\right)d{\sigma }_y,x\in E,}\hfill \\ \hfill & {\displaystyle ={\int }_{\partial D}{\partial }_{n_y}G(x,y)u\left(y\right)d{\sigma }_y-{\int }_{\partial D}G(x,y)g\left(y\right)d{\sigma }_y,x\in E,}\hfill \end{array}$$

(3)

where

$$G(x,y)=\left\{\begin{array}{cc}& {\displaystyle -\frac{1}{2\pi }log|x-y|\mathrm{for} d=2,}\hfill \\ \hfill & {\displaystyle \frac{1}{4\pi }\frac{1}{|x-y|}\mathrm{for} d=3,}\hfill \end{array}\right.$$

(4)

and the trace on the boundary satisfies the boundary integral equation of the second kind:

$${\displaystyle \frac{1}{2}u\left(x^{*}\right)-{\int }_{\partial D}{\partial }_{n_y}G(x^{*},y)u\left(y\right)d{\sigma }_y={\int }_{\partial D}G(x^{*},y)g\left(y\right)d{\sigma }_y,x^{*}\in \partial D.}$$

(5)

The fundamental solution G is singular when $y=x^{*}$. For $x\in {\mathbb{R}}^d\backslash \partial D$, assume we can write $x=x^{*}\pm \ell n_{x^{*}}$ with $n_{x^{*}}$ the unit outward normal at $x^{*}$, and $\ell >0$ the distance from the boundary. Then G is nearly singular at $y=x^{*}$ when $|x-y|=\ell \ll 1$ (i.e., when x is close to the boundary). A layer potential is said to be a weakly singular integral (resp. a nearly singular integral) when its kernel (G or ${\partial }_{n}G$ in the cases above) is singular at $y=x^{*}$ (resp. nearly singular at $y=x^{*}$). There exist high-order quadrature rules to approximate weakly singular integrals with very high accuracy (e.g., [36,37,38,39]). However, high accuracy is lost for nearly singular integrals: this is the so-called close evaluation problem. Assuming we have solved (5), we can modify (3) using (1) to address the close evaluation problem. Taking the difference we obtain

$$\begin{array}{cc}\hfill u\left(x\right)& {\displaystyle ={\int }_{\partial D}{\partial }_{n_y}G(x,y)[u\left(y\right)-\mathtt{v}\left(y\right)]d{\sigma }_y-{\int }_{\partial D}G(x,y)[g\left(y\right)-{\partial }_{n_y}\mathtt{v}\left(y\right)]d{\sigma }_y,x\in E.}\hfill \end{array}$$

(6)

If one finds $\mathtt{v}$ such that $\mathtt{v}\left(x^{*}\right)=u\left(x^{*}\right)$ and ${\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)=g\left(x^{*}\right)$, where $x^{*}\in \partial D$ denotes the closest boundary point of the evaluation point x ($x=x^{*}+\ell n_{x^{*}}$), then (6) does not suffer from the close evaluation problem.

Similarly, one can represent the solution of Problem (2) using a single-density representation given by the single-layer potential:

$${\displaystyle u\left(x\right)={\int }_{\partial D}G(x,y)\rho \left(y\right)d{\sigma }_y,x\in D,}$$

(7)

with $\rho$ a continuous density solution of the boundary integral equation of the second-kind:

$${\displaystyle -\frac{1}{2}\rho \left(x^{*}\right)+{\int }_{\partial D}{\partial }_{n_x^{*}}G(x^{*},y)\rho \left(y\right)d{\sigma }_y=g\left(x^{*}\right),x^{*}\in \partial D.}$$

(8)

Assuming we have solved (8) for $\rho$, subtracting (1) from (7) we obtain

$$\begin{array}{cc}\hfill u\left(x\right)& {\displaystyle ={\int }_{\partial D}G(x,y)[\rho \left(y\right)-{\partial }_{n_y}\mathtt{v}\left(y\right)]d{\sigma }_y+{\int }_{\partial D}{\partial }_{n_y}G(x,y)\mathtt{v}\left(y\right)d{\sigma }_y,x\in E.}\hfill \end{array}$$

(9)

If one finds $\mathtt{v}$ such that $\mathtt{v}\left(x^{*}\right)=0$ and ${\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)=\rho \left(x^{*}\right)$, then (9) does not suffer from the close evaluation problem.

Representations (6) and (9) are attractive representations, and several works have provided guidelines on how to build appropriate solutions $\mathtt{v}$. For (6) one can use Taylor-like functions $\mathtt{v}\left(x\right)=u\left(x^{*}\right)\tilde{g}\left(x\right)+{\partial }_{n_{x^{*}}}u\left(x^{*}\right)\tilde{f}\left(x\right)$, with $\tilde{g}$ and $\tilde{f}$ solutions of some Laplace boundary value problems [29,30,31]. This technique has been first developed in the context of Boundary Regularized Integral Equation Formulation (BRIEF) (namely to solve (5) using the same subtraction technique on boundary) and applied to evaluate the solution near the boundary. For (9) one can use density interpolation methods [21,23,32]: $\mathtt{v}=\mathtt{v}(x^{*},y)={\sum }_{j=0}^J{c}_j\left(y\right)H_j(x^{*}-y)$ where ${\left(H_j\right)}_j$ satisfy the PDE (in the above case ${\left(H_j\right)}_j$ are harmonic functions). In both methods the chosen auxiliary functions $\mathtt{v}$ necessarily depend on the trace u (and/or normal trace ${\partial }_{n}u$), or the density $\rho$ at the closest evaluation point. Furthermore they require to satisfy at least two conditions (two boundary value problems or two boundary conditions).

In this paper, we provide another construction of modified representations for single-density representations of Laplace and Helmholtz boundary value problems. The construction relies on auxiliary functions $\mathtt{v}$ that are independent of the density (solution of the boundary integral equation), and requires fewer constraints in the context of (7). As a consequence, our approach provides more freedom in choosing $\mathtt{v}$. The proposed modified representations are also simple to implement and do not add significant computational costs. In what follows we provide modified representations for Laplace and Helmholtz in 2D and 3D, and provide several examples to illustrate the efficiency of the method.

3. Modified Representations

We present modified representations for single-density representations of Laplace and Helmholtz boundary value problems. In particular, we consider the interior Dirichlet Laplace problem (where one can represent the solution using the double-layer potential), the exterior Neuman Laplace problem (2) (using the single-layer potential (7)), and the sound-soft scattering problem.

3.1. Modified Representation for the Laplace Double-Layer Potential

The interior Dirichlet problem for Laplace consists of finding $u\in {\mathcal{C}}^2\left(D\right)\cap {\mathcal{C}}^1\left(\overline{D}\right)$ such that

$$\left|\begin{array}{c}\Delta u=0\mathrm{in}D,u=f\mathrm{on}\partial D,\hfill \end{array}\right.$$

(10)

with some smooth data f. The solution of Problem (10) can be represented as a double-layer potential [34,35]:

$${\displaystyle u\left(x\right)={\int }_{\partial D}{\partial }_{n_y}G(x,y)\mu \left(y\right)d{\sigma }_y,x\in D,}$$

(11)

with G defined in (4), and $\mu$ a continuous density solution of the boundary integral equation:

$${\displaystyle -\frac{1}{2}\mu \left(x^{*}\right)+{\int }_{\partial D}{\partial }_{n_y}G(x^{*},y)\mu \left(y\right)d{\sigma }_y=f\left(x^{*}\right),x^{*}\in \partial D.}$$

(12)

We now make use of (1) to modify (11). One can show the following (see Appendix C.1 for details):

Proposition 1.

Given $x=x^{*}-\ell n_{x^{*}}\in D$ with $x^{*}\in \partial D$, let $\mathtt{v}$ be a solution of Laplace’s equation in $D\subset {\mathbb{R}}^d$, $d=2,3$, such that

$$\mathtt{v}\left(x^{*}\right)=1,{\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)=0.$$

(13)

The solution of the exterior Dirichlet Laplace problem (11) admits the modified representation:

$$\begin{array}{cc}& {\displaystyle u\left(x\right)={\int }_{\partial D}{\partial }_{n_y}G(x,y)\mu \left(y\right)\left[1-\mathtt{v}\left(y\right)\right]d{\sigma }_y+{\int }_{\partial D}{\partial }_{n_y}G(x,y)\left[\mu \left(y\right)-\right.\left.\mu \left(x^{*}\right)\right]\mathtt{v}\left(y\right)d{\sigma }_y}\hfill \\ & {\displaystyle -\mu \left(x^{*}\right)\mathtt{v}\left(x^{*}\right)+\mu \left(x^{*}\right){\int }_{\partial D}G(x,y)\left[{\partial }_{n_y}\mathtt{v}\left(y\right)-{\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)\right]d{\sigma }_y-\mu \left(x^{*}\right){\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right),x\in D.}\hfill \end{array}$$

(14)

The modified representation (14) has smoother integrands than (11), and it addresses the close evaluation problem, in the sense that nearly singular terms vanish as $y\to x^{*}$.

From Proposition 1 we can now build auxiliary functions $\mathtt{v}$ independent of $\mu$, and there exist plenty of candidates: constant, linear, based on Green’s function ($\mathtt{v}\left(y\right)=G(y,x_0)$ with $x_0\in E$), quadratic ($\mathtt{v}(y_1,y_2)=1+(y_1-x_1^{*})(y_2-x_2^{*})$, $\mathtt{v}(y_1,y_2)=1+{(y_1-x_1^{*})}^2-{(y_2-x_2^{*})}^2$), $\mathtt{v}(y_1,y_2,y_3)=e^{y_3}(sin{y}_1+sin{y}_2)$, etc. The solution $\mathtt{v}\equiv 1$ naturally satisfies the conditions (13), and the modified representation (14) boils down to

$${\displaystyle u\left(x\right)={\int }_{\partial D}{\partial }_{n_y}G(x,y)[\mu \left(y\right)-\mu \left(x^{*}\right)]d{\sigma }_y-\mu \left(x^{*}\right),x\in D.}$$

(15)

The modified representation (15) is well-known and widely used (e.g., [9,25,27]), it is the simplest representation that naturally addresses the close evaluation problem. Thus, we do not provide numerical results for this case. Rather, we concentrate on other layer potentials.

3.2. Modified Representation for the Laplace Single-Layer Potential

Going back to Problem (2), one can show the following (see Appendix C.2 for details):

Proposition 2.

Given $x=x^{*}+\ell n_{x^{*}}\in E$ with $x^{*}\in \partial D$, let $\mathtt{v}$ be a solution of Laplace’s equation in $D\subset {\mathbb{R}}^d$, $d=2,3$, such that

$${\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)=1.$$

(16)

The solution of the exterior Neumann Laplace problem (2) admits the modified representation:

$$\begin{array}{cc}\hfill {\displaystyle u\left(x\right)={\int }_{\partial D}G(x,y)\rho \left(y\right)}& {\displaystyle \left[1-{\partial }_{n_y}\mathtt{v}\left(y\right)\right]d{\sigma }_y+{\int }_{\partial D}G(x,y)\left[\rho \left(y\right)-\rho \left(x^{*}\right)\right]{\partial }_{n_y}\mathtt{v}\left(y\right)d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +\rho \left(x^{*}\right){\int }_{\partial D}{\partial }_{n_y}G(x,y)\rho \left(y\right)\left[\mathtt{v}\left(y\right)-\mathtt{v}\left(x^{*}\right)\right]d{\sigma }_y,\forall x\in E.}\hfill \end{array}$$

(17)

The modified representation (17) has smoother integrands than (7).

Contrary to auxiliary functions provided in Taylor-like methods and density interpolation methods (discussed in Section 2), auxiliary functions $\mathtt{v}$ do not depend on $\rho$ and rely on only one constraint (16). Therefore, there is a lot of freedom in choosing $\mathtt{v}$: given $\mathtt{u}$ a solution of Laplace’s equation, then one chooses $\mathtt{v}:=\frac{\mathtt{u}}{{\partial }_{n_x^{*}}\mathtt{u}\left(x^{*}\right)}$ (as long as ${\partial }_{n_x^{*}}\mathtt{u}\left(x^{*}\right)\ne 0$). Candidates may then include:

• The linear function $\mathtt{v}\left(y\right)=n_{x^{*}}·y$;
• The function $\mathtt{v}\left(y\right)=2^{d-1}\pi G(y,x^{*}+n_{x^{*}})$ based on Green’s function;
• The quadratic product function ${\displaystyle \mathtt{v}\left(y\right)=\frac{(y_1-x_{0,1})(y_2-x_{0,2})}{n_{x^{*},1}(x_2^{*}-x_{0,2})+n_{x^{*},2}(x_1^{*}-x_{0,1})}}$, $x_0\in D$;
• The quadratic difference function ${\displaystyle \mathtt{v}\left(y\right)=\frac{1}{2}\frac{{(y_1-x_{0,1})}^2-{(y_2-x_{0,2})}^2}{n_{x^{*},1}(x_1^{*}-x_{0,1})-n_{x^{*},2}(x_2^{*}-x_{0,2})}}$, $x_0\in D.$

Note that the above candidates are valid in ${\mathbb{R}}^d$, one can also consider any of the quadratic functions above in ${\mathbb{R}}^3$ as a function of $(y_i,y_j)$, $i,j=1,2,3$, $j\ne i$. In Section 4, we will test (17) using several candidates $\mathtt{v}$ and make comparisons. The modified representation (17) adds two terms to compute compared to (7), it is the price to pay to gain accuracy at close evaluation points. We will make comparative tests to quantify this aspect.

3.3. Modified Representation for the Helmholtz Double- and Single-Layer Potentials

We consider in this case the sound-soft scattering problem:

$$\left|\begin{array}{c}\mathrm{Find}u\in {\mathcal{C}}^2\left(E\right)\cup {\mathcal{C}}^1\left(\overline{E}\right)\mathrm{such}\mathrm{that}:\hfill \\ {\displaystyle \Delta u+k^{2}u=0\mathrm{in}E,u=f\mathrm{on}\partial D,\underset{R\to \infty }{lim}{\int }_{\left|y\right|=R}{|{\partial }_{n}u-iku|}^{2}d{\sigma }_y=0,}\hfill \end{array}\right.$$

(18)

with some smooth data f associated with the wavenumber k. Above, the last condition represents the Sommerfeld radiation condition. The solution of Problem (18) can be represented as a combination of double- and single-layer potentials [40]:

$${\displaystyle u\left(x\right)={\int }_{\partial D}\left[{\partial }_{n_y}{G}^H(x,y)-ik{G}^H(x,y)\right]\mu \left(y\right)d{\sigma }_y,x\in E,}$$

(19)

with $G^H$ defined by

$$G^H(x,y)=\left\{\begin{array}{cc}& {\displaystyle \frac{i}{4}{H}_0^{\left(1\right)}\left(k\right|x-y\left|\right),\mathrm{for} d=2,}\hfill \\ \hfill & {\displaystyle \frac{1}{4\pi }\frac{e^{ik|x-y|}}{|x-y|},\mathrm{for} d=3,}\hfill \end{array}\right.$$

(20)

with $H_0^{\left(1\right)}(·)$ the Hankel function of the first kind, and $\mu$ a continuous density satisfying:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}\mu \left(x^{*}\right)+{\int }_{\partial D}\left[{\partial }_{n_y}\right.}& G^H(x^{*},y)\left.-ik{G}^H(x^{*},y)\right]\mu \left(y\right)d{\sigma }_y=f\left(x^{*}\right),x^{*}\in \partial D.\hfill \end{array}$$

(21)

One obtain the following:

Proposition 3.

Given $x=x^{*}+\ell n_{x^{*}}\in E$ with $x^{*}\in \partial D$, let $\mathtt{v}$ be a solution of Helmholtz equation in $D\subset {\mathbb{R}}^d$, $d=2,3$, such that

$$\mathtt{v}\left(x^{*}\right)=1,{\partial }_{n_{x^{*}}}\mathtt{v}\left(x^{*}\right)=ik.$$

(22)

Then the solution of the sound-soft scattering problem (18) admits the modified representation:

$$\begin{array}{cc}\hfill {\displaystyle u\left(x\right)=}& {\int }_{\partial D}\left[{\partial }_{n_y}{G}^H(x,y)-{\partial }_{n_y}\mathtt{v}\left(y\right)G^H(x,y)\right]\left[\mu \left(y\right)-\mu \left(x^{*}\right)\right]d{\sigma }_y\hfill \\ \hfill & {\displaystyle +{\int }_{\partial D}{G}^H(x,y)\left[{\partial }_{n_y}\mathtt{v}\left(y\right)-ik\right]\mu \left(y\right)d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +\mu \left(x^{*}\right){\int }_{\partial D}{\partial }_{n_y}{G}^H(x,y)\left[1-\mathtt{v}\left(y\right)\right]d{\sigma }_y,\forall x\in E.}\hfill \end{array}$$

(23)

The modified representation (23) has smoother integrands than (19).

The proof can be found in see Appendix C.3. One can check in particular that plane waves $\mathtt{v}\left(y\right)=e^{ik{n}_{x^{*}}·(y-x^{*})}$ do satisfy (22), whereas Green-based functions like $\mathtt{v}\left(y\right)=G^H(y,x_{*}+n_{x^{*}})$ (up to some constant) cannot. We will use (23) with plane waves for the numerical examples.

4. Numerical Examples

The accuracy in approximating (11)–(15), (7)–(17), (19)–(23), respectively, relies on the resolution of the boundary integral Equations (12), (8) and (21), respectively. In what follows we assume that the boundary integral equations are sufficiently resolved. Given the density’s resolution, we compare the representations and their modified ones through several examples. All the codes can be found in [41].

4.1. Exterior Neumann Laplace Problem

4.1.1. Example 1: Exterior Laplace in Two Dimensions

Since $\partial D$ is a closed smooth boundary, we use the Periodic Trapezoid Rule (PTR) to approximate (7) and (17), where we will use several $\mathtt{v}$ according to Proposition 2. We consider an exact solution of Problem (2):

$$\begin{array}{cc}\hfill u_{\mathrm{exact}}\left(x\right)=u_{\mathrm{exact}}(x_1,x_2)& =\frac{x_1-x_{0,1}}{|x-x_0{|}^2},x_0=(x_{0,1},x_{0,2})\in D,\hfill \end{array}$$

which consists of choosing $g\left(x^{*}\right)={\partial }_{n_{x^{*}}}{u}_{\mathrm{exact}}\left(x^{*}\right)$, for any $x^{*}\in \partial D$.

All simulations are done outside of a kite-shaped domain using the Periodic Trapezoid Rule with $N=128$ quadrature points for the following representations:

• V0: standard representation (7);
• V1: modified representation (17) with the linear function ${\mathtt{v}}_1\left(y\right)=n_{x^{*}}·y;$
• V2: modified representation (17) with the Green’s function ${\mathtt{v}}_2\left(y\right)=2\pi G(y,x^{*}+n^{*});$
• V3: modified representation (17) with the quadratic function ${\displaystyle {\mathtt{v}}_3\left(y\right)=\frac{1}{2}\frac{y_1^2-y_2^2}{n_{x^{*},1}{x}_1^{*}-n_{x^{*},2}{x}_2^{*}};}$
• V4: modified representation (17) with the quadratic function

  $${\displaystyle {\mathtt{v}}_4\left(y\right)=\frac{(y_1-5)(y_2-5)}{n_{x^{*},1}(x_2^{*}-5)+n_{x^{*},2}(x_1^{*}-5)}.}$$

We solved (8) using the Nyström method based on the Periodic Trapezoid Rule (using Matlab classic backslash). The accuracy of all methods is limited by the accuracy of the resolution for $\rho$ (in particular when considering moderate N). This can be assessed by looking the density’s Fourier coefficients decay: in this case the coefficients decay is bounded by ${10}^{-5}$ for $N=128$. The results in Figure 1 and Figure 2 show that given $\rho$ resolved, the approximation of the modified representations provide better results overall. Far from the boundary, all methods approximate well the solution. As the evaluation point gets closer to the boundary ($\ell \to 0$), V0 approximated by PTR suffers from the close evaluation problem and the error increases (see [9]). Note that the single-layer potential commonly suffers less from this phenomenon than the double-layer potential (e.g., [24]). Using the modified representations (V1–V4) allows to reduce the error by a couple of orders of magnitude for the close evaluation problem. All modified representations provide a satisfactory correction overall. We use a naïve (straightforward) implementation of (7) and (17) in Matlab, computed on a Mac mini SSD 512Go. We provide run times in

Table 1. Laplace 2D single-layer. CPU times (in seconds) for various number of quadrature points and representations. Times account for computing the solution at $N\times 12$ grid points ($\ell ={10}^{-k}$, $k=〚0,11〛$) on a body-fitted grid.

Method	V0	V1	V2	V3	V4
$N=128$	0.014	0.044	0.055	0.045	0.05
$N=256$	0.056	0.07	0.112	0.08	0.081
$N=512$	0.12	0.192	0.263	0.2	0.19

for various number of quadrature points. Run times do not count the time to compute the boundary integral equation for $\rho$ (being the same for all methods).

Representation V0 is obviously cheaper (less terms to compute) than V1–V4, and V1 is cheaper than V2–V4 due to simpler terms: there are less operations to conduct to compute ${\mathtt{v}}_1\left(y\right)$ than the other provided auxiliary functions.

To better compare the methods, Figure 3 represents log plots of the maximum error with respect to the number of quadrature points N and for various distances ℓ from point A (indicated in Figure 1). The results show that modified representations allow to gain a couple of order of magnitude even for moderate N ($N<100$). Additionally, the error using V0 decreases linearly with the number of quadrature points whereas it is cubic using modified representations. While there is no significant difference between the considered modified representations V1–V4, one may consider run times (and simplicity of auxiliary function $\mathtt{v}$) to discuss competitiveness. Based on the above results, overall representation V1 seems to be the best choice for the best computational cost-accuracy trade-off. Let us emphasize that the focus of this paper is to highlight the efficacy and simplicity of the proposed modified representations, given a quadrature rule. Our results show that modified representations allow to naturally gain a couple of orders of magnitude in the error, addressing the close evaluation problem even for moderate computational resources. Additionally, the proposed auxiliary functions are independent of the density $\rho$. In the next section we investigate the efficacy of (17) in three dimensions.

4.1.2. Example 2: Exterior Laplace in Three Dimensions

Given a domain $D\subset {\mathbb{R}}^3$ with smooth boundary, we assume $\partial D$ to be an analytic, closed, and oriented surface that can be parameterized by $y=y(s,t)$ for $s\in [0,\pi ]$ and $t\in [-\pi ,\pi ]$. Then one can write (7) as

$$u\left(x\right)={\int }_{-\pi }^{\pi }{\int }_0^{\pi }G(x,y(s,t))J(s,t)\rho \left(y(s,t)\right)sin\left(s\right)\mathrm{d}s\mathrm{d}t,$$

(24)

with $J(s,t)=|y_s(s,t)\times y_t(s,t)|/sin\left(s\right)$ the Jacobian. We now work with a surface integral defined on a sphere, and we use a three-step method (see [26] for details) to approximate (7) and (17). This method corresponds to a modification of the product Gaussian quadrature rule (PGQ) [42], and it has been shown to be very effective for computing layer potentials in three dimensions at close evaluation points compared to other quadrature methods for nearly singular integrals [26]. It relies on (i) rotating the local coordinate system so that $x^{*}$ corresponds to the north pole, (ii) use Periodic Trapezoid Rule with $2N$ quadrature points to approximate the integral with respect to t, (iii) use Gauss–Legendre with N quadrature points mapped to $(0,\pi )$ (and not $(-1,1)$) to approximate the integral with respect to s. This leads to the approximation:

$$u\left(x\right)\approx \frac{{\pi }^2}{2N}\sum _{i=1}^N\sum _{j=1}^{2N}{w}_{i}sin\left(s_i\right)F(s_i,t_j),$$

with $F(s_i,t_j)=G(x,y(s_i,t_j))J(s_i,t_j)\rho \left(y(s_i,t_j)\right)$, $t_j=-\pi +\pi (j-1)/N$, $j=1,\cdots ,2N$, $s_i=\pi (z_i+1)/2$, $i=1,\cdots ,N$ with $z_i\in (-1,1)$ the N-point Gauss–Legendre quadrature rule abscissas with corresponding weights $w_i$ for $i=1,\cdots ,N$. One proceeds similarly for (17). We consider an exact solution of Problem (2):

$$u_{\mathrm{exact}}\left(x\right)=\frac{1}{|x-x_0|},x_0\in D,$$

which consists of choosing $g\left(x^{*}\right)={\partial }_{n_{x^{*}}}{u}_{\mathrm{exact}}\left(x^{*}\right)$, for any $x^{*}\in \partial D$. The efficacy of the three-step method for various geometries (including effects of curvature) is presented in [26]. Naively implementing this method has the same computational cost as the PGQ method. We do not focus in this paper on fast implementations but do believe that it is possible to speed up this method using ideas that have been previously developed including the fast multipole method [20]. Then for simplicity, results will be computed on a sphere where the resolution of $\rho$ does not require a lot of quadrature points. One can apply the technique for arbitrary closed smooth surfaces, but might be limited by the resolution of (8). All simulations are done outside of a sphere of radius 2 using the three-step method with $N=16$ for the following representations:

• V0: standard representation (7);
• V1: modified representation (17) with the linear function ${\mathtt{v}}_1\left(y\right)=n_{x^{*}}·y$;
• V2: modified representation (17) with the Green’s function ${\mathtt{v}}_2\left(y\right)=4\pi G(y,x^{*}+n^{*});$
• V3: modified representation (17) with the quadratic function ${\displaystyle {\mathtt{v}}_3\left(y\right)=\frac{1}{2}\frac{y_1^2-y_2^2}{n_{x^{*},1}{x}_1^{*}-n_{x^{*},2}{x}_2^{*}};}$
• V4: modified representation (17) with the quadratic product function

  $${\displaystyle {\mathtt{v}}_4\left(y\right)=\frac{(y_1-5)(y_2-5)}{n_{x^{*},1}(x_2^{*}-5)+n_{x^{*},2}(x_1^{*}-5)}.}$$

Note that there are other quadratic polynomials $\mathtt{v}$ (as a function of 2 variables instead of 3, see [22] where those polynomials serve as basis for interpolation method). We make here the choice to test using similar functions as in Section 4.1.1. We solve (8) using a Galerkin method and the product Gaussian quadrature rule [36,42,43,44,45] (see Appendix B for details).The accuracy in approximating V0–V4 is limited by the accuracy of the resolution for $\rho$. This can be assessed by looking at the coefficients’ decay of the density spherical harmonic expansion. In this case the coefficients’ decay has reached ${10}^{-15}$. The results in Figure 4 show that given $\rho$ resolved, the approximation of the modified representations provide better results overall, except for V2 where the error plateaus around ${10}^{-7}$ for small ℓ (providing less accurate results compared to standard representation V0). Note that the single-layer potential commonly suffers less from the close evaluation than the double-layer potential, and the chosen method provides already a good approximation. This is the reason why the error when considering V0 decays as ℓ decreases [26]. The modified representations allow to make it even better. To better assess the efficacy of the modified representations in three dimensions, Figure 5 represents log plots of the maximum error with respect to $N\in \{8,16,24,32\}$ (the method uses $2N\times N$ quadrature points) and for various distances ℓ from point B. The results show that as $\ell \to 0$, V1–V4 allow to gain a couple of orders of magnitude in the error, even for a small N. Note that the error produced by three-step method does not seem to depend on N, and in this case there are more variations with respect to the choice of auxiliary function $\mathtt{v}$ than in two dimensions. Here, V1 (the linear function) is the best representation, producing the smallest errors (and the fastest to compute as indicated in

Table 2. Laplace 3D single-layer. CPU times (in seconds) for various number of quadrature points and representations for computing the solution (as described in Figure 4) from points A and B, for $\ell ={10}^{-k}$, $k=〚0,11〛$.

Method	V0	V1	V2	V3	V4
N = 8	0.028	0.029	0.032	0.031	0.046
N = 16	0.143	0.146	0.148	0.150	0.142
N = 24	0.352	0.344	0.346	0.35	0.356

). Again, the three-step method has been designed to treat nearly-singular integrals. It is the reason why the method provides already satisfactory results (given the resolution of $\rho$). The modified representations allow to significantly gain even more accuracy in this case, even with limited computational resources.

4.2. Scattering Problem

Using Proposition 3, we compare (19) with the modified representation (23) obtained with $\mathtt{v}\left(y\right)=e^{ik{n}_{x^{*}}·(y-x^{*})}$:

$$\begin{array}{cc}\hfill u\left(x\right)=& {\displaystyle {\int }_{\partial D}\left[{\partial }_{n_y}{G}^H(x,y)-ik(n_y·n_{x^{*}})e^{ik(n_{x^{*}}·(y-x^{*}))}{G}^H(x,y)\right]\left[\mu \left(y\right)-\mu \left(x^{*}\right)\right]d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +ik{\int }_{\partial D}[(n_y·n_{x^{*}})e^{ik(n_{x^{*}}·(y-x^{*}))}-1]G^H(x,y)\mu \left(y\right)d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +\mu \left(x^{*}\right){\int }_{\partial D}{\partial }_{n_y}{G}^H(x,y)[1-e^{ik(n_{x^{*}}·(y-x^{*}))}]d{\sigma }_y,x\in {\mathbb{R}}^d\backslash \overline{D}.}\hfill \end{array}$$

(25)

4.2.1. Example 3: Scattering in Two Dimensions

We consider an exact solution of Problem (18):

$$u_{\mathrm{exact}}\left(x\right)=\frac{i}{4}{H}_0^{\left(1\right)}\left(k\right|x-x_0\left|\right),x_0\in D,$$

which consists of choosing $f\left(x^{*}\right)=u_{\mathrm{exact}}\left(x^{*}\right)$, for any $x^{*}\in \partial D$. All simulations are done outside of a star-shaped domain using the Periodic Trapezoid Rule with $N=256$ quadrature points and $k=15$ for the following representations:

• V0: standard representation (19);
• V1: modified representation (25) (i.e., (23) with the plane wave function ${\mathtt{v}}_1\left(y\right)=e^{ik{n}_{x^{*}}·(y-x^{*})}$).

We solved (21) using Kress product quadrature rule [40] (see Appendix A). The quadrature rule is well adapted to approximate kernels with a logarithmic singularity. The accuracy of both methods is limited by the resolution for $\mu$ (the Fourier coefficients’ decay of the density is bounded by ${10}^{-6}$ for $N=256$ and $k=15$). The results in Figure 6 and Figure 7 show that given $\mu$ resolved, the approximation of the modified representation provides better results overall. Similarly to Laplace’s examples, both methods approximate well the solution far from the boundary. As the evaluation point gets closer to the boundary ($\ell \to 0$), V0 approximated with PTR suffers from the close evaluation problem leading to large errors (see [9]). Using the modified representation V1 allows to reduce the error by a couple of order of magnitude for the close evaluation problem.

Figure 8 represents log plots of the maximum error with respect to the number of quadrature points $N\in 〚50,3000〛$ and for various distances ℓ from point A (indicated in Figure 6). The results show that for any number of quadrature points, the error when considering V0 explodes as we approach the boundary (error larger than ${10}^5$) while the error with V1 remains bounded (of the order of ${10}^{-2}$ in the case presented above). In this case standard rerpresentation V0 strongly suffers from the close evaluation problem, however the modified representation V1 significantly reduces the error. Even when standard quadrature rules fail to compute the standard representation, the proposed modified one regularizes the solution and provides satisfactory results without significant additional computational time (as shown in

Table 3. Helmholtz 2D. CPU times (in seconds) for various number of quadrature points and representations. Times account for computing the solution for $N\times 12$ grid points (for $\ell ={10}^{-k}$, $k=〚0,11〛$) on a body-fitted grid.

Method	$\mathit{N}=128$	$\mathit{N}=256$	$\mathit{N}=512$
V0	0.18	0.27	0.71
V1	0.21	0.33	0.89

).

4.2.2. Example 4: Scattering in Three Dimensions

We consider an exact solution of (10):

$$u_{\mathrm{exact}}\left(x\right)=\frac{1}{4\pi }\frac{e^{ik|x-x_0|}}{|x-x_0|},x_0\in D,$$

which consists of choosing $f\left(x^{*}\right)=u_{\mathrm{exact}}\left(x^{*}\right)$, for any $x^{*}\in \partial D$.

All simulations are done outside of an ellipsoid parameterized by$y(s,t)=(2cos(t)sin(s),sin(t)sin(s),2cos(s\left)\right)$, $(s,t)\in [0,\pi ]\times [-\pi ,\pi ]$, and using the three-step method with various N. This is in order to investigate the technique in the context of limited resolution, namely the coefficients’ decay of the density spherical harmonic expansion does not reach machine precision. We consider $k=5$ and the following representations:

• V0: standard representation (19);
• V1: modified representation (25).

We solved (21) using Galerkin method and the product Gaussian quadrature rule (see Appendix B for details). The accuracy of both methods is limited by the accuracy of the resolution for $\mu$. This limitation can be checked for instance by looking at the density spherical harmonics coefficients’ decay: for $k=5$, the resolution will be capped around ${10}^{-2}$ for $N=16$, ${10}^{-4}$ for $N=24$, and ${10}^{-7}$ for $N=32$. The results in Figure 9 show that given $\mu$ resolved, standard representation incurs bigger errors at close evaluation points while the modified representation provides better results overall. Here, the resolution of the boundary integral equation was fairly limited. Figure 10 represents log plots of the maximum error with respect to $N\in 〚8,32〛$ (the method uses $2N\times N$ quadrature points) and for various distances ℓ from the boundary from point A. While the three-step method has been designed to treat nearly-singular integrals and provided satisfactory results for Laplace’s problems, the method here requires more quadrature points to achieve accuracy due to the wavenumber (see Section 4.2.3 for more details). The standard representation V0 suffers from both the close evaluation problem and the poor density resolution. The modified representation V1 allows to gain accuracy even with limited resolution (without significant additional computational time as indicated in

Table 4. Helmholtz 3D. CPU times (in seconds) for various number of quadrature points and representations. Times account for computing the solution from points A and B, for $\ell ={10}^{-k}$, $k=〚0,11〛$.

Method	$\mathit{N}=8$	$\mathit{N}=16$	$\mathit{N}=20$
V0	0.027	0.15	0.313
V1	0.03	0.15	0.314

).

4.2.3. High Frequency Behavior

It is well-known that for a fixed number of quadrature points N, accuracy is lost for larger wavenumbers k. Figure 11 and Figure 12 represent the high frequency behavior for the Examples 3 and 4, for various k and N. We consider the same quadrature rules, exact solution $u_{\mathrm{exact}}$, boundary shapes, as in Section 4.2.1 and Section 4.2.2, but we vary k and/or N. The modified representation annihilates some oscillatory behavior by subtracting plane waves along the normal of the evaluation points. It then allows a better approximation for a wider range of wavenumbers (until the number of quadrature points is not enough), and results in a greater wavenumber stability. The results in Figure 11 and Figure 12 confirm this phenomenon.

5. Modified Boundary Integral Equations

We have used (1) to modify the representation of solution of boundary value problems close to (but not on) the boundary. One could also use (1) to avoid weakly singular integrals in the boundary integral equation as done in BRIEF [31]. In the section we present a modified representation of (21).

Proposition 4.

Given $x^{*}\in \partial D$, let $\mathtt{v}$ be a solution of Helmholtz equation in $D\subset {\mathbb{R}}^d$, $d=2,3$, satisfying conditions (22). Then the boundary integral Equation (21) admits the modified representation:

$$\begin{array}{cc}\hfill & {\displaystyle {\int }_{\partial D}\left[{\partial }_{n_y}{G}^H(x^{*},y)-{\partial }_{n_y}\mathtt{v}\left(y\right)G^H(x^{*},y)\right]\left[\mu \left(y\right)-\mu \left(x^{*}\right)\right]d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +{\int }_{\partial D}{G}^H(x^{*},y)\left[{\partial }_{n_y}\mathtt{v}\left(y\right)-ik\right]\mu \left(y\right)d{\sigma }_y+\mu \left(x^{*}\right){\int }_{\partial D}{\partial }_{n_y}{G}^H(x^{*},y)\left[1-\mathtt{v}\left(y\right)\right]d{\sigma }_y}\hfill \\ \hfill & =f\left(x^{*}\right),\forall x^{*}\in \partial D.\hfill \end{array}$$

(26)

The modified representation (26) has smoother integrands than (21).

The proof can be found in Appendix C.3. Using again $\mathtt{v}\left(y\right)=e^{ik{n}_{x^{*}}·(y-x^{*})}$, Proposition 4 gives us the modified boundary integral equation:

$$\begin{array}{cc}\hfill & {\displaystyle {\int }_{\partial D}\left[{\partial }_{n_y}{G}^H(x^{*},y)-ik(n_y·n_{x^{*}})\right.\left.e^{ik{n}_{x^{*}}·(y-x^{*})}{G}^H(x^{*},y)\right]\left[\mu \left(y\right)-\mu \left(x^{*}\right)\right]d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +ik{\int }_{\partial D}{G}^H(x^{*},y)\left[(n_y·n_{x^{*}})e^{ik{n}_{x^{*}}·(y-x^{*})}-1\right]\mu \left(y\right)d{\sigma }_y}\hfill \\ \hfill & {\displaystyle +\mu \left(x^{*}\right){\int }_{\partial D}{\partial }_{n_y}{G}^H(x^{*},y)\left[1-e^{ik{n}_{x^{*}}·(y-x^{*})}\right]d{\sigma }_y=f\left(x^{*}\right),x^{*}\in \partial D.}\hfill \end{array}$$

(27)

Equation (27) has no singular integrals (in the sense its integrands have vanishing singularities), in particular it could be approximated using standard quadrature rules such as PTR in two dimensions. Going back to Examples 3 and 4 presented in Section 4.2.1 and Section 4.2.2, we now compare the approximation of the representations (19)–(25) where the density $\mu$ has been computed via (21)–(27). We then have four representations:

• V0: standard representation (19) with previous approximation of (21);
• V1: modified representation (25) with previous approximation of (21);
• V2: standard representation (19), approximation of (27) using PTR as Nyström method (2D), using product Gaussian quadrature rule (3D);
• V3: modified representation (25), approximation of (27) using PTR as Nyström method (2D), using product Gaussian quadrature rule (3D).

Figure 13 represents the results in two dimensions and illustrates how the resolution of $\mu$ limits the approximation of the solution of (18). Far from the boundary the error made using V2–V3 cannot be better than order ${10}^{-6}$. This limitation is due to the poor resolution of $\mu$ using Nyström method based on PTR to approximate (25). This can be assessed by looking at the density Fourier coefficients’ decay, which caps at ${10}^{-6}$ for $N=256$. However, as the evaluation point gets closer to the boundary ($\ell \to 0$), V3 yields competitive (sometimes better) results. Additionally, the use of Nyström PTR allows to reduce CPU times as indicated in

Table 5. Helmholtz 2D. CPU times (in seconds) for various number of quadrature points to compute the solution of the boundary integral equation.

Method	$\mathit{N}=128$	$\mathit{N}=256$	$\mathit{N}=512$
(21) with Kress product rule	0.12	0.45	1.70
(27) with PTR	0.09	0.302	1.16

. The modified boundary integral Equation (27) can be approximated using standard quadrature rules such as Periodic Trapezoid Rule (note that Nyström PTR was not possible to use to solve for (21) due to singular integrals). Its resolution may be limited but it offers interesting corrections for the close evaluation problem using simple quadrature rules as well as faster solvers.

The results in Figure 14 show that the resolution of the solution using both methods yields the same accuracy in three dimensions. The product Gaussian quadrature rule is an open quadrature at the singular point $y=x^{*}$ (see Appendix B). Thus, the modification introduced in (25) does not affect the approximation. The product Gaussian quadrature rule is a well-used, efficient, easy to implement method, but one could consider a closed quadrature rule to study the effect of (27) more closely.

6. Conclusions

In this paper, we have provided modified representations for Laplace and Helmholtz layer potentials to address the close evaluation problem in several boundary value problems. Similar to Gauss’ law, we take advantage of one auxiliary function, satisfying the partial differential equation at stake. A similar technique has been used in the context of BRIEF and density interpolation. Our approach provides guidelines on how to develop them independently of the density, and valuable insights into the layer potentials inherent nearly singular behavior. Several examples in two and three dimensions have been presented and demonstrated the efficiency of the modified representations. Given a quadrature rule, the modified representation of the solution provides a better approximation by several orders of magnitude even with limited computational resources. This assumes that the density, solution of the boundary integral equation, is sufficiently well-resolved. The modified boundary integral equation has no singular behaviors anymore, and allows us to use standard quadrature rules that do not treat singularities.

We have provided general modified representations, one can use them with any solution of their choice as long as they follow the provided guidelines to address the close evaluation. One can use this technique to modify any other wave problems, including sound-hard, penetrable obstacles. Future work includes applying those techniques to plasmonic scattering problems [46,47], deriving an asymptotic analysis to quantify the limit behavior of the error as the evaluation point approaches the boundary, as well as extensions to other partial differential equations such as Stokes problems and others.