Quadrature by Two Expansions: Evaluating Laplace Layer Potentials using Complex Polynomial and Plane Wave Expansions

The recently developed quadrature by expansion (QBX) technique accurately evaluates the layer potentials with singular, weakly or nearly singular, or even hyper singular kernels in the integral equation reformulations of partial differential equations. The idea is to form a local complex polynomial or partial wave expansion centered at a point away from the boundary to avoid the singularity in the integrand, and then extrapolate the expansion at points near or even exactly on the boundary. In this paper, in addition to the local complex Taylor polynomial expansion, we derive new representations of the Laplace layer potentials using both the local complex polynomial and plane wave expansions. Unlike in the QBX, the local complex polynomial expansion in the new quadrature by two expansions (QB2X) method only collects the far-field contributions and its number of expansion terms can be analyzed using tools from the classical fast multipole method. The plane wave type expansion in the QB2X method better captures the layer potential features near the boundary. It is derived by applying the Fourier extension technique to the density and boundary geometry functions and then analytically utilizing the Residue Theorem for complex contour integrals. The internal connections of the layer potential with its density function and curvature on the boundary are explicitly revealed in the plane wave expansion and its error is bounded by the Fourier extension errors. We present preliminary numerical results to demonstrate the accuracy of the QB2X representations and to validate our analysis.



page 4

page 15

page 16

page 17

page 18

page 19


High-order close evaluation of Laplace layer potentials: A differential geometric approach

This paper presents a new approach for solving the close evaluation prob...

On the regularization of Cauchy-type integral operators via the density interpolation method and applications

This paper presents a regularization technique for the high order effici...

On the Approximation of Local Expansions of Laplace Potentials by the Fast Multipole Method

In this paper, we present a generalization of the classical error bounds...

Isoparametric singularity extraction technique for 3D potential problems in BEM

To solve boundary integral equations for potential problems using colloc...

Fast Computation of Electromagnetic Wave Propagation and Scattering for Quasi-cylindrical Geometry

The cylindrical Taylor Interpolation through FFT (TI-FFT) algorithm for ...

Quadrature error estimates for layer potentials evaluated near curved surfaces in three dimensions

The quadrature error associated with a regular quadrature rule for evalu...

Efficient Fast Multipole Accelerated Boundary Elements via Recursive Computation of Multipole Expansions of Integrals

In boundary element methods (BEM) in ℝ^3, matrix elements and right hand...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

When the integral equation method is applied to solve a given partial differential equation, one numerical challenge is the accurate and efficient evaluation of the singular, weakly singular, or hyper singular integrals representing different potentials in the integral equation reformulations. For example, the solutions of a homogeneous elliptic equation (e.g., Laplace, Helmholtz, or Yukawa equations) with different types of boundary conditions are often re-expressed as combinations of the single layer and double layer potentials with density functions and


where is the free-space Green’s function for the underlying elliptic PDE, is the source point located on the boundary , is any target point located in the computational domain, and

is the outward normal vector at

. The Green’s function is usually a smooth function when is away from on the boundary , but becomes singular when . Therefore, different numerical strategies have to be designed for cases when is far away from the boundary, exactly on the boundary, and close to the boundary.

The research topics of developing different numerical integration schemes for evaluating the layer potentials at a particular point for different cases have been extensively studied. When is far away from the boundary, classical Newton-Cotes or Gaussian quadratures for a general smooth integrand can be applied; when is located on the boundary, special quadrature rules can be designed, for example, the trapezoidal rule with end-point corrections in [1, 2, 22, 25, 28] or the generalized Gauss quadrature rules in [6, 7, 29]; and when

is close to the boundary, existing techniques include the change of variables to remove the principal singularity and the regularized kernel and corrections using asymptotic analysis

[4, 9, 14, 19]. We particularly mention the pioneering work in [20]

which applies the Barycentric Lagrange polynomial interpolation formula to derive a globally compensated spectrally accurate quadrature rule for evaluation at points close to the boundary (also see

[3]); and the pioneering “quadrature by expansion” (QBX) scheme in [23] which derives a partial wave (harmonics) expansion valid in a region close to (or even containing points on) the boundary. The QBX scheme has been combined with the fast multipole method (FMM) in [27] for solving the integral equation reformulation of PDEs.

In this paper, we introduce new representations of the Laplace layer potentials that are valid in the entire leaf (childless) box in the FMM hierarchical tree structure. As the representation can be evaluated at any point in the box, the numerical scheme also belongs to the class of “quadrature by expansion” (QBX) schemes to evaluate the layer potential integrals. In Fig. 1, the leaf boxes in a uniform FMM tree with levels are categorized into three groups: The green box is well separated from the boundary source points, from established fast multipole method (FMM) theory [16, 17], the layer potential at each target point in the box can be represented by a complex Taylor polynomial expansion which is referred to as the local expansion of the green box in the FMM algorithm. Both the red and yellow boxes are not well separated from the boundary, and each of the red boxes contains target points located both inside and outside the boundary, hence two separate solution representations become necessary, one for the interior and one for the exterior. For the red and yellow boxes, contributions from the well-separated curved line segments of the layer potential can still be represented using the complex local Taylor polynomial expansion, which can be efficiently computed using the FMM through the upward and downward passes for the “far-field” contributions of the sources. The numerical difficulty is an accurate and efficient representation of the near-field source contributions.

Using both the complex local Taylor polynomials and plane wave (exponential) functions in the basis, we propose a new representation of the 2 layer potentials for the red and yellow boxes due to the near-field layer potential source contributions. Combining both the far-field and near-field (local) density contributions, the main contribution of this paper is that the Laplace layer potentials inside each 2 leaf box of the FMM hierarchical tree structure can be represented as the sum of two expansions


Here, is a target point in the leaf box centered at , the boundary is described by , is a point close to determined by solving the equation when is inside either the red or yellow boxes (assuming after proper translations and rotations), the operator takes the real part of a complex number, and are the complex coefficients of the polynomial and plane wave expansions, respectively, the complex number is referred to as the node for the exponential (plane wave) expansion, and and are respectively the numbers of terms in the local Taylor polynomial and plane wave expansions ( for the green boxes).

The local complex polynomial expansion only collects the far-field contributions and its number of expansion terms can be analyzed using the established error analysis from the classical fast multipole method. The plane wave type expansion is derived by applying the Fourier extension technique to the density and boundary geometry functions followed by analytically utilizing the Residue Theorem for complex contour integrals, and the number of terms is the same as the number of terms required in the Fourier extensions.

Figure 1: Different expansions for the leaf boxes in a uniform FMM hierarchical tree structure. Green: complex polynomial expansion; Yellow: one QB2X for the leaf node; Red: two QB2X required, one for the interior and one for the exterior.

Note that two different types of basis functions are used in the representation of the layer potential, hence we refer to our approach as the quadrature by two expansions (QB2X). In harmonic analysis, the redundant basis functions form a frame [10, 13]. Compared with classical QBX, the QB2X representation of the layer potential is valid in a much larger region and allows easier analysis of the error and its dependency on the numbers of expansion terms. Another nice feature of the new representation is that the nonlinear impact of the boundary on the layer potential becomes explicit in Eq. (2), providing an analytical tool useful for other applications.

We organize this paper as follows. In Sec. 2, we review the classical QBX and the well-established Fourier extension technique which form the foundation of the QB2X technique. In Sec. 3, we derive the new representations for the single and double layer potentials in the red and yellow boxes using both the complex Taylor polynomial and plane wave basis functions. In Sec. 4, we present preliminary numerical experiments to demonstrate the accuracy of the new representations in the leaf box and to validate our analytical results. Finally in Sec. 5, we summarize our results and discuss our current work to generalize QB2X to layer potentials for other types of equations in both two and three dimensions.

2 Preliminaries

The new quadrature by two expansions (QB2X) technique uses two different basis functions, the complex polynomial expansion (also referred to as the local expansion in the fast multipole method) and the plane wave expansion using exponential functions. In this section, we present (a) the original QBX [23] which introduces the complex polynomial expansion (or partial waves for the Helmholtz and Yukawa equations) to evaluate layer potentials and (b) the Fourier extension technique [5, 8, 21] which will provide explicit formulas for the plane wave expansion in the new QB2X technique.

2.1 Quadrature by Expansion: Evaluating Layer Potential Using Complex Polynomial Expansion

Assuming both the density function and boundary curve are sufficiently smooth, to evaluate the singular, near-singular, or hyper-singular layer potentials, in [23], it was observed that the layer potentials are smooth functions on either side of the boundary, and the integrand singularity is only associated with the non-smoothness across the boundary. Therefore the polynomial expansion of the Laplace layer potential centered at a point either in the interior or exterior of the boundary is valid at least locally. In Fig. 2, we consider a Laplace layer potential explicitly given by in a box where is the complex variable. The box contains part of the boundary given by the equation . We assume the polynomial expansion is in the form centered at . We neglect the numerical errors when evaluating the expansion coefficients, i.e., the coefficients are derived exactly. In (a), we plot the analytical layer potential. In (b) and (c), we present the polynomial approximation of the layer potential using and terms in the expansion, and in (d), (e), and (f), we plot the errors of the representation when , , and , respectively. Clearly, when the number of expansion terms increases, the error decreases and the representation becomes valid in a much larger region.

(a) Analytical layer potential
(b) terms expansion
(c) terms expansion
(d) terms error
(e) terms error
(f) terms error
Figure 2: An implementation of QBX. (a): Analytical layer potential. (b) and (c): approximation using and terms, respectively. (d), (e), and (f): errors for , , and .

Using standard error analysis from the FMM theory, assuming the coefficients are computed accurately for the well-separated green boxes in Fig. 

1, the local expansion in the form can achieve -digits accuracy when is approximately and -digits accuracy when . The task of evaluating the layer potentials at points in the red and yellow boxes becomes more complicated. In [23], the choice of the local complex polynomial expansion center, the degree of the polynomial, and the quadrature schemes for computing the expansion coefficients are numerically studied to provide guidelines on these parameter selections. In [15]

, estimates for the rate of convergence of these local expansions are derived, which can be used to analyze the approximation error of the local Taylor expansions in a leaf box. In existing FMM+QBX implementation

[27], for evaluation points in one leaf node, several complex polynomial expansions may have to be formed with different expansion centers and degrees, unless the FMM hierarchical tree oversamples the density on the boundary or the solution. We also mention that in existing QBX implementations, the curvature of the boundary does not explicitly appear in the formulas. This will be addressed in the new QB2X in Sec. 3.

2.2 Fourier Extension: Approximation Using Exponentials

Compared with a polynomial basis, the exponential expansions, if it can be derived accurately and efficiently, may show better numerical properties in efficiency. One example is the translations in the FMM algorithms. When the exponential (plane wave) expansions are used, the translations become diagonal and the number of operations is reduced from the polynomial expansion’s to the plane wave expansion’s when terms are used in both expansions [11, 12, 18]. Unfortunately, deriving the optimal exponential expansion for a general function requires nonlinear optimization, and the uniqueness of the solution is not guaranteed. However, in some particular cases, a good exponential expansion approximation can be derived. For example, when the function is smooth and periodic, then the Fourier series can be computed efficiently and the expansion converges rapidly. Another example is when the inverse integral transform of the function is available, then the exponential expansion problem becomes an integration problem, and the weights and nodes of the exponential expansions can be computed using the generalized Gauss quadrature method [24, 29].

Figure 3: (a) and (c): Fourier extension (solid line) of the given function (dashed line), . (b) and (d): approximation errors on .

Using the Fourier series to approximate a non-periodic function is also a well-studied topic. Consider a non-periodic function on , applying the Fourier extension technique, a suitable periodic function on a larger domain is computed stably, so the Fourier series expansion of matches on the interval , see [5, 8, 21] and references therein. We have implemented the scheme in [21], which solves the least square optimization problem to compute an accurate Fourier series representation of a smooth function defined on . In Fig. 3, we present the computed Fourier series with fundamental period for the two Chebyshev basis polynomials and . The Fourier series approximation on achieves machine precision accuracy for both cases. Using the Fourier extension technique, a translation matrix can be precomputed to map the commonly used Chebyshev or other orthogonal polynomial basis functions to the Fourier basis functions that are periodic in a larger domain.

3 Quadrature by Two Expansions: Combining Complex Polynomial and Plane Wave Expansions

3.1 Problem Setup

As all the far-field layer potential density contributions can be accurately and efficiently computed using the fast multipole method, in this section, we focus on the near-field (local) density contributions to the red and yellow boxes in Fig. 1. We assume the boundary is parametrically described by , and after proper scaling, translation, and rotation as in Fig. 4. One of the yellow leaf box is shown and we assume the two end points and are well-separated from this leaf box with center . Following previous work in [27], we assume both the density function and boundary are well resolved by polynomials for .

Figure 4: A leaf box close to the boundary.

Note that the degree of the polynomials to approximate and can be different. Next, we present the detailed representation formulas for the single and double layer potentials


where is the target point in the yellow or red leaf box, is the source on the boundary, and .

3.2 Laplace Double Layer Potential

We start from the double layer potential to avoid the branch cut analysis of the function in the kernel. To evaluate the double layer potential at the target point in the yellow or red leaf box, the contribution from the source density function defined on the boundary segment , becomes


where , . The curvature of the boundary at is given by . In order to apply the complex contour integral theory and Residue Theorem, using

we get

The double layer potential becomes


As and are resolved by polynomials, therefore both terms are in the form of the complex integral


where is a polynomial defined for .


We first consider the case when to simplify the discussions and formulas, i.e., is on the straight line segment connecting and . In this case, Eq. (6) becomes

As is a function defined on the real line segment, the Fourier extension technique can be applied using the precomputed translation operator from the polynomial basis to the Fourier series, so

and Eq. (6) becomes


where only contains the non-negative frequencies and contains the negative ones.

Figure 5: The upper (left) and lower (right) contours for and , respectively.

We first study as part of the contour integral

where the contour is shown in the left plot of Fig. 5. It consists of the line segment from to and the semi-circle (denoted by ) on the upper half plane. As the integrand is analytic inside the contour and the pole is located at outside the contour, by the Residue Theorem, we have



Unlike the line segment from to , the semi-circle is well-separated from the leaf box containing and the contribution from the density defined on the semi-circle can be collected into a local expansion as in


where is the center of the leaf box, are the local expansion coefficients and the number of terms is controlled by the decay rate of which can be easily estimated for the given contour and target leaf box following the standard fast multipole method analysis. Introducing and noting that when is non-negative, decreases when on the contour moves away from the line segment, a very loose estimate of the truncation error is given by for some constant . When , it is straightforward to verify that , , , and will provide results with at least , , , and accurate digits, respectively. Also, the numerical stability issues associated with exponentially growing values can be avoided and the local expansion coefficients can be computed accurately using standard quadrature rules. Therefore, the non-negative modes can be represented as a local complex Taylor polynomial expansion. Unfortunately, the negative frequencies cannot be computed using this contour, as the function grows exponentially when on the contour moves away from the real axis for .

To compute , a different contour on the lower half complex plane has to be chosen, a sample contour is show on the right plot of Fig. 5. It consists of a rectangle with one side being the line segment from to , and a sufficiently large (can be ) constant is introduced to determine the length of the other side. We denote the part of the contour consisting of the three other sides of the rectangle by . Note that decays exponentially when on moves away from the real axis for . As is inside the contour, applying the Residue Theorem, for each negative frequency , we have

Therefore, we can compute using


where the local expansion (second summation in the formula) coefficients are given by

which are derived using the same separation of variables as in . As the leaf box is well-separated from , the number of the local polynomial expansion can be determined using the same FMM error analysis as in and we skip the details.

Combining , , and the far-field density contributions for this special case, we conclude that the double layer potential in the leaf box can be represented as a combination of the local complex Taylor polynomial expansion and plane wave expansion as in Eq. (2). The number of terms in the local polynomial expansion is determined by standard FMM error analysis as all the involved contributions are well-separated from the leaf box. The number of terms in the plane wave expansion is the same as that in the Fourier extension of the density and boundary functions for .

Comment on : Smaller values are possible by including a larger portion of the boundary when computing and , e.g., by also including contributions from the second nearest neighbors [26], at the cost of more terms in the Fourier extension. The balance of the numbers of terms in the exponential expansion and in the polynomial expansion is related with the optimal discretization strategies when generating the FMM adaptive tree, which is currently being studied.

Comment on the contours: We mention that the choice of the contour is not unique. Other contours can also be used. However changing the contours will not change the values of the expansion coefficients. It will only change the accuracy and efficiency of the numerical integration scheme for computing these values and the estimated number for truncating the local Taylor expansion. Our numerical experiments in Sec. 4 show that the current choice allows accurate and efficient computations of these coefficients and provides acceptable bounds for . Finding the “optimal” contours is still a challenging task and is being studied.

3.2.2 Curved Line: Role of Boundary Geometry

Next we consider the role of the boundary geometry in the representation. We assume is a polynomial of , at , and . The meaning of the last assumption is that the boundary is sufficiently resolved by , so that the boundary curve in each leaf box is reasonably close to flat. In this case, the double layer potential representation will depend on the boundary geometry described by the polynomial . When the degree of the polynomial is less than , the roots of can be derived analytically, and the one close to will be denoted as . The assumption that for small also implies that all other roots will be very far away from the leaf box and the region enclosed by the contours. Similar to the case, we replace by its Fourier series expansion, and separate the integral in Eq. (6) into two parts


To evaluate , we use the same contour on the left of Fig. 5. Note that the semi-circle (denoted by ) is well-separated from the leaf box, the integrand is analytic in the region enclosed by the contour, and decays exponentially when and . Therefore no numerical stability issues will appear and can be approximated using the same strategy as in the case by a complex polynomial expansion


where is the center of the leaf box,

are the local expansion coefficients, and the number of terms can be determined by standard FMM error analysis using the ratio for on the contour and in the leaf target box. Note that when , the ratio only changes slightly when compared with the case.

We also apply the same contour on the right of Fig. 5 to evaluate . As the root is inside the contour so we factor . Simple algebra will show that and the Residue Theorem becomes

We can therefore represent using the following formula,


where the local expansion coefficients are given by

and the number of terms can be estimated using the ratio and standard FMM error analysis. This formula shows that the representation depends nonlinearly on the geometry in three different ways, the additional term in the denominator of the (slightly modified) plane wave expansion, the nonlinear dependency when finding the root using the polynomial equation , and the ratio which depends on the function . When the degree of the polynomial is less than , an analytical formula is available to express explicitly as a function of . When the degree is higher, an asymptotic expansion can be derived to approximate using when assuming . The resulting mathematical formulas in the QB2X therefore reveal the solution dependency on the geometry and can be useful tools in PDE analysis.

Comment on : It is worth mentioning that even without this assumption, most parts in the analysis are still valid. It is therefore possible to apply the combined (slightly modified) plane wave and local complex Taylor polynomial expansions for much larger leaf boxes in the numerical discretization. However there are several numerical challenges: It becomes possible to have new poles (roots of ) moving inside the contour for or ; more polynomial expansion terms become necessary for a prescribed accuracy requirement; computing the local expansion coefficients may require new contours to avoid any numerical difficulties; and the discretization scheme to generate the FMM hierarchical tree structure also requires further study for optimal performance of the algorithm. These new challenges are being studied.

3.3 Laplace Single Layer Potential: Integration by Parts

To derive the representation for the single layer potential


where , , and We first apply the Fourier extension technique to represent the real function as and define

which is a particular anti-derivative of as . Using and integration by part, we have


For the term, as both end points are well-separated from the leaf box containing , a local expansion can be derived. For , simple algebra shows that


Both terms are in the form of Eq. (6) with different density functions, therefore results from previous section for double layer potentials can be applied directly. We skip the details.

3.4 Comparing Plane Wave Expansion with Complex Polynomial Expansion

The new QB2X technique uses both the local complex Taylor polynomial expansion and (the slightly modified) plane wave expansion. Different error sources in the QB2X can be easily analyzed: when the density function for is approximated by an orthogonal polynomial, the truncation error is a well-studied topic; the error from the Fourier approximation is controlled by the Fourier extension precomputation process which generates the linear mapping from the polynomial basis to Fourier basis; the error from the truncated local complex polynomial approximation of the far-field density contribution on the well-separated boundary segments or on the contour follows the standard FMM analysis when a proper contour is chosen; the (slightly modified) plane wave in the representation is derived analytically using the Residue Theorem; finally finding the roots of a polynomial is a well studied topic.

The QB2X also provides an alternative approach to analyze the error in the classical QBX method (see [15] for existing results), by studying the truncation error when a plane wave term is re-expanded as a local Taylor polynomial expansion


where is the expansion center. We assume , and study how many terms are required in the polynomial expansion in order to achieve machine precision accuracy, i.e., we need to find such that Using the incomplete gamma function and gamma function , we have

In Fig. 6, we numerically solve the inequality to get an estimate of . When is used in the plane wave expansion, to achieve machine precision the estimated is about , and this number becomes larger if the nonlinearity due to the boundary geometry is included. Therefore, introducing two different bases in the representation improves the efficiency when approximating the layer potentials as we may gain accuracy from the plane wave approximation that would take many terms from a complex Taylor polynomial expansion alone.

Figure 6: Estimated (y-axis) for different wave number (x-axis).

4 Numerical Experiments

We present preliminary numerical results to validate the analytical formulas and demonstrate the achieved accuracy for different (for the polynomial expansion) and (for the plane wave expansion) values.

4.1 Double Layer Potential

We first consider the straight line segment connecting and when . The double layer potential of interest is then given by


In existing implementations which combine QBX with FMM, the density function is often approximated by an orthogonal polynomial expansion, e.g., the Chebyshev polynomial expansion where is the Chebyshev basis polynomial given by The first three basis polynomials in terms of are explicitly given by , , and . In the first numerical test, we choose different density functions: (a) ; (b) ; (c) ; and (d) . For (a), as it is already in the form of an exponential expansion, so . We use in the Fourier expansion for (b) to guarantee machine precision accuracy. For (c) and (d), we apply the precomputed mapping from the polynomial basis to the Fourier basis to compute the Fourier extensions with . The approximation errors are also around machine precision. We assume the target point and the center of the leaf box is given by . For the and terms given explicitly in Eq. (7), we choose the contours in Fig. 5, where is used. The values are for the upper contour and for the lower contour, respectively. We choose which guarantees at least -digits accuracy in the complex local Taylor polynomial expansion using standard FMM error analysis. The approximation errors are shown in the Fig. 7. For all cases, the QB2X representations achieve -digits accuracy.

(a) ,
(b) ,
(c) ,
(d) ,
Figure 7: Approximation errors of QB2X representations for double layer potentials with for different density functions. The plot legends are .

In classical FMM error analysis, when , the complex local Taylor polynomial expansion is guaranteed to achieve -digits accuracy, and the accuracies increase to , , and digits when , , and , respectively. The same error estimates can be derived using the values in this example. In Fig. 8, we show the approximation error for different number of expansion terms for case (c) when