Zeta Correction: A New Approach to Constructing Corrected Trapezoidal Quadrature Rules for Singular Integral Operators

by   Bowei Wu, et al.

We introduce a new quadrature method for the discretization of boundary integral equations (BIEs) on closed smooth contours in the plane. This quadrature can be viewed as a hybrid of the spectral quadrature of Kress (1991) and the locally corrected trapezoidal quadrature of Kapur and Rokhlin (1997). The new technique combines the strengths of both methods, and attains high-order convergence, numerical stability, ease of implementation, and compatibility with the “fast” algorithms (such as the Fast Multipole Method or Fast Direct Solvers). Important connections between the punctured trapezoidal rule and the Riemann zeta function are introduced, which enable a complete convergence analysis and lead to remarkably simple procedures for constructing the quadrature corrections. The paper reports a detailed comparison between the new method and the methods of Kress, of Kapur and Rokhlin, and of Alpert (1999).



There are no comments yet.


page 12


Corrected Trapezoidal Rules for Boundary Integral Equations in Three Dimensions

The manuscript describes a quadrature rule that is designed for the high...

High order corrected trapezoidal rules for a class of singular integrals

We present a family of high order trapezoidal rule-based quadratures for...

Corrected trapezoidal rules for singular implicit boundary integrals

We present new higher-order quadratures for a family of boundary integra...

High-order Corrected Trapezoidal Rules for Functions with Fractional Singularities

In this paper, we introduce and analyze a high-order quadrature rule for...

FMM-LU: A fast direct solver for multiscale boundary integral equations in three dimensions

We present a fast direct solver for boundary integral equations on compl...

Fourier smoothed pre-corrected trapezoidal rule for solution of Lippmann-Schwinger integral equation

For the numerical solution of the Lippmann-Schwinger equation, while the...
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

This paper describes techniques for discretizing boundary integral equations (BIEs) of the form


where is a smooth closed contour in the plane and the arc length measure on , where is a given smooth function, and where is a given kernel with a logarithmic singularity as . Equations such as (1) commonly arise as reformulations of boundary value problems from potential theory, acoustic and electromagnetic wave propagation, fluid dynamics and many other standard problems in engineering and science. When a PDE can be reformulated as an integral equation that is defined on the boundary of the domain, there are several advantages to doing so, in particular when the BIE is a Fredholm equation of the second kind.

A key challenge to using (1) for numerical work is that upon discretization, it leads to a system of linear equations with a dense coefficient matrix. Unless the problem is relatively small, it then becomes essential to deploy fast algorithms such as the Fast Multipole Method (FMM) [7] or Fast Direct Solvers (FDS) [15]. A second challenge is that the singularity in the kernel function

means that if a standard quadrature rule is used when discretizing the integral, then only very slow convergence is attained as the number of degrees of freedom is increased.

This paper introduces a new family of quadrature rules for discretizing (1) that are numerically stable even at high orders. For instance, a rule of order 42 is included in the numerical experiments. It is perfectly stable and capable of computing solutions to 14 correct digits with as few degrees of freedom as spectrally convergent quadratures such as the method of Kress [12]. Moreover, unlike the Kress quadrature, it can easily be used in conjunction with fast solvers such as the FMM or FDS.

1.1 Nyström discretization and corrected trapezoidal rules

Upon parameterizing the domain over an interval , a BIE such as (1) can be viewed as a one dimensional integral equation of the form


In (2), the new kernel encodes both the parameterization and the original kernel . Observe that all functions in (2) are -periodic, that and are smooth, and that is smooth except for a logarithmic singularity as .

To discretize (2) using the Nyström method [13, §12.2], we consider the -point periodic trapezoidal rule (PTR)


where and . When is smooth and periodic, the PTR converges super-algebraically as [19]. The Nyström method first collocates (2) to the quadrature nodes of the PTR, and then approximates the integral by a quadrature supported on the same nodes with unknowns , yielding a linear system


The coefficient matrix should be formed so that the approximation


holds to high accuracy. If were to be smooth, this task would be easy, since we could then use the PTR (3) without modifications and simply set


In this unusual case, the solution of (4) will converge super-algebraically to as .

In the more typical case where is logarithmically singular at , some additional work is required to attain high order convergence in (5). Let us start by describing two existing methods that resolve this problem — the Kress quadrature and the Kapur-Rokhlin quadrature — that are closely related to the new quadrature that we will describe.

Kress [12] introduced a quadrature that is spectrally accurate for any periodic function of the form


where and are smooth functions with known formulae. Kress integrates the first term of (7) by Fourier analysis and the second term using the PTR, resulting in a corrected trapezoidal quadrature where all the PTR weights are modified. It is not obvious how a BIE scheme based on the Kress quadrature can be accelerated by the existing fast algorithms. We will use a “localized” analog of the analytic split (7) to develop our quadrature.

Kapur and Rokhlin [10], on the other hand, constructed a family of quadratures for a variety of singular functions by correcting the trapezoidal rule locally near the singularity. These quadratures have correction weights that are essentially independent of the grid spacing and possess an essential benefit in that nearly all entries of the coefficient matrix are given by the simple formula (6); only a small number of entries near the diagonal are modified. This local nature of Kapur-Rokhlin quadrature makes it very easy to combine it with the FMM and other fast algorithms. For the logarithmic singularity, two different quadratures are developed. The first quadrature is for functions of the “nonseparable” form


where the formulae of the smooth functions and can be unknown. This “nonseparable” quadrature ignores the data at completely and modifies a few trapezoidal weights on both sides of the singular point. The magnitude of the correction weights (tabulated in [10, Table 6]) grow rapidly with the order of the correction, and moreover, some of the weights are negative. These properties mean that Kapur Rokhlin quadrature becomes less useful at higher orders (say order higher than 6), since the resulting coefficient matrix can be far worse conditioned than the underlying BIE [8, Sec. 7.3]. The second Kapur-Rokhlin quadrature is for functions of the “separable” form


Unlike the first rule, this “separable” quadrature also uses the data at the singular point for its correction; the correction weights (tabulated in [10, Table 7]) are uniformly bounded regardless of the correction order, and decay rapidly away from the singular point. Despite the excellent stability properties of the second Kapur-Rokhlin rule, it has received little attention due to the simple fact that the kernels arising from BIEs typically are of the nonseparable type (8). (In fact, the first rule is often referred to as “the Kapur-Rokhlin” rule.)

To avoid the issue of large correction weights, Alpert [3] developed a hybrid Gauss-trapezoidal quadrature that uses an optimized set of correction points that are off the uniform trapezoidal grid, whose weights are very well-behaved. For additional details on high order accurate techniques for discretizing (2), as well as a discussion of their relative advantages, we refer to the survey [8].

1.2 Contributions

This paper describes a quadrature rule that is closely related to the neglected second Kapur-Rokhlin rule for separable functions. The new rule works almost exactly the same in practice in that it involves a small correction stencil that includes a correction weight at the origin. Both rules display excellent numerical stability and lead to discretized systems that are as well conditioned as the original equation. The key innovation is that the new rule that we present is applicable to functions of the form (8) (provided that formulae for and are known), which makes these rules applicable to a wide range of BIEs.

The new quadrature is derived from a local kernel split analogous to Kress’ analytic split (7). We analyze the error of applying the punctured trapezoidal rule to the singular component of this kernel split based on the lattice sum theory. This results in an error expansion with coefficients that are explicitly computable using the Riemann zeta function, a fact we called the “zeta connections.” From these error coefficients, we construct local correction weights in the fashion of Kapur and Rokhlin. As it turns out, the correction weights constructed this way have a component that is the weights of the “separable” Kapur-Rokhlin quadrature mentioned above, while the remaining component depends on the explicit kernel split. (Remarkably, the zeta connection simplifies the construction of the separable Kapur-Rokhlin correction weights to the extent that Table 7 of [10] can be computed with three lines of code, as shown in Figure 1.)

The zeta connection associated with the singular function , was first introduced by Marin et. al. [14]. In this paper, we extend this connection to a “differential zeta connection” (Theorem 3) associated with the logarithmic singularity. We would also like to point out that the zeta connection has recently been generalized to higher dimensions by the authors in [20]. Combining the “differential zeta connection” in this paper with the higher dimensional zeta connection of [20], we expect that a rigorous theory can be developed for higher-dimensional logarithmic quadratures such as the one developed by Aguilar and Chen [1]. On the other hand, there is also the connection between the zeta function and the endpoint corrections of the trapezoidal rule for regular functions, which has been established much earlier, see [16, 3].

1.3 Organization

In Section 2, we introduce the theory for the local correction of the trapezoidal quadrature and its connection to the zeta function. We extend this connection to construct a quadrature for the logarithmic singularity, which recovers the “separable” Kapur-Rokhlin quadrature but with much simpler computations. In Section 3, we generalize this Kapur-Rokhlin rule to logarithmic kernels on planar curves using a localized version of Kress’ kernel split, developing quadratures for the Laplace and Helmholtz layer potentials. Finally in Section 4, we present numerical examples of solving BIEs associated with the Helmholtz and Stokes equations and compare our quadrature method with existing state-of-the-art methods.

2 Corrected trapezoidal rules and the zeta connections

In this section, we introduce the theory for the local correction of the trapezoidal rule for singular functions. In particular, we introduce the “zeta connections,” based on which simple and powerful procedures are presented (see Remark 1) for constructing quadratures for functions with algebraic or logarithmic branch-point singularities.

2.1 Singularity correction by moment fitting

To set up the notation for our discussion, we let denote the integral of a function on the interval , where , and where may be singular at . We denote the punctured trapezoidal rule discretization of as


for some integer , where and , and where the prime on the summation sign indicates that is omitted. (Note that by our definition of , the endpoints are not included as quadrature nodes, thus the usual factor for the weights at the endpoints is not needed. However, this choice is just for convenience and using the usual trapezoidal rule does not affect our analysis.)

We are interested in singular integrals of the form


where has an isolated, integrable singularity at and is smooth. We assume that either the integrand is periodic or that is compactly supported in so that the only obstruction to high-order convergence is the singularity of at . In general, boundary corrections to the trapezoidal rule can be introduced near independent of the singularity at ; see [2], for example, for more detail.

To analyze the error of the approximation , we decompose the integrand into a regular component and a local component:


where is a smooth and compactly supported function that is at least -time continuously differentiable and which satisfies and . For the regular component , the trapezoidal discretization


holds for any (note that the integrand is at ). Thus the overall convergence of is restricted by the error in the local component .

Using the idea of moment fitting

[11, 10], this local error can be corrected by fitting a set of moment equations for monomials up to a sufficiently high degree as follows


where is an integer and the factor depends explicitly on the singularity at , and where the double prime on the summation sign indicates that the term is multiplied by 2. (For convenience we let .) There are equations in (14) for the unknowns . Then combining (1314) yields a locally corrected trapezoidal quadrature


which is high-order accurate as .

2.2 The singularity and converged correction weights

The fact that the weights in (15) depend on is inconvenient in practice. Fortunately, this can be remedied by what we called the “zeta connection.” When , , Marin et. al. [14] showed that by letting only in the term in the moment equations (14), the right-hand side of (14) in this limit can be represented as Riemann zeta function values, and the corresponding limiting weights are independent of ; more importantly, they proved that using these converged weights in the quadrature (13) in place of does not affect the order of accuracy. We summarize this result of [14] in this section, which will serve as the foundation for our extensions to other quadratures.

We first introduce the important concept of “converged correction weights.”

Substitute in (14) and let , we have


where the last equality holds because is compactly supported. Note that both and are even, so by requiring that

both sides of (16) vanish for all odd. Using this symmetry and eliminating on both sides yield


From here we define the converged correction weights to be


such that are the solution of (17). To further simplify the equations, we need the following theorem.

Theorem 1 (Zeta connection).

For all ,


Consequently, the converged weights , as defined by (18), are the solution of the system


The zeta connection (19) is proved in [14, Lemma A2]. Then taking in (17) yields (20). (Note that we used the fact that .) ∎

Based on the zeta connection, the next theorem constructs a high-order corrected trapezoidal rule using converged correction weights.

Theorem 2.

For , one has the locally corrected trapezoidal rule


where the correction weights are the solution of (20), and where , must have at least vanishing derivatives at , i.e.


See Theorem 3.7 and Lemma 3.8 of [14]. ∎

2.3 Logarithmic singularity and the differential zeta connection

One can also construct a quadrature for the logarithmic singularity by extending the zeta connection (19) based on the following simple observation:

This leads to the next two theorems that are completely analogous to Theorems 1 and 2.

Theorem 3 (Differential zeta connection).

For all ,


Consequently, if we define the converged weights such that are the solution of the system


then are the solution of the system


Taking the derivative with respect to under the limit sign on both sides of (19) yields (23), which is justified since the expression under the limit sign is analytic in . Then taking in (24) yields (25). ∎

Theorem 4.

For , one has the locally corrected trapezoidal rule


where the correction weights are the solution of (25), and where satisfies the same conditions as in Theorem 2.


The local error of the punctured trapezoidal rule for is


Notice that in the second parentheses, the integral is smooth so the regular trapezoidal rule converges super-algebraically, i.e. using the fact that we have


holds for any . On the other hand, using the idea of moment fitting and following a similar derivation as (1617), the terms in the first parentheses of (27) are approximated by


where are the solution of (24). Then substituting (2829) into (27) gives

The above equation implies (26) once it is shown that

which in turn can be proved by showing that the limit (23) converges as uniformly for given the condition of . This last statement can be proved following almost verbatim the proofs of Theorem 3.1 and Lemma 3.3 in [14] by replacing with therein, hence we omit the detail here. ∎

The logarithmic quadrature (26) is equivalent to the “separable” Kapur-Rokhlin quadrature developed in [10, §4.5.1], hence the correction weights are identical (up to a minus sign) to those given in [10, Table 7]; however, the differential zeta connection has greatly simplified the construction of these weights.

Remark 1.

In practice, the value , , can be approximated using “complex step differentiation” [18] as


where , . This formula is free of cancellation errors that plagued typical finite difference methods. For instance, using will yield an approximation of full double-precision accuracy.

On the other hand, as mentioned in [10], the Vandermonde system (25) is ill-conditioned for large . Thus when precomputing the weights , (25) should be solved symbolically or under extended precision. Simple code snippets that compute for any given are given in Figure 1.

MATLAB: [frame=single] rhs = -imag(zeta(-vpa(0:2:2*K)’+1i*eps)/eps); V = vpa(0:K).^((0:2:2*K)’); w = double(V); Julia: [frame=single] using SpecialFunctions rhs = -imag.(zeta.(-(0:2:2*K).+im*eps())./eps()) V = (BigFloat.(0:K))’.^(0:2:2*K) w = V Mathematica: [frame=single] Unprotect[Power]; Power[0,0]=1; Protect[Power]; (* set 0^0=1 *) rhs = -Im[Zeta[-2(Range[K+1]-1)+I*MachineEpsilon]; V = Array[(#2-1)^(2(#1-1))&, K+1,K+1]; w = Inverse[V].rhs;

Figure 1: Code that given (such that the correction order is ) constructs the correction weights for integration against , as described in Remark 1 (where is chosen to be the machine precision eps for simplicity). The code generates [10, Table 7] up to a minus sign.

3 Logarithmic kernels on curves

In this section, we extend the “separable” Kapur-Rokhlin rule (26) to construct our “zeta-corrected quadrature.” We will combine the differential zeta connection with local kernel splits (that are analogous to Kress’ global analytic split (7)) to construct quadratures for some important logarithmic kernels on closed curves, including the Laplace and Helmholtz layer potentials; the quadrature for the Laplace single-layer potential will also be applied to integrate the Stokes potential in Section 4.

3.1 Laplace kernel

Consider a smooth closed curve parameterized by a -periodic function . We consider the Laplace single-layer potential (SLP) from to (for the general case, one simply replace with any other target point on ),


where and . The next theorem extends Theorem 4 to construct a corrected quadrature for (31).

Theorem 5.

For the Laplace SLP (31), one has the locally corrected trapezoidal rule


where and is smooth, and where the correction weights are exactly the same as in Theorem 4.

Note that the only difference of (32) from (26) is that is replaced with and is replaced with .


First we analyze the singularity of . Note that


which can be rewritten as


We will show that (34) is smooth, thus the only singular term in (33) is . To this end, first expand as a Taylor-Maclaurin series



is smooth near , which implies that (34) is indeed smooth. Next, we use the decomposition

to analyze the error of the punctured trapezoidal rule being applied to , as follows


where, because (34) is smooth, the terms in the first curly brackets of (35) happen to be the error of the regular trapezoidal rule applied to a smooth function (notice that the integrand is zero at ), which vanishes super-algebraically; the terms in the second curly brackets, analogous to (28), converge to super-algebraically; finally for the terms in the last curly brackets, one simply applies the quadrature (26) of Theorem 4, with replaced by

. Combining all these estimates, as well as (

13), one concludes that (35) implies (32). ∎

3.2 Helmholtz kernels

We now apply Theorem 5 to construct formulae for the Helmholtz layer potentials. Consider the Helmholtz SLP on a smooth closed curve evaluated at ,


where and is the wavenumber, and where the kernel has the form [5, §3.5]


where and are, respectively, the Hankel and Bessel functions of the first kind of order , where is some smooth function of such that , and where such that is Euler’s constant. Analogous to Kress’ analytic split (7), we introduce the kernel split

where , such that the component is smooth. Therefore we can split (36) as


where and


are smooth function. Notice that the singular integral can be approximated by (32) with replaced by , giving


where we have used the fact that . Then combining (40) with the PTR for (which converges super-algebraically), we finally have


where, again, is given by (39).

Finally, we can also obtain the formulae for the Helmholtz double-layer potential (DLP), , and the normal derivative of the SLP, , using similar derivations. We will just state the formulae and omit the derivations. Using similar notations from (36) and , these layer potentials are given by


where with being the unit outward normal at , and where with and . The corresponding corrected trapezoidal rules for and are, respectively,


where the curvature at scaled by , and where

with being the Bessel function of the first kind of order .

4 Numerical experiments

In this section, we present numerical examples of solving BIEs associated with the Stokes and Helmholtz equations. In each case, we obtain a linear system of the form (4), where the matrix is filled using a particular quadrature. Then the linear system is solved either directly by inverting the matrix or iteratively by GMRES.

We compare our quadrature method with the three singular quadratures mentioned in the introduction: Kapur and Rokhlin’s locally corrected trapezoidal quadrature [10], Alpert’s hybrid Gauss-trapezoidal quadrature [3], and Kress’s spectral quadrature [5, §3.6]. The correction weights for our quadrature are precomputed by solving the equations (25) and using the techniques described in Remark 1. When implementing the Kapur-Rokhlin, Kress, and Alpert quadratures, we followed the survey [8].

(a) (b) (c)
Figure 2: Problem setup for the tests in Figures 3 and 4. the Stokes problem (45) and the Helmholtz problem (46). In all cases, the star-shaped geometry is parameterized by the polar function , while the diamonds represent testing locations. (a) Streamlines of a shear flow around an island with no-slip boundary condition. (b) Real part of a wave field generated at the source locations indicated by the dots. The wavenumber is . (c) Same as in (b), except that the wavenumber is now , so the wave decays exponentially

Stokes problem.

As shown in Figure 2(a), consider a viscous shear flow around an island whose boundary is a smooth closed curve , with no-slip boundary conditions on . Let be the true velocity field and its associated pressure field, then is described by the exterior Dirichlet problem for the Stokes equation [9, §2.3.2]

The integral equation formulation for this problem is obtained using the mixed potential representation for the velocity [17, §4.7], where the integral operators

are the Stokes SLP and DLP in 2D, where and is the unit outward normal to at , and where

denotes the tensor product. Then the vector-valued unknown density function

is the solution of the following BIE


As , the only singular component in the linear operator is the term in , which can be efficiently handled by the corrected trapezoidal rule (32).

Figure 3 compares the convergence results of solving the Stokes problem using different quadrature methods. We see that all the quadratures have the expected convergence rates, with the Kapur-Rokhlin quadrature yielding higher absolute errors due to its larger correction weights. The Kress quadrature is the most accurate at virtually any given , but the new quadrature of order 16 is remarkably close.

Figure 3: Comparison of singular quadratures for solving the Stokes problem (45) using the setup in Figure 2(a). The reference solution is obtained using the Kress quadrature with points on .

Helmholtz problem.

As shown in Figure 2(b) or (c), consider the Helmholtz Dirichlet problem exterior to the curve , with boundary data and satisfying the Sommerfeld radiation condition,