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)  or Fast Direct Solvers (FDS) . 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 . 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 .
where and . When is smooth and periodic, the PTR converges super-algebraically as . 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  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 , 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  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 .
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  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. . 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 . Combining the “differential zeta connection” in this paper with the higher dimensional zeta connection of , we expect that a rigorous theory can be developed for higher-dimensional logarithmic quadratures such as the one developed by Aguilar and Chen . 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].
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 , 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 (13–14) 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.  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  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
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
Based on the zeta connection, the next theorem constructs a high-order corrected trapezoidal rule using converged correction weights.
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 . ∎
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:
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
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
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  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.
In practice, the value , , can be approximated using “complex step differentiation”  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.
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 ),
First we analyze the singularity of . Note that
which can be rewritten as
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
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 , Alpert’s hybrid Gauss-trapezoidal quadrature , 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 .
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 whereis 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.
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,