Adaptive quadrature by expansion for layer potential evaluation in two dimensions

TL;DR

This paper introduces AQBX, an adaptive extension of QBX, for accurate layer potential evaluation in 2D boundary integral methods, with automated parameter selection and error estimation, improving practicality and precision.

Contribution

The paper develops AQBX, an adaptive scheme combining QBX with automated parameter selection and error estimation for 2D layer potential evaluation.

Findings

Achieves target error tolerance within an order of magnitude.

Simplifies parameter selection compared to original QBX.

Demonstrates effectiveness on Laplace and Helmholtz problems.

Abstract

When solving partial differential equations using boundary integral equation methods, accurate evaluation of singular and nearly singular integrals in layer potentials is crucial. A recent scheme for this is quadrature by expansion (QBX), which solves the problem by locally approximating the potential using a local expansion centered at some distance from the source boundary. In this paper we introduce an extension of the QBX scheme in 2D denoted AQBX - adaptive quadrature by expansion - which combines QBX with an algorithm for automated selection of parameters, based on a target error tolerance. A key component in this algorithm is the ability to accurately estimate the numerical errors in the coefficients of the expansion. Combining previous results for flat panels with a procedure for taking the panel shape into account, we derive such error estimates for arbitrarily shaped boundaries in 2D that are discretized using panel-based Gauss-Legendre quadrature. Applying our scheme to numerical solutions of Dirichlet problems for the Laplace and Helmholtz equations, and also for solving these equations, we find that the scheme is able to satisfy a given target tolerance to within an order of magnitude, making it useful for practical applications. This represents a significant simplification over the original QBX algorithm, in which choosing a good set of parameters can be hard.

Figures

Tables

Table 1. Table 1: Results for varying tolerance and r / h = 1 / 4 𝑟 ℎ 1 4 r/h=1/4 , comparing AQBX and direct QBX on our reference problem. Error is measured on ∂ Ω Ω {\partial\Omega} . Reported parameters are for AQBX an average over all expansion centers (avg), and for direct QBX the optimal fixed values used at all centers (opt). Note that the smallest error that we can obtain is around 10 − 12 superscript 10 12 10^{-12} . This limitation holds also when using direct QBX.

Tol.	Eval. error	Exp. order $p$	Upsamp. rate $\kappa$	Work $W$	Speedup (avg)
$ϵ$	${\parallel u_p-u\parallel }_{\infty }$	avg (opt)	avg (opt)	avg (opt)	$W_{\text{QBX}}/W_{\text{AQBX}}$
${10}^{-4}$	$1.4\cdot {10}^{-4}$	5.6 (  4)	1.1  (  2)	6.0 (  8)	1.3
${10}^{-6}$	$1.7\cdot {10}^{-6}$	7.0 (  5)	1.5  (  2)	10.4 ( 10)	1.0
${10}^{-8}$	$1.5\cdot {10}^{-8}$	8.8 (  7)	1.9  (  3)	17.0 ( 21)	1.2
${10}^{-10}$	$2.2\cdot {10}^{-10}$	10.3 (  9)	2.3  (  3)	23.2 ( 27)	1.2
${10}^{-12}$	$2.0\cdot {10}^{-12}$	12.2 ( 11)	2.6  (  4)	32.2 ( 44)	1.4
${10}^{-13}$	$1.1\cdot {10}^{-12}$	13.3 ( 11)	2.8  (  4)	37.6 ( 44)	1.2

Table 2. Table 2: Results for varying r / h 𝑟 ℎ r/h and tolerance ϵ = 10 − 10 italic-ϵ superscript 10 10 \epsilon=10^{-10} , computed in the same way as in Table   1 .

Dist.	Eval. error	Exp. order $p$	Upsamp. rate $\kappa$	Work $W$	Speedup (avg)
$r/h$	${\parallel u_p-u\parallel }_{\infty }$	avg (opt)	avg (opt)	avg (opt)	$W_{\text{QBX}}/W_{\text{AQBX}}$
$0.10$	$1.7\cdot {10}^{-10}$	8.0 (  6)	4.1  (  7)	32.8 ( 42)	1.3
$0.25$	$2.2\cdot {10}^{-10}$	10.3 (  9)	2.3  (  3)	23.2 ( 27)	1.2
$0.50$	$6.2\cdot {10}^{-9}$	13.8 ( 12)	1.7  (  2)	23.0 ( 24)	1.0
$0.75$	$5.4\cdot {10}^{-10}$	17.7 ( 20)	1.5  (  2)	27.1 ( 40)	1.5
$1.00$	$1.1\cdot {10}^{-9}$	22.0 ( 30)	1.8  (  2)	40.3 ( 60)	1.5

Table 3. Table 3: Results when solving the problem of Fig.   6 using AQBX, with 200 panels. The iteration count is the number of GMRES iterations needed to converge, and the reference value is the count when using the method of [ 14 ] . The error in the solution is measured on a circle of radius 2, and the upsampling error estimate is reported as the maximum over all panels.

GMRES tol.	AQBX tol.	Iter. count	Rel. error	Upsamp. err. est.
${ϵ}_{\text{GMRES}}$	${ϵ}_{\text{AQBX}}$	AQBX (ref)	$\frac{{\parallel u_p-u\parallel }_{\infty }}{{\parallel u\parallel }_{\infty }}$	$\max \frac{\left|{\hat{\sigma }}_{15}\right|}{{\parallel \sigma \parallel }_{\infty }}$
${10}^{-2}$	${10}^{-4}$	5 (  5)	$9.0\cdot {10}^{-03}$	$7.1\cdot {10}^{-03}$
${10}^{-4}$	${10}^{-6}$	11 ( 11)	$7.1\cdot {10}^{-05}$	$2.7\cdot {10}^{-05}$
${10}^{-6}$	${10}^{-8}$	17 ( 17)	$5.6\cdot {10}^{-07}$	$3.5\cdot {10}^{-07}$
${10}^{-8}$	${10}^{-10}$	22 ( 22)	$9.0\cdot {10}^{-09}$	$4.6\cdot {10}^{-09}$
${10}^{-10}$	${10}^{-12}$	28 ( 28)	$6.4\cdot {10}^{-11}$	$7.8\cdot {10}^{-11}$
${10}^{-12}$	${10}^{-14}$	34 ( 33)	$3.9\cdot {10}^{-13}$	$7.9\cdot {10}^{-13}$

Equations

$$u(z)=D\sigma (z)+\alpha S\sigma (z)={\int }_{\partial \Omega }\frac{\partial G(z,w)}{\partial n_w}\sigma (w)d{s}_w+\alpha {\int }_{\partial \Omega }G(z,w)\sigma (w)d{s}_w,z\in \Omega ,$$

$$\left(\frac{1}{2}I+D+\alpha S\right)\sigma (z)=f(z),z\in \partial \Omega .$$

$$u(z)={\int }_{\partial \Omega }G(z,w)\sigma (w)d{s}_w.$$

$$G(z,w)=\sum _{m=0}^{\infty }{A}_m^r(w,z_0)\cdot B_m^r(z,z_0),$$

$$|z-z_0|<r\le |w-z_0|.$$

$$\begin{array}{rl}\left|B_m^r(z,z_0)\right|& \le 1,\\ \underset{m,z}{\max }\left|B_m^r(z,z_0)\right|& =1.\end{array}$$

$$r=\underset{w\in \partial \Omega }{\min }|w-z_0|,$$

$$u(z)=\sum _{m=0}^{\infty }{a}_m\cdot B_m^r(z,z_0),|z-z_0|\le r,$$

$$a_m={\int }_{\partial \Omega }{A}_m^r(w,z_0)\sigma (w)d{s}_w.$$

$$u_p(z)=\sum _{m=0}^p{\tilde{a}}_m\cdot B_m^r(z,z_0).$$

$$u(z)-u_p(z)=\underset{e_T}{\underbrace{u(z)-\sum _{m=0}^p{a}_m\cdot B_m^r(z,z_0)}}+\underset{e_C}{\underbrace{\sum _{m=0}^p(a_m-{\tilde{a}}_m)\cdot B_m^r(z,z_0)}}.$$

$$\left|e_T\right|=\left|\sum _{m=p+1}^{\infty }{a}_m\cdot B_m^r(z,z_0)\right|\le \sum _{m=p+1}^{\infty }\left|a_m\right|.$$

$$e_T\le M(p,r)r^{p+1}{\Vert u\Vert }_{{\mathcal{C}}^{p+1}(|z-z_0|\le r)},$$

$$\left|e_T\right|\approx \left|a_{p+1}\cdot B_{p+1}^r(z,z_0)\right|\le \left|a_{p+1}\right|.$$

$$\left|e_T\right|\approx \left|a_p\cdot B_p^r(z,z_0)\right|\le \left|a_p\right|.$$

$$E_C(n,m)\approx \left|a_m-{\tilde{a}}_m\right|\ge \left|(a_m-{\tilde{a}}_m)\cdot B_m^r(z,z_0)\right|.$$

$$\Delta u$$

$$u$$

$$\varphi (z,w)=\frac{1}{2\pi }\log \left|z-w\right|.$$

$$u(z)=D\sigma (z)={\int }_{\partial \Omega }\frac{\partial \varphi (z,w)}{\partial n_w}\sigma (w)d{s}_w.$$

$$u(z)$$

$$v(z)$$

$$\left(\frac{1}{2}I+D\right)\sigma =f.$$

$$\frac{-1}{z-w}=\sum _{m=0}^{\infty }\frac{(z-z_0{)}^m}{(w-z_0{)}^{m+1}},$$

$$A_m^r(w,z_0)$$

$$B_m^r(z,z_0)$$

$$\Delta u+k^{2}u$$

$$u$$

$$\underset{r\to \infty }{\lim }{r}^{1/2}\left(\frac{\partial u}{\partial r}-iku\right)=0,$$

$${\varphi }_k(z,w)=\frac{i}{4}{H}_0^{(1)}(k|z-w|).$$

$$u(z)={\int }_{\partial \Omega }\left(\frac{\partial {\varphi }_k(z,w)}{\partial n_w}-\frac{ik}{2}{\varphi }_k(z,w)\right)\sigma (w)d{s}_w=\left(D_k\sigma -\frac{ik}{2}{S}_k\sigma \right)(z).$$

$$\frac{\partial {\varphi }_k(z,w)}{\partial n_w}$$

Full text

{AMS}

65R20, 65D30, 65D32, 65G99

1 Introduction

Integral equation methods are a class of numerical methods that are based on the reformulation of an elliptic partial differential equation (PDE) as a boundary integral equation. When applicable, this solution approach has several attractive features. Among these are: high-order discretization methods, well-conditioned linear systems after discretization, reduced number of unknowns compared to volume methods, and straightforward handling of moving boundaries.

A suitable starting point for our discussion on integral equation methods is the representation of the solution $u$ using layer potentials, which are evaluated by integrating the PDE’s fundamental solution $G$ and a layer density $\sigma$ over the domain boundary ${\partial\Omega}$. We represent $u$ as a linear combination of the double layer potential $D$ and the single layer potential $S$,

[TABLE]

where ${n}_{w}$ denotes the unit normal pointing into $\Omega$. For a Dirichlet problem with boundary condition $u=f$ we get, considering the limit $\Omega\ni z\to{\partial\Omega}$ of (1) leads to a second kind integral equation in $\sigma$,

[TABLE]

This is a generic form for a Dirichlet problem, and we will later define the forms for the interior Laplace and the exterior Helmholtz equations. The constant $\alpha$ is selected such that the integral equation has a unique solution and is well-conditioned [7]. Nyström discretization of this equation using a suitable quadrature rule generates a dense linear system that can be solved rapidly by exploiting the fact that the off-diagonal blocks are of low rank. Solution methods include the fast multipole method (FMM) [12] and fast direct methods [11, 17, 19].

The main difficulty in the procedure outlined above lies in finding a quadrature rule that can evaluate the layer potentials $S\sigma(z)$ and $D\sigma(z)$ when the target point $z$ is close to or on the boundary ${\partial\Omega}$. In this case, the integrals of the layer potentials are nearly singular or singular, requiring specialized quadrature methods. For two-dimensional (2D) problems, several efficient such methods are available; see the summaries [13] and [14]. Some of the 2D methods can be extended to axisymmetric surfaces in three dimensions (3D) [15]. However, for general surfaces in 3D the available methods are not as mature as those in 2D. We refer to [23] for a recent summary of the current state of the art.

1.1 QBX

In this paper, we will focus on the relatively recent method of quadrature by expansion (QBX) [5, 18]. The method provides a way of evaluating both singular and nearly singular integrals by representing the layer potential as a local expansion, centered at a point some distance away from the boundary. While originally proposed for Helmholtz in 2D, it can be generalized to other PDEs in 2D and 3D; see [1, 3, 24] for its successful use in several applications. A strength of QBX is that it uses the same type of local expansions as the FMM, which allows the two methods to be combined into a fast method for evaluating layer potentials at arbitrary locations. This is a topic of ongoing research [22]. As with most methods, however, QBX solves one problem and introduces another. While the problem of evaluating layer potentials on or very close to the boundary is solved, it is replaced with the new problem is how to efficiently compute the local expansion of the layer potential. In particular, computing the expansion coefficients entails evaluating a series of integrals with increasing order of near singularity. Effectively one has traded a hard problem for several easier problems. Nevertheless, QBX is still a competitive method for these problems, especially since it has a solid analytical groundwork [9].

One of the difficulties of QBX is that of parameter selection. The convergence of the local expansion is governed by the order $p$ at which it is truncated, and by the distance $r$ between the boundary and the expansion center. The expansion coefficients are computed using a quadrature rule for smooth integrands, and for them to be accurate it is necessary to upsample the boundary points by some factor $\kappa$. These three parameters together affect the two competing errors of QBX: the truncation error and the coefficient error (often referred to as the quadrature error). The truncation error increases with $r$ and decreases with $p$, while the coefficient error increases with $p$ and decreases with $\kappa$ and $r$ (see [1, Fig. 3] for an example). Together $r$, $p$ and $\kappa$ form a large parameter space, and how to best set these parameters is not clear. Instead, experimentation must be used to determine good parameter ranges for a specific application. This is in itself not unfeasible, but it would be preferable to reduce the number of free parameters.

1.2 Contribution

In this paper we propose a scheme for adaptively setting the order $p$ and upsampling factor $\kappa$ at the time of computation, such that a the error is maintained below a target tolerance. The key ingredient for this scheme to be successful is the ability to accurately estimate the magnitude of the coefficient error, which is the quadrature error in the expansion coefficients. Such estimates were derived in [2] for simple geometries in two dimensions, namely flat Gauss-Legendre panels and the trapezoidal rule on the unit circle. Here we extend these estimates, by locally using a polynomial to represent the mapping between a flat panel and a panel of general shape. This greatly increases the accuracy of the estimates, and allows us to build an adaptive QBX scheme based on them. We here restrict ourselves to analyzing Gauss-Legendre panel quadrature, but the methodology could equally well be applied to a discretization using the trapezoidal rule.

Taking into account the mapping of the parametrization when analyzing nearly singular quadrature errors is by itself not new; a discussion on nearly singular quadrature errors similar to ours can be found in [5]. Our main contribution in this regard is the construction of an explicit representation of the mapping, which we then invert in order to compute an error estimate.

This paper is organized as follows: In Section 2 we introduce the general structure of our scheme, and derive the details for the Laplace double layer potential and the Helmholtz combined field potential. In Section 3 we show how to compute the coefficient error estimates necessary for the scheme to be useful. In Section 4 we briefly discuss how to combine the scheme with a fast method. In Section 5 we show a selection of numerical results that illustrate the performance of our method.

2 Foundations of AQBX

In this section we begin by introducing the foundations of our method using a generic notation, before we give the specific details for the Laplace and Helmholtz equations. We start from a given layer potential representation,

[TABLE]

To evaluate this using AQBX, we first need to split the fundamental solution using a suitable addition theorem,

[TABLE]

where $A_{m}^{r}$ and $B_{m}^{r}$ are either vectors or scalars, depending on the fundamental solution. For a given $r$, the above addition theorem is valid for

[TABLE]

For the specific formulas for Laplace and Helmholtz, see (24) and (35). We assume that (4) is normalized such that the following holds for $|z-z_{0}|\leq r$ and integer $m\geq 0$:

[TABLE]

Here, and throughout the paper, we let $|\cdot|$ denote the $\ell^{2}$-norm for vectors, and the complex modulus for scalars. If the second condition holds for each $m$, i.e. $\max_{z}\mathinner{\!\left\lvert B_{m}^{r}(z,z_{0})\right\rvert}=1$, then we get slightly sharper error estimates in Algorithm 1. This is however not necessary, and indeed only holds for the Laplace formulation used in this paper. The addition theorem allows us to pick an expansion center $z_{0}$, determine $r$ as

[TABLE]

and evaluate the layer potential as a local expansion in the neighborhood of $z_{0}$,

[TABLE]

where

[TABLE]

The fact that (8) holds also at the equality $|z-z_{0}|=r$ was shown in [9] for the Laplace and Helmholtz equations, but can be generalized to other kernels, e.g. Stokes equations [1]. This allows us to evaluate the expansion also on ${\partial\Omega}$, at the single point which is closest to $z_{0}$. In fact, in [5] it is shown that the expansion even converges in a disc about $z_{0}$ with radius $R>r$, as long as the density $\sigma$ is analytic inside that disc.

In a computational scheme the local expansion is truncated at a maximum order $p$, and the coefficients, which we now denote $\tilde{a}_{m}$, are computed using a suitable quadrature rule. This gives us the QBX approximation of the layer potential,

[TABLE]

The error in this approximation can, assuming exact arithmetic, be separated into a truncation error $e_{T}$ and a coefficient error $e_{C}$ [9],

[TABLE]

2.1 Truncation error

The truncation error $e_{T}$ of the scheme arises because we limit the local expansion to $p+1$ terms. In our normalized form (6) it can be bounded by

[TABLE]

It was shown in [9] that the truncation error for several kernels satisfies

[TABLE]

for some positive $M(p,r)$. In our experience, the truncation error typically decays exponentially in $p$, with a rate that is proportional to $r$. Finding a usable a priori estimate for the error is hard, since it depends on both the local geometry of ${\partial\Omega}$ and the regularity of the density $\sigma$ in a nontrivial way. However, assuming that the expansion coefficients decay exponentially, a good estimate for the truncation error is

[TABLE]

In practice we only have the coefficients up to $a_{p}$. Therefore, we can define a useful (and usually conservative) a posteriori estimate as

[TABLE]

2.2 Coefficient error

The coefficient error $e_{C}$ in (11) is a result of the numerical evaluation of the coefficient integrals (9) for $m=0,\dots,p$, which for a given boundary ${\partial\Omega}$ is evaluated using some $n$-point quadrature rule. Assuming that the density $\sigma$ and the boundary ${\partial\Omega}$ are well-resolved by that quadrature (which is a prerequisite for the underlying Nyström discretization), the main source of the error is the near singularity in $A_{m}^{r}$ when evaluated at $z_{0}$. The order of this near singularity typically grows with $m$, and to counter this the density must be upsampled (by interpolation) to a grid which is fine enough to resolve $A_{m}^{r}$. The amount of upsampling required can be determined through $E_{C}(n,m)$, which is an accurate a priori estimate of the coefficient error in term $m$ when using $n$ quadrature nodes,

[TABLE]

We will in section 3 show how to derive such error estimates for the cases of the Laplace and Helmholtz equations.

Remark 1** (Resolution error).**

*Here we only discuss errors in the coefficients $\tilde{a}_{m}$ due to near singularities in the integrals (9). However, there is also a lower bound on the accuracy of the coefficients, imposed by how well the underlying grid resolves the boundary ${\partial\Omega}$ and the layer density $\sigma$. For a given panel, this error could be estimated by analyzing a modal expansion of the grid point coordinates and density values, such as the Legendre polynomial expansion used in Section 3. This is briefly explored in section 5.3, though a more in-depth analysis is beyond the scope of the present paper. *

2.3 The AQBX scheme

Let us now assume that we have a discretization of ${\partial\Omega}$ characterized by $n$, which for a global quadrature denotes the total number of points, and for a panel-based quadrature denotes the number of points on each panel. Denoting that quadrature $\operatorname{Q}_{n}$, we define a combined interpolation and quadrature operator $\operatorname{Q}_{\kappa n}$, which computes the quadrature by first upsampling the density to $\kappa n$ points (for simplicity we assume $\kappa$ integer). Furthermore, we assume that the error when computing a coefficient $a_{m}$ is well estimated by a function $E_{C}(\kappa n,m)$. Then the adaptive algorithm for evaluating $u(z)$ in the neighborhood of $z_{0}$ to a tolerance $\epsilon$ can be summarized using Algorithms 1 and 2.

Note that if Algorithm 1 has produced $p+1$ coefficients, then Algorithm 2 is guaranteed to terminate within $p+1$ iterations, since $|\delta_{p}|\leq|a_{p}|<\epsilon$. For more conservative termination criteria, one can use $\max(|a_{m-1}|,|a_{m}|)<\epsilon$ and $\max(|\delta_{m-1}|,|\delta_{m}|)<\epsilon$ to take into account any even/odd behavior in the expansion.

2.4 Specific applications

We will now proceed with formulating AQBX for two different applications: the Laplace double layer potential and the Helmholtz combined field potential.

2.4.1 Laplace equation

We first consider the Laplace Dirichlet problem in a domain $\Omega$ bounded by a boundary ${\partial\Omega}$,

[TABLE]

The fundamental solution to this PDE is

[TABLE]

The interior Dirichlet problem can be represented using the double layer potential,

[TABLE]

In complex notation this can be compactly represented as

[TABLE]

On ${\partial\Omega}$ the density $\sigma$ satisfies the integral equation

[TABLE]

It should be noted that this particular layer potential actually has a smooth limit on the boundary, so no special quadrature is needed for solving the integral equation, unless different parts of the boundary are close to each other. Special quadrature is however still needed for evaluating the solution close to the boundary, once $\sigma$ has been computed.

Starting from the Taylor expansion of the Cauchy kernel,

[TABLE]

it is straightforward to derive the AQBX formulation (4) for $v(z)$,

[TABLE]

Note that, by these definitions together with (9), the real part of the coefficient $a_{0}$ will simply hold the value of the double layer potential evaluated at $z_{0}$, i.e. $\operatorname{Re}a_{0}=u(z_{0})$.

2.4.2 Helmholtz equation

We now consider the Helmholtz Dirichlet problem in an unbounded domain $\Omega$ exterior to a boundary ${\partial\Omega}$, which for a wavenumber $k$ is stated as

[TABLE]

The solution must satisfy the Sommerfeld radiation condition for $r=|z|$,

[TABLE]

which gives a fundamental solution that is essentially the zeroth-order Hankel function of the first kind,

[TABLE]

It is possible to represent the solution using the combined field integral representation111The general form of the representation is $D_{k}-i\eta S_{k}$. We have here set $\eta=k/2$, following [14].,

[TABLE]

Here the double layer kernel is

[TABLE]

On the boundary we get the integral equation

[TABLE]

To formulate AQBX for the combined field representation, we start from the Graf addition theorem [20, §10.23(ii)],

[TABLE]

Here $(r_{w},\theta_{w})$ and $(r_{z},\theta_{z})$ are the polar coordinates of $w-z_{0}$ and $z-z_{0}$,

[TABLE]

We can thus form an expansion for the kernel of the combined field representation (31) as

[TABLE]

Here $s_{m}$ is an immediate result of (35),

[TABLE]

while $d_{m}$ is obtained by differentiation of (39). A compact form for $d_{m}$ was derived in [5],

[TABLE]

The Bessel functions $J_{m}$ decay rapidly with $m$. It however easy to experimentally verify that a normalization satisfying (6) is obtained by using the first term of the power series for $J_{m}(kr)$ [20, §10.4,§10.8],

[TABLE]

We have defined our general AQBX formulation (4) for indices $m\geq 0$. To fit the expansion (37) into this we let $A_{m}^{r}$ and $B_{m}^{r}$ be vector-valued for $m>0$, representing the terms with indices $\pm m$ in (37). The AQBX formulation for the combined field representation is then given by

[TABLE]

and, for $m>0$,

[TABLE]

such that the coefficients $a_{m}\in\mathbb{C}^{2}$ for $m>0$.

We again note that the coefficient $a_{0}$ (9) will hold the value of the potential at $z_{0}$, since (trivially) $s_{0}=\phi_{k}$ in (30) and (through (33))

[TABLE]

3 Coefficient errors

In this section we derive a central piece of AQBX: the coefficient error estimate $E_{C}$ required in Algorithm 1. We will consider the layer potentials of section 2.4 combined with a panel-based quadrature, where the boundary curve is subdivided into smaller segments, each of which is discretized using an $n$-point Gauss-Legendre quadrature. This is sometimes referred to as a composite Gauss-Legendre quadrature.

We begin by considering the error in the contribution to a coefficient from a single panel; the total error can then be computed as the sum of errors from all adjacent panels. For this, let $\Gamma$ be an open curve (i.e. a panel) parametrized by an analytic function $\gamma(t)\in\mathbb{C}$, $t\in[-1,1]$, with the normal defined as ${n}(t)=i\gamma^{\prime}(t)/|\gamma^{\prime}(t)|$, and $\gamma$ oriented such that ${n}$ points into the domain $\Omega$. For simplicity we denote by $\sigma(t)$ the pullback of $\sigma$ under $\gamma$, i.e. $\sigma(t)=\sigma(\gamma(t))$.

3.1 Laplace coefficient error

Beginning with QBX for the Laplace double layer potential (sec. 2.4.1), the expansion coefficients are computed as

[TABLE]

The coefficient errors are introduced when this integral is evaluated using a discrete quadrature rule. In a boundary integral equation context, we assume a panel $\Gamma$ such that $\gamma$ and $\sigma$ are well-resolved by an $n$-point Gauss-Legendre quadrature rule. We will here focus on the standard choice $n=16$.

For our discussion on quadrature errors, we introduce the following notation: Let $\operatorname{I}$ denote the integral over $[-1,1]$, and let $\operatorname{Q}_{n}$ denote the Gauss-Legendre approximation of that integral,

[TABLE]

The quadrature error $\operatorname{R}_{n}$ is then defined as

[TABLE]

and $\operatorname{I}$, $\operatorname{Q}_{n}$ and $\operatorname{R}_{n}$ are all linear functionals on $C(-1,1)$.

If the above assumptions on the resolution hold, then we can expect the quadrature error to be small for the integrand (48), provided that $z_{0}$ is far away from $\Gamma$. If on the other hand $z_{0}$ is close to $\Gamma$, then the quadrature error will be dominated by the nearly singular quadrature error that arises when the integrand is evaluated close to its pole.

To estimate nearly singular quadrature errors, we can can proceed in the same way as in [2]. The central property which we will use is the following: Let $f(t)$ be a function which is analytic on $[-1,1]$ and everywhere inside a contour $\mathcal{C}$ enclosing $[-1,1]$, except at a finite number of poles $\{t_{j}\}_{j=0}^{N}$ enclosed by $\mathcal{C}$. If $f$ is integrated on $[-1,1]$ using the $n$-point Gauss-Legendre rule, then the quadrature error (51) is given by [8]

[TABLE]

where $k_{n}$ is the characteristic remainder function,

[TABLE]

While not available in closed form, it can in the limit $n\to\infty$ be shown to satisfy [6, appx.]

[TABLE]

with the sign in the denominator given by the sign of $\operatorname{Re}t$. This result, though asymptotic, is accurate already for small $n$ [2]. Since $n\geq 1$, we can bound the remainder function as

[TABLE]

for some $C(t)>0$. If $t$ lies on the Bernstein ellipse $B_{\rho}$, defined as the ellipse with foci $\pm 1$ where the semimajor and semiminor axes sum to $\rho>1$, then $|t\pm\sqrt{t^{2}-1}|=\rho$ and

[TABLE]

where $C_{\rho}=\max_{t\in B_{\rho}}C(t)$. If $|f|\leq M_{\rho}$ on some Bernstein ellipse $B_{\rho}$, and we let the contour integral in (52) follow that ellipse, then the integral can be bounded as

[TABLE]

The last above step is motivated by the geometry of a Bernstein ellipse; the circumference of $B_{\rho}$ has its maximum as $\rho\to 1$, in which case the ellipse collapses onto $[-1,1]$. The contribution from the contour integral in (52) vanishes completely if $f(t)k_{n}(t)t\to 0$ as $|t|\to\infty$, and is generally small compared to the residue for poles close to $[-1,1]$, as was noted in [2].

In the case where $f$ has a simple pole $t_{0}$, we can bound the error as

[TABLE]

Now consider the more general case, when $f$ has a single pole $t_{0}$ of order $m+1$ (the case of several poles follows trivially). For the residue we can write

[TABLE]

since $k_{n}(t_{0})$ is smooth for $t\notin[-1,1]$. Defining

[TABLE]

the residue can be bounded as

[TABLE]

This allows us to bound the quadrature error as

[TABLE]

The above bound will overestimate the error by a large factor. In practice a good approximation of the error is achieved if we neglect the contribution from the contour, under the assumption that the pole is close to the interval, and assume that $k_{n}$ varies much more rapidly with $t$ than the other factors. Only keeping the term in (59) with the highest derivative in $k_{n}$, we get the error approximation

[TABLE]

This, with the remainder function evaluated using (54), is the approximation that we will be using henceforth. The derivatives of $k_{n}$ can (at least for small $m$) be derived analytically from (54), or, as shown in [2], approximated as

[TABLE]

Returning to the QBX coefficients of the Laplace double layer potential, the integral (47) clearly has a pole of order $m+1$ at $z_{0}$. In parametrized form (48) the pole instead lies at the point $t_{0}\in\mathbb{C}$, such that

[TABLE]

Note that evaluating $\gamma(t_{0})$ corresponds to evaluating an analytic continuation of $\gamma(t)$ in some neighborhood of $[-1,1]$. Denoting by $\tilde{a}_{m}$ the approximation of $a_{m}$ (as defined in (48)), and also assuming that there exists an analytic continuation of $\sigma$, then we have from (63) that the error can be approximated as

[TABLE]

This we can simplify by noting that since $z_{0}=\gamma(t_{0})$,

[TABLE]

Thus a coefficient error estimate $E_{C}(n,m)\approx|a_{m}-\tilde{a}_{m}|$ is given by combining (54), (64), (66) and (67),

[TABLE]

This result is essentially a generalization of [2, Thm. 1] to curved panels.

The estimate (68) holds well in the limit $n\to\infty$, but will lose accuracy with increasing $m$ for a fixed $n$. We have in practice found that it starts losing its reliability around $m>n/2$, so that it then is safer to trigger an upsampling in Algorithm 1, instead of continuing to rely on the estimate.

Though the above derivation is rather straightforward, some more work is required if we want to evaluate (68) in practice, since that requires finding $t_{0}$ and being able to evaluate $\sigma(t_{0})$ and $\gamma(t_{0})$. For this, one would ideally like to have access to analytic expressions for $\sigma$ and $\gamma$, but what we typically have is instead the values of the functions at the quadrature nodes. To be able to evaluate (68) using this pointwise data, we first need to construct continuations of $\gamma$ and $\sigma$. Then we need to solve (65) using the continuation of $\gamma$, in order to find $t_{0}$ and finally evaluate the estimate. This may sound hard, but can actually be done in an efficient way using polynomial extrapolation.

Beginning with the continuation of $\gamma$, let $t_{i}$ and $w_{i}$, $i=1\dots n$, be the nodes and weights of the n-point Gauss-Legendre quadrature. We can then let $\operatorname{P}_{n}[\gamma](t)$ be the polynomial of degree $n-1$ that interpolates $\gamma(t)$ at the nodes $t_{i}$. For this high-order ($n=16$) interpolation to be accurate, we need to compute $\operatorname{P}_{n}[\gamma]$ in a way that is both well-conditioned and stable [27, Ch. 14]. The distribution of the Legendre points $t_{i}$ ensures that our interpolation problem is well-conditioned. For stability, we use as our basis the Legendre polynomials $P_{\ell}(t)$,

[TABLE]

These are orthogonal on $[-1,1]$, allowing us to explicitly compute the coefficients $\hat{\gamma}_{\ell}$ as

[TABLE]

where, for a given $n$, $L$ is a constant matrix with elements

[TABLE]

To improve convergence in the following step, we assume that $\Gamma$ has endpoints at $-1$ and $1$; this can be achieved for any open curve by first applying a simple scaling and rotation to both $\Gamma$ and $z_{0}$. Letting the interpolant $\operatorname{P}_{n}[\gamma](t)$ be the analytic continuation of $\gamma(t)$, we can now find $t_{0}$ in (65) by instead solving

[TABLE]

using a numerical root-finding algorithm. This does, for the purposes of error estimation, work very well in the near neighborhood of $\Gamma$, see example in Fig. 1. The Legendre polynomials $P_{\ell}(t)$ are evaluated using the standard recurrence relation [20, §14.10.3], which holds also for complex $t$ [20, §14.21(iii)]. The polynomials grow exponentially in $\ell$ for $t\notin[-1,1]$, but we have not observed any stability issues related to this. In particular, the coefficients $\gamma_{\ell}$ of a well-resolved panel decay faster than the polynomials $P_{\ell}(t)$ grow for the range of $t$-values that are of interest to us.

Equation (72) can be solved efficiently using Newton’s method, but the choice of starting guess can be important, especially near concave regions of the curve (e.g. below the curve of Fig. 1). In such regions the inverse mapping $t_{0}=\gamma^{-1}(z_{0})$ is no longer single-valued, and for our estimate to be accurate we want the solution $t_{0}$ that predicts the largest error, when the estimate is evaluated at that point. This “worst solution” can often be found by solving with $t_{0}=\pm 1$ as starting guesses and comparing the results; this strategy was used for the results shown to the right in Fig. 2. As an alternative, one can use $t_{0}=z_{0}$ as as starting guess (given the above assumptions on the endpoints of $\Gamma$). This simpler strategy works well for practical purposes, and is what we use throughout the remainder of this paper.

Once $t_{0}$ is found, we are able to evaluate the coefficient error estimate as given in (68). The value of $\gamma^{\prime}(t_{0})$ can be found by differentiation of the interpolant,

[TABLE]

We can also evaluate $\sigma(t_{0})$ through an interpolating polynomial,

[TABLE]

This works very well if $\sigma$ is well-resolved on $\Gamma$, and can be used to obtain the fine-scale correspondence between error and estimate seen in Fig. 3. A less expensive alternative is to use the max norm of $\sigma$ on the interval/panel,

[TABLE]

This works well in practice, as we mainly need to get the order of magnitude right, and also appears to be slightly more robust whenever $\sigma$ is not fully resolved on $\Gamma$, especially for Helmholtz (see Figs. 4 and 5). We use this approach in the results of section 5.

The above results can also be used for estimating the nearly singular quadrature error when evaluating the double layer potential (20) near $\Gamma$, and hence to determine when AQBX needs to be used. To that end, we simply use the observation that $u(z_{0})=\operatorname{Re}a_{0}$. Denoting by $\tilde{u}(z_{0})$ the potential evaluated using direct quadrature, we thus have from (66) that

[TABLE]

This estimate, when evaluated using the above procedure, is very accurate. As an example, in Fig. 3, we consider the Laplace Dirichlet problem in a starfish domain. The density is computed by solving (23) using the Nyström method and a right-hand side given by a collection of point sources (marked by + in the figure). The solution is then evaluated inside the domain directly using the composite Gauss-Legendre quadrature, and compared to the reference solution given by the boundary condition. Comparing these results to those in [2, fig. 7], which used a flat panel approximation, shows the importance of taking into the account the inverse mapping of the target point.

3.2 Helmholtz coefficient error

We now turn to AQBX for the Helmholtz combined field potential (sec. 2.4.2). When evaluating the expansion coefficients, the Hankel functions in the integrands have a singularity as $r_{w}\to 0$, which to leading order behaves like [20, §10.4,§10.8]

[TABLE]

The quadrature error due to near singularity is dominated by that from the highest-order pole, so for the purposes of error estimation it is suitable to only keep the highest-order Hankel function in the expression for (44), and to approximate that Hankel function using only the first term of (77). We can then approximate (39) as

[TABLE]

and (40) as

[TABLE]

and, for $|m|>0$ and using (36),

[TABLE]

Inserted into (42) and (44), this gives

[TABLE]

Note that we in the above series of simplifications (79–84) have removed constant factors which are irrelevant to the magnitude of the error (i.e. $-1$ and $i$). Estimating the error of these simplified integrands on a panel with parametrization $\gamma$, it follows that

[TABLE]

Note that (85) is bounded by (86) with $m=0$, so it is sufficient to use only (86), if we can accept being on the conservative side for $m=0$. Also note that the integrand in (86) is exactly the same as in the error estimate for the Laplace double layer potential (66). The conclusion is that, remarkably, the coefficient error estimate $E_{C}$ for Helmholtz is identical to that previously derived for Laplace (68),

[TABLE]

This correspondence between Laplace and Helmholtz also holds for the quadrature error when evaluating the underlying potential itself. This can be seen from (85) and the observation from (46) that $a_{0}=u(z_{0})$,

[TABLE]

This estimate has been derived by simplifying the Helmholtz double layer kernel (79), and noting that the singularity to leading order is identical to that of the Laplace double layer. However, in experiments we can observe that the small-scale oscillations predicted by (88) (and which are clearly noticeable in the Laplace case) only appear for small wavenumbers, and it is therefore generally better to use (89). The reason for the disappearance of the oscillations is unknown to us, though an interesting feature is that they only seem to disappear when $\sigma$ is the solution to (34). Setting $\sigma=1$ and computing the error by comparing to a finer grid produces the oscillations also for large wavenumbers.

As a demonstration, in Fig. 4 we repeat the experiment of Fig. 3, but for the Helmholtz exterior Dirichlet problem. We set a number of point sources (marked $*$) inside a starfish domain and solve the integral equation (34) using the discretization scheme of [14]. The correspondence between error and estimate is still very good, though not as excellent as in the Laplace case. Fig. 5 shows that the estimate appears to work also for a source point very close to the boundary, when the accuracy of the solution $\sigma$ has started to deteriorate. Here, when $\sigma$ is no longer well-resolved, one can clearly see that evaluating $\sigma(t_{0})$ using the panel max norm (75) is more stable than using polynomial extrapolation (74). Using (75) results in a useful error estimate close the boundary, where the nearly singular quadrature error dominates. Further away the error is due to $\sigma$ not being accurate, something which our estimate does not take into account.

It is worth noting that the error estimate (89) is independent of the wavenumber $k$. This might come as a surprise, as one usually needs to increase the grid resolution with increasing wavenumber. However, here we only take into account the nearly singular quadrature error, under the assumption that far-field interactions are well-resolved. The result simply reflects the fact that the singularity in the kernel is independent of $k$.

4 Local AQBX

Since QBX is a special quadrature scheme for target points that are close to or on the boundary ${\partial\Omega}$, it makes sense to only use QBX for those parts of the boundary that are close to a given target point. This is known as “local QBX” [21] (as opposed to “global QBX”), and can be particularly straightforward to combine with a fast method. For panel-based quadrature on a simple curve this is easy to implement; only the panels that are near a given expansion center are used in the local expansion. Selecting panels to include can be done using an error estimate of the layer potential, such as (76) in combination with a tolerance, or by simply including a fixed number of neighboring panels (this works well if all panels are of equal length). When evaluating the potential, the contribution from the near panels is computed through the local expansion, while the contribution from the remaining panels is computed directly using the underlying Gauss-Legendre quadrature.

Let the boundary be composed of a set of panels $\Gamma_{i}$,

[TABLE]

We can then denote by ${\mathcal{N}}$ the near panels that are included in the local expansion, and by ${\mathcal{F}}$ the far panels that are evaluated directly. Note that this division must be made such that the panels in ${\mathcal{F}}$ are well-separated from all target points at which the expansion will be evaluated. If we write the layer potential as

[TABLE]

then the numerical approximation of $u$ using local QBX can be written as

[TABLE]

To combine this with a fast method that directly evaluates the interactions between all source points (such as the FMM), one can simply subtract the direct contribution from the near panels,

[TABLE]

The last two terms in this expression can together be viewed as a correction term to the direct quadrature. Since this correction only has local support, the local QBX scheme is FMM compatible, in the sense introduced in [13]. This “black box” way of using the fast method could potentially introduce cancellation errors, when the direct contribution from the near panels is added and subtracted, though we have not observed any such problems in practice. The alternative is to modify the FMM to ignore those local contributions in the first place, which is non-trivial for complex geometries.

A subtle feature of local QBX is that the width of the segment ${\mathcal{N}}$ affects the convergence rate of the local expansion. This is due to artificially induced endpoint singularities at the interfaces between ${\mathcal{N}}$ and ${\mathcal{F}}$, and has been discussed to some extent in both [5] and [21]. The choice of width of ${\mathcal{N}}$ is therefore a balance; widening ${\mathcal{N}}$ means more points contributing to each expansion, while narrowing ${\mathcal{N}}$ gives a slower convergence of the expansion. Actually, not only is the convergence slower, the exponential convergence of the truncation error (13) also tends to be less regular. This in turn makes it harder for AQBX to correctly determine when to terminate. We have found that a good balance is struck by using the five panels that are closest to the expansion center.

5 Numerical experiments

We have implemented the above algorithms in Matlab, and used the FMM as implemented in FMMLIB2D [10] for fast far-field evaluations. Timings will not be reported here, as our code is a proof of concept rather than a production implementation. We will in the following numerical experiments only report on the Helmholtz problem, as that is the more challenging one. Carrying out the same experiments for the Laplace problem just results in similar, though slightly better, results.

In our numerical experiments we will mainly use AQBX as a post-processing tool, meaning that we use it to evaluate the layer potential given by a known density. This allows us to study the behavior of AQBX in isolation, without taking into account the method used to compute the density. However, in section 5.4 we show that AQBX can be used also for solving the integral equation.

For our experiments we set up the reference problem shown in Fig. 6: The Helmholtz problem in the domain exterior to the starfish curve $\gamma(t)=(1+0.3\cos(10\pi t))e^{-2\pi it}$, $t\in[0,1]$, with Dirichlet boundary condition given by the potential from five point sources in the interior domain. The point sources are randomly positioned on a circle with radius 0.2 centered at the origin, with strengths that are randomly drawn and then normalized such that $\mathinner{\!\left\lVert u\right\rVert}_{\infty}=1$ on ${\partial\Omega}$. This way all errors reported below are both relative and absolute.

We discretize the boundary using 200 Gauss-Legendre panels of order 16 and equal arc length $h$, and position one expansion center at a distance $r=h/4$ in the normal direction from each point on the boundary. The wavenumber is set to ${k}=2/h=44.36$. The density $\sigma$ is computed by solving the integral equation (34) using the Nyström method of [14], which gives us a solution that has a relative error of approximately $10^{-14}$ when using direct quadrature away from the boundary (measured on a circle of radius 2). Increasing the number of panels or decreasing ${k}$ does not significantly improve the accuracy of the solution, nor the accuracy of the QBX evaluation.

Once we have computed $\sigma$, we can evaluate the layer potential using AQBX, and compare the result to the exact solution given by the potential from the five point sources. The error in AQBX has in our tests always been largest when evaluating the layer potential on the boundary (where the integral is singular), so we mainly report the errors as measured there.

5.1 Performance of the algorithm

We first perform a few experiments that illustrate how AQBX works. Figure 7 shows the error along the entire boundary when evaluating the solution using AQBX and a tolerance of $10^{-10}$. One can clearly see how the magnitude of the coefficients $a_{m}$ provides a good overestimate of the truncation error, while the coefficient error is closely tracked by the estimate $E_{C}$. Figures 8 and 9 show example results from when AQBX is used for evaluating the potential in the domain, where the integral is nearly singular. It can be seen that AQBX is only activated at the points where it is needed, and that the potential is then evaluated to the desired accuracy at those points.

The lowest error that we can achieve using QBX (both adaptive and direct) for this problem is around $10^{-12}$, which can also be seen in Table 1. This presents a loss of two digits of accuracy compared to the error of $10^{-14}$ achieved using the direct quadrature away from the boundary. While we can not immediately explain this loss, it is reminiscent of the results in [18, Tab. 1]. There, the error when using QBX for evaluating the double layer is two orders of magnitudes larger than for the single layer.

5.2 Comparison to direct QBX

We believe that the main benefit of using AQBX rather than a direct QBX implementation is that the parameter choice is greatly simplified; given an expansion distance $r$ and a tolerance $\epsilon$, the upsampling rate $\kappa$ and expansion order $p$ are set on the fly as needed at each expansion center. A second benefit is that setting $\kappa$ and $p$ on the fly can save some work, compared to using fixed values everywhere. In an attempt to quantify this, we now introduce a measure of the work (${W}$) needed to form a local expansion, in terms of source evaluations per original source point. If direct QBX is used with order $p$ and a fixed upsampling rate $\kappa$, then the work is given by

[TABLE]

If AQBX is used to compute $p$ coefficients, with upsampling rate $\kappa_{m}$ used to evaluate the $m$th coefficient, then the work is given by

[TABLE]

As a comparison, in Tables 1 and 2 we measure the work when computing all expansion coefficients in our reference problem using AQBX, and compare that to the work needed if $p$ and $\kappa$ are fixed everywhere to the minimum values required to achieve the same accuracy as AQBX. These fixed values are tuned by hand to the optimal values for this specific problems, but our algorithm still gives a slight speedup in our definition of work. More importantly, our algorithm is in most cases able to keep the error at the desired order of magnitude without any manual intervention. Hand-tuning parameters is on the other hand strictly unfeasible in real applications, and will generally result in an unnecessarily conservative choice of parameters.

Our measure of work does not take into account the extra effort needed to evaluate the estimate of the AQBX scheme. The reported speedup should therefore be viewed as an upper limit compared to the optimal parameter set, and as an indicator that automatic parameter selection does not necessarily have to be more expensive than using fixed parameters. To minimize the overhead of the scheme one can evaluate the error estimates Eqs. 68 and 87 recursively, and the multiple levels of upsampling can for each panel be computed and reused as needed using a caching algorithm.

5.3 Source points close to the boundary

In order to push the limits of the algorithm, we run a sequence of successively harder test cases. In each case, the Dirichlet boundary condition is given by the potential from a source point that lies at a distance $d$ from the boundary, with $d$ becoming successively smaller. Geometry and results are shown in Fig. 10. Setting the AQBX tolerance to $10^{-12}$, which is the lowest that we can achieve, we see that the algorithm starts losing accuracy once the source point is within a distance $d=h$ from the boundary. Meanwhile, the error when using direct quadrature, measured away from the boundary, does not start growing until $d=h/2$.

The computation of the QBX coefficients relies on upsampling the density $\sigma$ on each panel using polynomial interpolation. If the density is not well represented by its interpolant, then that will limit the accuracy of the QBX coefficients. One way of approximating the accuracy of this interpolation is by considering the relative magnitude of the highest-order coefficient in the Legendre expansion (69) of the density on each panel, $|\hat{\sigma}_{15}|/\mathinner{\!\left\lVert\sigma\right\rVert}_{\infty}$. We take the maximum of this measure over all panels, and denote it the “upsampling error estimate”, shown in Fig. 10. While not a strictly defined error measure, this quantity gives us some indication of how well the density is resolved on the grid. It is clear from the figure that this error estimate starts growing for source points closer than $3h$ from the boundary. Additionally, it seems that the AQBX error grows at the same rate for $d\leq h$, though we can not say for certain that this is the mechanism governing the AQBX error.

The results for this particular test case are encouraging, since AQBX does not exhibit any adverse behavior due to the nearby singularities. Instead, the error increase appears to be due to a lack of resolution.

5.4 Solving the integral equation

The above tests indicate that AQBX works well for evaluating layer potentials, both close to and on the boundary. Since the method is accurate on the boundary, it can also be used to solve the underlying integral equation (2). AQBX is then used to evaluate the left-hand side of the discretized integral equation (i.e. the matrix-vector product), and a solution is found iteratively using GMRES [25]. The details of this are for direct QBX discussed to some extent in [18]. There, they recommend that the principal value integral of the double layer is computed using an average of QBX expansions on both sides of the boundary. This was further studied in [23], where they found that two-sided expansions had better convergence properties. We here follow this recommendation.

Table 3 shows the results when solving our test problem for a range of tolerances. We find that the tolerance set for the GMRES iterations is matched by both the error in the solution $u$ and the smoothness of the density $\sigma$. However, we find that the AQBX tolerance must be set two orders of magnitude smaller than the GMRES tolerance, otherwise GMRES stagnates. We believe that this can partly be explained by the fact that AQBX does not strictly satisfy its given tolerance. Though beyond the scope of the present paper, it would be interesting to further study the interplay between the tolerances of AQBX and GMRES. In particular, we believe that AQBX could be successfully combined with the Inexact GMRES method [26], which is designed to work with a matrix-vector product that has been deliberately made inexact.

6 Conclusions

We have in this paper formulated a scheme for adaptive quadrature by expansion (AQBX), which allows for the evaluation of singular and nearly singular layer potential integrals on a curve discretized using composite Gauss-Legendre quadrature. The scheme automatically sets parameters in order to satisfy a given tolerance. This is a simplification compared to the original QBX scheme, which has a large parameter space. Given a target tolerance, the only free parameter is here the expansion radius $r$. Since the remaining parameters are set on the fly, varying $r$ will mainly affect the speed of the algorithm. The optimal value for $r$ with respect to speed will be implementation-dependent, though values of $W$ in Table 2 suggest that $r/h$ in the range $0.25$–$0.50$ is a good choice.

The key component of our scheme is the ability to accurately estimate the magnitude of the quadrature errors in the QBX coefficients. To do this we have built on the results of [2], where such estimates were reported for a flat panel. We have extended these to curved panels by taking into account the mapping between a flat panel and a curved panel. This mapping can be locally constructed using only the locations of the quadrature nodes on each panel, and therefore requires no additional analytical information. A side benefit of our estimates is that the nearly singular quadrature error of the underlying layer potential can be accurately estimated, which provides an excellent criterion for when to activate special quadrature also when other methods are used [16, 4]. This could also prove useful for QBX schemes where the expansion is formed by multiple layer potential evaluations in a neighborhood of the expansion center [3, 23].

The focus on a target accuracy is, in our experience, uncommon in the context of singular and nearly singular quadrature. We do however believe that this is important if integral equation methods are to be used in large-scale simulations, where the focus is on achieving a target accuracy at the lowest possible computational cost.

While several excellent special quadrature methods exist in two dimensions, methods in three dimensions have not yet reached the same maturity. The QBX method has been successfully used on simple geometries in 3D [1], while development for more general use is ongoing. If accurate parameter selection is important in 2D, it is absolutely essential in 3D, as costs are higher and the impact of suboptimal parameter choices more severe. The principles of the present 2D scheme can be extended to three dimensions. Estimates for the coefficient errors in 3D were in [2] developed for the special case of spheroidal surfaces. Developing estimates for general surfaces in 3D is a topic of ongoing work, and results will be reported at a later date.

Acknowledgments

The authors wish to thank Prof. Johan Helsing for providing an implementation of the singular integration scheme of [14].