Estimation of quadrature errors in layer potential evaluation using quadrature by expansion

In boundary integral methods it is often necessary to evaluate layer potentials on or close to the boundary, where the underlying integral is difficult to evaluate numerically. Quadrature by expansion (QBX) is a new method for dealing with such integrals, and it is based on forming a local expansion of the layer potential close to the boundary. In doing so, one introduces a new quadrature error due to nearly singular integration in the evaluation of expansion coefficients. Using a method based on contour integration and calculus of residues, the quadrature error of nearly singular integrals can be accurately estimated. This makes it possible to derive accurate estimates for the quadrature errors related to QBX, when applied to layer potentials in two and three dimensions. As examples we derive estimates for the Laplace and Helmholtz single layer potentials. These results can be used for parameter selection in practical applications.


Introduction
At the core of boundary integral equation (BIE) methods for partial differential equations (PDEs) lies the representation of a solution u as a layer potential, where Γ is a smooth, closed contour (in R 2 ) or surface (in R 3 ), σ(y) is a smooth density defined on Γ, and G(x, y) is a Green's function associated with the current PDE of interest. The Green's function is typically singular along the diagonal x = y, which leads to the integrand being singular for x ∈ Γ. We will refer to this as a singular integral. The field that it produces is however smooth on each side of Γ. When x is close to Γ, but not on it, the layer potential is difficult to evaluate accurately using a numerical method, even though G is smooth. In this case we say that the integral is nearly singular, since we evaluate G close to its singularity. In this paper we discuss the use of residue calculus for estimating the error committed when computing a nearly singular integral using a quadrature method. The error estimates are derived in the limit n → ∞, n being the number of discrete quadrature points, but turn out to be accurate also for moderately large n. Throughout we shall use the symbol to mean "asymptotically equal to", such that The discussion is limited to the Gauss-Legendre rule and the trapezoidal rule, which are perhaps the two most common quadrature rules in the BIE field. The kernels that we consider are related to a new method for singular and nearly singular integration, called quadrature by expansion (QBX) [11,2,8]. Specifically, we consider two classes of kernels that appear when applying QBX to the single layer potential of Laplace's equation in two and three dimensions, also referred to as the single layer harmonic potential, This paper is organized as follows: In section 2 we introduce QBX and the relevant kernels for studying the quadrature errors associated with the method. In section 3 we discuss the required framework for estimating quadrature errors, and use it to derive error estimates for our kernels. Finally, in section 4 we show how our results can be used to compute quadrature error estimates that are useful in practical applications. The results provided in section 4 are for model geometries; in appendix A we include results for a more complex geometry.

Quadrature by expansion (QBX)
The central principle of QBX [11] is the observation that a layer potential is smooth away from Γ, such that it locally can be represented using a Taylor expansion around an expansion center x 0 . Given a quadrature rule and a target point x on Γ, the method can be summarized as: 1. Pick an expansion center x 0 at a distance r from x in the normal direction, such that x 0 lies in the quadrature rule's region of high accuracy.
2. Compute a local expansion of the potential around x 0 , truncated to some expansion order p.
3. Evaluate the local expansion at x.
Γ x0 x r Figure 1: QBX geometry. The local expansion formed at x 0 is valid inside the ball of radius r and at x.
The expansion is convergent inside the ball of radius r centered at x 0 , and at the point of intersection of the ball and Γ [8], as illustrated in Figure 1. We can therefore use QBX to compute the potential both when it is singular and nearly singular, as long as our target point lies inside the domain of convergence of a local expansion.
Two-dimensional single layer potential In two dimensions it is convenient to let R 2 = C, and work with the complex version formulation for the single layer potential, u(z) = Re Γ σ(w(t)) log(z − w(t))dt, where w(t) is an arc length parametrization of the contour Γ ∈ C. Forming the local expansion at a center z 0 using the series expansion of log(1 − ω), |ω| < 1, we get [8] where the expansion coefficients a j are given by Truncating the expansion to order p and denoting byã j the coefficients computed using a quadrature rule, we define our QBX approximation of the potential as Three-dimensional single layer potential In three dimensions an expansion of the Green's function about a center x 0 is formed using the spherical harmonic addition theorem [9], where Y m l is the spherical harmonic of degree l and order m, Here P m l is the associated Legendre function, and (θ x , ϕ x ) and (θ y , ϕ y ) are the spherical coordinates of x − x 0 and y − x 0 . The local expansion is then with the expansion coefficients α m l given by Again truncating the expansion at order p and lettingα m l denote coefficients approximated by quadrature, we get our QBX approximatioñ

Error analysis of QBX
The error when computing a layer potential using QBX can be divided into two components: the truncation error and the quadrature error. The truncation error comes from limiting the local expansion to a finite number of terms, and was analyzed by Epstein et al. [8]. The quadrature error comes from evaluating the integrals (7) and (12) using a discrete quadrature rule. To see the separation of errors we add and subtract the exact expansion coefficients to the QBX approximation. For the single layer in two dimensions we then get the separation as 4 In [8,Thm 2.2] it is shown that there is a constant M p,Γ such that if σ ∈ C p (Γ), then where r is the distance from Γ to the expansion center. In three dimensions we similarly have (assuming for the moment x 0 = 0) Here a generalization of the results in [8,Thm 3.1] gives that there is a constant M Γ,δ , δ > 0, such that if σ belongs to the Sobolev space H 3+p+δ (Γ), then If we let h be a characteristic length of the discretization of Γ, and furthermore keep the ratio r/h constant under grid refinement, then a consequence of the truncation error estimates is that e T converges with order p + 1 under refinement, In [8] and [11] it is argued that if the quadrature error can be maintained at a fixed level , then such that QBX is convergent with order p + 1 until hitting the "floor" given by .
Turning to the quadrature error e Q , one can in practical applications observe that it grows with the expansion order p, such that there for a given problem exists an optimal p where the total error e T + e Q has a minimum. To control e Q (and thereby the minimum error) one can interpolate σ to a finer grid before computing the coefficients, a technique usually referred to as "upsampling" [2] or "oversampling" [11]. This works because at large p the difficulties lie in accurately resolving the kernels in (7) and (12), not σ itself.
The quadrature error e Q has for the two-dimensional case been discussed in [8] and [2]. In [8] an upper bound, which did not explicitly show the dependence on p, was derived for the case when Γ is discretized using Gauss-Legendre panels. In [2, Thm 3.2], a bound including the p-dependence was derived for Γ discretized using the trapezoidal rule.
In what follows of this paper we shall develop quadrature error estimates, not bounds, that accurately predict the order of magnitude of |a j −ã j | and |α m l −α m l | in the asymptotic region of n → ∞. We will do this for the n-point trapezoidal and Gauss-Legendre quadrature rules on certain geometries. To that end, we will consider two different classes of kernels (depicted in Figure 2), where the respective singularities z 0 and (x 0 , y 0 ) of the kernels are assumed to lie close to, but not on, the interval of integration. The kernel f p is relevant when analyzing the computation of a j , which we will refer to as the complex kernel. The corresponding kernel for α m l is g p , and we will refer to it as the Cartesian kernel.

Estimating quadrature errors
The purpose of quadrature is to numerically approximate the definite integral of a function f over some interval Γ, where Γ is typically a subset of the real line or a closed curve in the complex plane. This approximation is computed through an n-point quadrature rule, using a set of nodes x i ∈ Γ and weights w i , For an n-point interpolatory quadrature rule, the weights w i are determined by integrating the polynomials of degree n − 1 that interpolate f at the nodes x i . The error of the approximation is given by the remainder term which converges to zero as n → ∞ unless the function is ill-behaved in some way, e.g. if f has a singularity on Γ. As practitioners of quadrature, we typically want to know the relationship between R n and f in order to make our computations efficient and reliable. Luckily, most quadrature rules come equipped with a number of error bounds involving f in one way or another. However, most such bounds are useful only if f is analytic in the whole complex plane, or in a large region around Γ. In section 3.1.1 we give an example where an error bound for the Gauss-Legendre rule wildly overestimates the actual error.
We now consider what happens when f is meromorphic, i.e. analytic everywhere except for at a set of poles. The focus of our present work is the case when f has poles close to Γ. Then the best way of obtaining accurate estimates of R n is to use contour integrals in the complex plane, using the theory of Donaldson and Elliott [6]. The key principle is that the quadrature rule Q n can be connected to a rational function q n (z), which has simple poles at the quadrature nodes x i , and residues at those poles equal to the quadrature weights w i . If C is a contour enclosing Γ, on and within which the complex continuation of f is analytic, then There also exists a characteristic function m(z) such that we can express the integral over Γ as a contour integral, We define the remainder function which is analytic in the complex plane with Γ deleted. Using k n we can express the remainder as a contour integral, As |z| tends to infinity, k n tends to zero at least like O(|z| −n ) for interpolatory quadrature [7], so by taking C to infinity we see that interpolatory quadrature is exact for polynomials of degree < n. If f has one or more poles z j enclosed by C, we deform the contour integral to include small circles enclosing those poles (see e.g. the illustration in [7]). Letting the radius of the circles go to zero, the remainder is given by the integral over C minus the residues at the poles, For reference we here state the definition of the residue of a function g(z) which has an order p pole at z 0 : If the poles of f are close to Γ, then the remainder R n [f ] is typically dominated by the corresponding residues, and the contribution from the contour integral is negligible. If f (z)k n (z) goes to zero as C is taken to infinity for some n ≥ N , then the remainder is equal to the sum of the residues for those n.
Given a meromorphic integrand f and a quadrature rule Q n , our ability to estimate R n depends on our knowledge of k n (z) and our ability to compute the residues of f (z)k n (z) at the poles. For some quadrature rules we have closed form expressions for k n (z), while for others we only have asymptotic estimates valid for large n. In what follows we shall summarize the relevant formulas for the Gauss-Legendre and trapezoidal quadrature rules and apply them to our model kernels.

The Gauss-Legendre quadrature rule
The n-point Gauss-Legendre quadrature rule [1, ch. 25] belongs to the wider class of Gaussian quadrature, and is extensively used in applications where the integrand is not periodic. It is by convention defined for the interval [−1, 1], The weights and nodes of the quadrature rule Q n are associated with P n (x), the Legendre polynomial of degree n. The nodes are the roots of P n , P n (x i ) = 0, i = 1, . . . , n, and ordered, such that x i < x i+1 . The weights are given by In our analysis of the kernels f p and g p we will stay on the standard interval [−1, 1] on the real axis. Setting our target kernels are then

Classic error estimate
There exists a classic error estimate for the Gauss-Legendre rule, available in e.g. Abramowitz and Stegun [1, eq. 25.4.30], stating that on an interval of length L the error is given by The most important interpretation of this result is that n-point Gauss-Legendre quadrature will integrate polynomials of degree up to 2n − 1 exactly. In practice the estimate is only useful for very smooth integrands, due to the high derivative of f in the estimate.
Consider the example of f = g 1 with a = 0, which can be found in Brass and Petras [5], The norm of the derivative is given by Inserting this into the estimate and applying Stirling's formula, we get that This bound goes to infinity exponentially fast for b < 1/2, while in practice Gauss-Legendre quadrature exhibits exponential convergence for this integrand. There exists a large number of error estimates involving lower order derivatives of the integrand, available in the work by Brass et al. [4,5]. Alternatively, one can estimate the error through a contour integral, which is what we will do in the following section.

Contour integral
An expression for estimating the error of Gaussian quadrature as a contour integral was found by Barrett [3] for certain cases, and later generalized by Donaldson and Elliott [6]. The below derivation follows that of Barrett [3].
It can be shown for Gauss-Legendre quadrature that the weights are given by where x i are the zeros of P n (x). The weights can also be computed as the residues of the function at the nodes x i , It follows that where C now contains [−1, 1]. It can also be seen that since we can write where It follows that the remainder function k n (z) in (28)-(30) is given by While we do not have a closed form expression for k n (z), it can in the limit n → ∞ be shown to satisfy [3,6] k n (z) where c n = 2π(Γ(n + 1)) 2 Γ(n + 1/2)Γ(n + 3/2) 2π.
While this is an asymptotic result valid for n → ∞, we shall see that it provides an accurate approximation of k n also for moderately large n.
In the following analysis we will need derivatives of k n , for the residue at high order poles. The first two derivatives are From (54) and (55) we can induce that the qth derivative of k n will have the form such that we for large n can use the approximation

Complex kernel
We now wish to study the Gauss-Legendre rule applied to the complex kernel (37) on The integral is then which has a pole of order p at z 0 . Letting C in (30) enclose [−1, 1] and z 0 , the integrand f p (z)k n (z) vanishes as we take C to infinity. The remainder is given by the residue, Using the estimate (57) for the derivatives of k n , we get the following estimate for the remainder: The magnitude of the remainder when using the n-point Gauss-Legendre rule to compute 1 −1 f p (x)dx is then in the asymptotic limit n → ∞ given by Proof. The proof follows from (52), (57) and (59).
We have (see e.g. [3]) that the factor |z + √ z 2 − 1| 2n+1 is constant on an ellipse with foci at ±1, and consequently that it has a minimum when Re z = 0. The decay rate for a general z 0 = a + ib is therefore bounded by the decay rate given by z 0 = ib, Assuming b 1 and discarding O b 2 terms, we can simplify the result of Theorem 1 to where a(n) b(n) now denotes "approximately less than or equal to" in the limit n → ∞, in the sense that for a K(n) such that a(n) ≤ K(n), then lim n→∞ K(n)/b(n) ≈ 1.
Although just an approximation, the result in (62) compares very well to numerical experiments, as shown in Figure 3.

Cartesian kernel
Let us now consider the Cartesian kernel (38) and the quadrature error when computing using Gauss-Legendre quadrature. The case of p = 1 was studied by Elliott, Johnston & Johnston [7], and following their example we write the complex extension of g p as which has poles at z 0 and z 0 , Letting the contour C enclose [−1, 1] and the poles, we see that the integrand vanishes as we take C to infinity, and we are left with a remainder determined by the residues, with k n as defined in (52). These residues are easily evaluated as (67) Carrying out the full computations for p = 1, 2, 3, we get For large n, the dominating term will be the one with the highest k n -derivative, so we can make the approximation We combine (57) denote the remainder when using the n-point Gauss-Legendre rule to compute 1 −1 g p (x)dx. Then, in the asymptotic limit n → ∞, where z 0 = a + ib and Proof. The proof follows from (66), (71) and (57).
Repeating the argument of section 3.1.3, we can by assuming b 1 estimate the largest error for all a as which provides a relatively clear view of how the error depends on p, b and n. Figure  4 shows experimental results using both this estimate and Theorem 2, as well as the results obtained when including all terms of the differentiation in (67). The left column of Figure 4 shows results for a = 0, in which case the error has an oscillatory behavior in n. These oscillations are bounded by the case a = 0, shown to the right. The top row shows results for p = 1, where our estimates are very accurate. The center row shows that our estimates lose accuracy for small n at p = 5, though that loss can be recovered by including all derivatives, as in the bottom row.  67), "Estimate" using (72) and "Simple estimate" using (74). 14

The trapezoidal rule
For periodic integrands, the trapezoidal rule is typically the quadrature rule of choice. It is easy to implement, has an even point distribution and converges exponentially fast. In our analysis we will assume the interval of integration to be the unit circle, parametrized in an arc length parameter t. Without loss of generality we can then assume the singularity to lie on the x-axis at a distance b from the boundary, with b > 0. Our target kernels are then The kernel f 1 corresponding to p = 1 is essentially that of the Laplace double layer potential. A bound on the quadrature error for this potential was derived in [2, Thm 2.3] for Γ a general curve discretized using the trapezoidal rule. This bound corresponds to the result in Theorem 3 (below) for p = 1 and Γ the unit circle.

Contour integral
We here state the necessary results for error estimation using contour integrals for two cases: integral over a periodic interval and integral over a circle in the complex plane. These results and more can be found in the thorough review by Trefethen and Weideman [14]. It is worth noting that the remainder functions (79)  f (x k ), x k = 2πk/n.
To compute R n [f ], we let C be the rectangle [0, 2π] ± ia, a > 0, traversed in the positive direction. The sides of rectangle cancel due to periodicity, so we need only consider the top and bottom lines. The appropriate remainder function is given by [14] k n (z) = 2πi Integral over circle in the complex plane For a function f (z) on the unit circle we have z = e it , such that and For the contour integral we let C be the circle |z| = r > 1, and the remainder function is then given by [14] k n (z) = −2π z(z n − 1) .

Complex kernel
Integrating f p (76) on the unit circle in the complex plane, we wish to compute with z(t) = e it and z 0 = 1 + b. Letting C enclose the unit circle and z 0 , we can compute the remainder of the trapezoidal rule as with k n as defined in (82). Taking C to infinity the integral vanishes, and the error is given by the residue at z 0 , where there is a pole of order p such that If we evaluate the derivative analytically, this expression is exact. For large n we can estimate k n as, such that We can simplify the product in the numerator through (n + 1) · · · (n + p − 1) (n + p) p−1 .
Putting it all together, we get the following result: be the quadrature error when computing 2π 0 f p (e it )dt using the n-point trapezoidal rule. In the limit n → ∞ we then have that Proof. The proof follows from (85), (87) and (88).
The expression in Theorem 3 gives a very good estimate of the error. In Figure 5 we compare it to numerical results for p = 1 and p = 10, and in both cases it captures the region of exponential convergence very well.

Cartesian kernel
Integrating g p (77) on the unit circle with (x 0 , y 0 ) = (1 + b, 0), our target integral is where t is an arc-length parameter describing the unit circle. Letting g p (z) be the complex extension of g p (t), we can write This function has poles at z 0 and z 0 , It is also periodic on the interval [0, 2π], so we let C be the rectangle [0, 2π] ± ia, a > 0, traversed in the positive direction, and take the remainder function as defined in (79).
Letting a go to infinity the contribution from the contour vanishes, and the remainder is given by the residues at the poles, To compute the residue at z 0 we begin with the definition (31) Taking the limit for the first few p we note that the residues from z 0 and z 0 are equal, such that we can express the remainder as When n is large we can estimate the remainder function as k n (z) 2πi −e inz Im z > 0, such that For large n the remainder will be dominated by the highest derivative of k n , so we can simplify to get the following result: be the quadrature error when computing 2π 0 g p (cos t, sin t)dt using the n-point trapezoidal rule. For n → ∞ the error is then asymptotically given by  The exponential convergence is well captured, though at low n and large p some accuracy is sacrificed by our simplifications. If the point (x 0 , y 0 ) is not on the x-axis, but still on the circle of radius 1 + b, then the error has wiggles similar to those present in Figure 4. The result in Theorem 4 then follows the maximum of those wiggles.

Comments
We now have asymptotically accurate estimates of R n [f p ] and R n [g p ] for the trapezoidal and Gauss-Legendre rules on simple model geometries (unit circle and line segment). The derivations of these four estimates all follow the same basic recipe, which for an integrand f and remainder function k n can be summarized as follows: (i) Find the poles {z i } of f .
(ii) Take the error to be the residues of k n f at {z i }.
(iii) If necessary, simplify the result by keeping only the term that dominates as n → ∞.
In the examples we have considered, this has been the one with the highest derivative of k n .
A few comments are also in order before we move on to applying our results to more general problems.
Non-integer poles Our derivations of error estimates for f p and g p are only valid for integer p. In practice we want to estimate the error for kernels like |x − x 0 | −q , q ∈ N, which means that p also takes half-integer values. Instead of carrying out new derivations for p half-integer, we simply note that the estimates we already have work well in practice also for half-integers, as long as the factorial is computed using the gamma function, which is true for integer p.
Taking the density into account When computing the QBX coefficients (7) and (12), we typically have kernels like f p or g p times a non-constant density σ. For the 2D coefficients a j in (7) we for example have the integrand σ(z)f j (z). With a pole of order j in f j at z 0 the residue contains all the derivatives of σ up to σ (j−1) evaluated at z 0 , so a minimum requirement is that those derivatives are bounded. Let us now assume σ to be smooth everywhere, which in the BIE setting is quite reasonable. The only residue to consider is then the one at z 0 , which is dominated by the highest derivative of k n , Depending on σ, the error may also have a contribution from the contour C, as the contour integral may be non-vanishing at infinity. This contribution can however be safely neglected in the asymptotic region, as the contribution from the poles then dominates the error, at least for poles close to the domain of integration. As a concrete example, let us now consider the density σ(x) = x k e imx multiplied with g p , integrated on [−1, 1] using Gauss-Legendre quadrature. (A density like this with k or m nonzero would appear when using a discretization that corresponds to a polynomial or a Fourier series that is integrated term by term.) The integrand is analytic everywhere expect at the poles z 0 = a + ib and z 0 , so we can compute the error as the contour integral of k n h over C plus the residues. As C tends to infinity the integral does not vanish for any n, since then |k n (z)h(z)| |z| k−2p−2n−1 e m|z| . Taking C to be a circle of radius R 1, Minimizing the rightmost bound with respect to R gives R min = 2n−k m , such that We have σ smooth everywhere, so we use the simplification (102) of keeping only the highest order derivative in k n , such that The last term has faster than exponential decay, so the contribution from the poles will dominate the error. Inserting (71), We can thus get an error estimate for h by taking our previous results for the kernel f p , multiplied by the density σ evaluated at the poles of f p . Our experience is that the above reasoning holds true also in the general case, such that we can estimate the quadrature error for an arbitrary, smooth density using the error estimate for the kernel times the density evaluated at the poles. In practice we can simplify this even further by using the values of the density on Γ, since that is what we have access to in a numerical implementation. If x c is the point on Γ closest to z 0 , we then use that σ(z 0 ) ≈ σ(x c ) if r small and σ smooth. For a nearly singular kernel f this allows us to use existing error estimates by writing (109)

Applications
In the previous section we showed how to develop quadrature error estimates for the kernels f p and g p , when integrated on model geometries using the trapezoidal and Gauss-Legendre quadrature rules. We will in this section show how these error estimates can be used to estimate the quadrature errors of QBX. Before we do that, however, we will show an example of how we can use our results to estimate the error when evaluating a nearly singular layer potential using regular quadrature.

Double layer potential in two dimensions
To see how our results can be used when working with the double layer potential, we now consider an example from Helsing & Ojala [10, sec. 10.1]. We then consider the two-dimensional double layer potential in complex form where σ is the solution to an interior Dirichlet problem with a known reference solution 1 u ref (z). The boundary is starfish-shaped with parametrization and is divided into 35 panels Γ i of equal length in t. Each panel is discretized using a 16point Gauss-Legendre quadrature rule, for a total of 560 discretization points. Computing u(z) in the interior using the Gauss-Legendre quadrature results in large errors close to the boundary. This can be seen in Figure 7a, which shows the relative pointwise error An accurate estimate of e(z) for Γ discretized using the trapezoidal rule can be observed in [2, Fig. 1]. To estimate e(z) when using Gauss-Legendre panels, we can use our results from section 3.1.3 to estimate the error e i from each panel Γ i , and then sum them together, (In practice it suffices to use the contribution from the two closest panels, as the closest panel completely dominates the error except for when z is close to the edge between two panels.) Replacing each panel with a corresponding flat panel, we can generalize the results of Theorem 1 to estimate the error from panel Γ i as where z 0 is the location of the pole z under the transformation that takes Γ i to [−1, 1]. We approximate the imaginary part of z 0 as Im z 0 = d/L, where L is the length of Γ i and d is the shortest distance from z to Γ i . The real part Re z 0 is approximated as the real part of z after applying the scaling and rotation that takes the endpoints of Γ i to −1 and 1.
Evaluating the estimate (114) in the interior produces an error plot which in "eyeball norm" is identical to that in Figure 7a. If we are more careful and compare the level sets of the error and the estimate (Figure 7b), we see that the correspondence is extremely good for panels that are close to flat, while the estimate suffers some inaccuracy for curved panels. Increasing the number of panels improves the accuracy of the estimate ( Figure  8), as the individual panels then are less curved. By removing the absolute value in the denominator of (114) and instead taking the imaginary part of the whole expression we can also reproduce the small-scale oscillations of the true error, though that has small practical relevance.    Figure 8: Same plots as in Figure 7, but with 70 Gauss-Legendre panels.
These results suggest that our error estimates are useful for estimating the magnitude of quadrature errors due to the near singularity of the integral. This allows for a cheap way of determining when one needs to use a special quadrature method, such as QBX or the one outlined in [10].

QBX in two dimensions
Let us now return to the QBX quadrature error for the single layer potential in two dimensions, which we introduced in section 2. We let Γ be a curve divided into N panels of length h, and compute the coefficientsã j at an expansion center z c from (7) using n-point Gauss-Legendre quadrature on each panel. We then have from (14) that the quadrature error is which is bounded by This error is discussed in Epstein et al. [8], where they use the standard Gauss-Legendre error estimate (39) to get from which it is concluded that r > h/4 is required for the quadrature to converge. This does however not give any information about the rate at which the quadrature error grows with p. Nor does it give any practically useful information about the quadrature error, since it is easy to numerically verify that r < h/4 works just fine as long as n is large enough. In fact, the standard Gauss-Legendre error estimate is overly pessimistic for this type of integrand, as discussed in section 3.1.1. Using the result in (62), as follows from Theorem 1, and generalizing as discussed in section 3.3, we can estimate the quadrature error on the interval [−1, 1] for a function of type h(z) = σ(z)(z − z 0 ) −j as under the assumption Im z 0 1 and σ smooth. Using a change of variables from an interval of length h to [−1, 1] and setting z 0 = 2ir/h, we can estimate the quadrature error for a single coefficient as where Γ r is the strip of width r stretching from Γ into the domain. In practice we can use σ L ∞ (Γr) ≈ σ L ∞ (Γ) if r small and σ smooth. For Γ a perfectly flat panel the constant would be C(Γ, h) = 2π, and in practice C is close to 2π if the panels on a general Γ are close to flat, since the error is dominated by that from the panel closest to z c . Putting (116) and (119) together allows us to estimate the quadrature error as Discarding the term corresponding to j = 0 (which is small) and using Stirling's formula (42), we can simplify this to This expression gives a reasonably clear of view of how the error depends on the involved variables, and Figure 9 shows that it captures the behavior of the quadrature error quite well. The quotient 4r/h appears here too, but without any type of bound; the convergence in n will just stall as r/h → 0. Interestingly, the estimate (121) is similar to the bound on the quadrature error derived in [2, Thm 3.2] for the double layer potential using the trapezoidal rule. The sum over j in (120) and (121) makes the expressions a bit cumbersome to interpret, though we can deduce that the error will in some parameter regions grow exponentially in p. Recognizing that the sum in (120) is the truncated exponential sum, we can write such that, surprisingly, the quadrature error has an approximate upper bound independent of p, Alternatively, we can use the definition of the incomplete gamma function Γ(n, x) for integer n, to put the estimate of (120) in a form that may not be easier to interpret, but at least is more compact,   (14) of 2D QBX on the unit circle, using density σ = (sin θ) 10 , measured at all nodes by comparing to a highly resolved reference solution. Numerical setup is N = 20, n = 100 and r/h = 0.1. The quadrature error e Q and its upper bound are estimated using (120) and (123) with C = 2π, while the decay of the truncation error e T roughly follows (15).

QBX in three dimensions
For the single layer potential in three dimensions, we need to estimate the quadrature error in (12). The surface quadrature on Γ is assumed to be a tensor product quadrature rule, such that the surface is sliced into one-dimensional cross sections. We will here show how to develop estimates for a simple square Gauss-Legendre patch, but results for the more complicated case of a spheroid can be found in appendix A. Before doing that, however, we need to work out two preliminaries: (i) How to relate the remainder of a composite surface quadrature to the remainders of the corresponding one-dimensional quadrature rules. (ii) How to account for the magnitude and nearly singular behavior of the spherical harmonics component of the kernel in (12).
Surface quadrature remainder Let I 2 denote a composite surface integral that is independent of the order of integration, and let the subscript I s denote an integral carried out in the variable s, Applying a tensor product quadrature rule over the surface, we get Assuming the quadratic remainder term to be negligible and that Q n,α ≈ I α , we can approximate the remainder of the surface quadrature as So the surface remainder is approximately equal to the integrals of the one-dimensional remainders, which is what one would expect.
Spherical harmonics kernel For 3D QBX, our goal is to compute the quadrature error e Q (17), where We can simplify this by using the Legendre polynomial addition theorem [9], where θ is the angle between x and y. Using this, we can simplify the quadrature error to |x − x 0 | l R 2 n,y |y − x 0 | −l−1 P l (cos θ)σ(y) .
Our task is then to estimate the quadrature error of the kernel which we shall refer to as the Legendre kernel. We are mainly interested in the quadrature error when QBX is used for singular integration, so we will assume x ∈ Γ, such that In order to estimate R 2 n [K l ], let us now consider the integral along a curve that is the intersection of Γ and a plane containing the expansion center x 0 = (x 0 , y 0 , z 0 ). We name the curve γ and assume without loss of generality that it lies in the xz-plane, such that y = (x, 0, z) and x 0 = (x 0 , 0, z 0 ). Then, Also assuming x 0 − x = rẑ, we have that The Legendre polynomials P l (cos θ) are degree l polynomials in cos θ, and the quadrature error will be determined by the leading coefficients of those polynomials, for reasons that will become clear shortly. From Rodrigues' formula we can derive that Inserting (136), we get Now we see the reason why the leading term in the polynomial will dominate the quadrature error; cos θ has a pole at x 0 , and the highest order of that pole will dominate the quadrature error. The quadrature error for K l can therefore be found by studying the quadrature error of the function since on γ The magnitude of the spherical harmonics B l can be simplified using Stirling's formula n! ≈ √ 2πn n+ 1 2 e −n , An error estimate for ψ l follows from our results for the Cartesian kernel g p (22) in Theorems 2 and 4, since In fact, a good estimate of the error is obtained by analyzing the simpler form where r = min x∈Γ |x − x 0 |. Figure 10: An n × n Gauss-Legendre patch.

Gauss-Legendre patch
When forming a quadrature rule for a general surface Γ, a straightforward method that is often used in BIE methods is to divide the surface into approximately square patches, and then use an n × n tensor product Gauss-Legendre rule on each patch (which then looks like Figure 10). Denoting the patches Γ i , the QBX expansion coefficients at a center x 0 are then computed as with the approximate coefficientsα m l computed using the associated quadrature rule for each patch. To be able to estimate the QBX quadrature error e Q (17) we need to be able to estimate the quadrature error from each patch when evaluating (146) using Gauss-Legendre patches.
We now focus on the error from a single Gauss-Legendre patch, in the special case of it being square and flat. For that we let Γ be the patch (x, y, z) ∈ [−1, 1] × [−1, 1] × {0}, and let x 0 = (x 0 , y 0 , r) be a point close to Γ. We consider the simplified form of the quadrature error (132), where K l is the Legendre kernel (133). To estimate the quadrature error of K l on the patch, we consider the quadrature error of the equivalent kernel ψ 2 l , ψ 2 l (x, y) = B l r l g 2 which is the 2D analogue of (145) and satisfies R 2 n K l ≈ R 2 n ψ 2 l . Here g 2 p is the Cartesian kernel over the patch, defined as Evaluating the integral of g 2 p on Γ using the tensor product Gauss-Legendre rule, we can expect (and verify) that the error for a given r is largest for x 0 above the center of the patch, x 0 = y 0 = 0. By symmetry we then have that I x R n,y = I y R n,x , so we can use the results in (74) and (101) with b 2 = y 2 + r 2 to estimate R n,x [g 2 p ], which gives us We compute an approximation of this integral by expanding the square root around y = 0, e −2n(r+y 2 /(2r)) dy = r −p e −2nr πr/n erf( n/r).
We are considering large n and small r, so erf( n/r) ≈ 1, such that This means that the result of the integration in y is approximately equal to multiplying the value of the integrand at y = 0 with δ = πr/n, independent of the integration bounds (as long as the interval is wider than δ). An interpretation of this is that most of the error comes from a strip of width δ centered around y = 0. Finally combining (152), (150) and (148), we get the estimate for ψ 2 l on a patch, We now consider a slightly more general case, where the patch is of size h × h. The corresponding change of variables from a unit square patch in the integral of ψ 2 l allows us to use the results in (153) with the additional factor (2/h) l−1 and the change r → 2r/h. It follows that the quadrature error for the Legendre kernel K l on a flat Gauss-Legendre patch with sides h can be estimated as where r is the distance from the expansion center to the patch. We now let x t = (x t , y t , 0), −1 ≤ x t , y t ≤ 1, be a target point on the patch Γ and x 0 = (x, y, r) be the corresponding expansion center. Using (109) gives us Inserting this, (141) and (154) into (147) gives us the final expression for the QBX quadrature error,  Figure 11: The error components (16) and (17) of 3D QBX when applied to a Gauss-Legendre patch and compared against a reference solution. The quadrature error e Q is well approximated by the estimate (156), while the decay rate of the truncation error e T follows (18).
The accuracy of this estimate is demonstrated in Figure 11, which shows the QBX errors for x t at the center of the patch, x t = (0, 0, 0). The grid has 96 × 96 points, r = 0.2 and σ is a two-dimensional polynomial of degree 15 with random coefficients. This would correspond to using a 16th order patch with factor 6 oversampling.
To put our result in (156) in a more general form, we simplify it one step further. From Stirling's formula (42) we have that Γ(l + 1/2) ≈ Cl l e −l . Combining this with (143) and discarding the l = 0 term allows us to estimate the 3D QBX quadrature error on a Gauss-Legendre patch as This is completely analogous to the 2D QBX result for Gauss-Legendre panels (121).

Complex geometries
So far we have developed error estimates for simple geometries, but the framework can be extended to more complex geometries in a straightforward fashion, though the computations may be cumbersome. As an example of how one can proceed, in appendix A we develop QBX quadrature error estimates for Γ being a spheroid, discretized using both the trapezoidal and Gauss-Legendre quadrature rules. Another possibility would be to extend the previous section's estimates for the Gauss-Legendre patch, to account for the patch having curvature. This might be useful when estimating quadrature errors on a general surface which has been divided into patches.

Conclusions
The model kernels which we have considered, f p (z, w) = |z − w| −p and g p (x, y) = |x − y| −2p , can be found in two and three dimensional BIE applications in general (for small p), and in QBX in particular. Using the method of contour integrals, we have in section 3 derived accurate estimates (Theorems 1-4) of the quadrature errors when integrating these kernels close to a singularity using the n-point Gauss-Legendre and trapezoidal quadrature rules (on [−1, 1] and the unit circle, respectively). These estimates are not in the form of bounds, which is what one classically seeks in numerical analysis. Instead, they are asymptotic equalities valid in the limit n → ∞. Their key feature is that they predict the magnitude of the errors surprisingly well also for small n, as we have seen throughout this paper. Extracting the dominating features of our estimates, we arrive at the summary presented in Table 1, which shows how the quadrature errors depend on n, b and p.   By applying suitable generalizations, we have in section 4 of this paper demonstrated how to use our results when working with BIEs and QBX. These results are approximations rather than asymptotic equalities, but nevertheless provide error estimates which are accurate enough to be used for parameter selection in practical applications. Explicit estimates have been developed for Gauss-Legendre panels in two dimensions, and for flat Gauss-Legendre patches and spheroids (appendix A) in three dimensions. Though these results may prove useful by themselves, the more important result is that the methodology used to derive them can be generalized to other kernels and geometries in a straightforward manner.
The QBX quadrature error e Q for the Laplace single layer potential in two (14) and three (12) dimensions is well captured by our estimates (121) and (157), when using Gauss-Legendre as the underlying quadrature. For the Helmholtz single layer potential in two dimensions, the corresponding estimate is (170). A common feature of these estimates is that most of their dependence on the parameters 2 n, p, r and h is captured This is in turn very similar to the results for the spheroid (213) derived in appendix A, and to the bound for the double layer potential in two dimensions derived by Barnett [2,Thm 3.2]. Taken together, these expressions provide a better understanding of the QBX quadrature error, which appears to have similar behavior independent of kernel and dimension. Previous results by Epstein et al. [8] for the truncation error (15,18) provide insight into the truncation error and establish the analytic foundation for the method.
Putting these together with the results presented here, it is now possible to understand the complete error spectrum when working with QBX both in two and three dimensions. We have in this paper focused on estimates for the harmonic single layer potential in two and three dimensions. We have also shown how to derive an estimate for the Helmholtz single layer potential in two dimensions, and derivation of similar estimates for other kernels should follow along the same lines.