An adaptive kernel-split quadrature method for parameter-dependent layer potentials

by   Ludvig af Klinteberg, et al.
Simon Fraser University

Panel-based, kernel-split quadrature is currently one of the most efficient methods available for accurate evaluation of singular and nearly singular layer potentials in two dimensions. However, it can fail completely for the layer potentials belonging to the modified Helmholtz, biharmonic and Stokes equations. These equations depend on a parameter, denoted α, and kernel-split quadrature loses its accuracy rapidly when this parameter grows beyond a certain threshold. The present report describes an algorithm that remedies this problem, using per-target adaptive sampling of the source geometry. The refinement is carried out through recursive bisection, with a carefully selected rule set. This maintains accuracy for a wide range of α, at an increased cost that scales as α. Using this algorithm allows kernel-split quadrature to be both accurate and efficient for a much wider range of problems than previously possible.



page 4


Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping

The method of Helsing and co-workers evaluates Laplace and related layer...

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

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

Highly accurate special quadrature methods for Stokesian particle suspensions in confined geometries

Boundary integral methods are highly suited for problems with complicate...

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

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

Neglecting discretization corrections in regularized singular or nearly singular integrals

A method for computing singular or nearly singular integrals on closed s...

Continuity estimates for Riesz potentials on polygonal boundaries

Riesz potentials are well known objects of study in the theory of singul...

Split representation of adaptively compressed polarizability operator

The polarizability operator plays a central role in density functional p...
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 report presents an extension to the panel-based, kernel-split quadrature scheme by Helsing and Ojala [9]

, for evaluating singular and nearly singular layer potentials in two dimensions. This scheme represents one of the current state of the art methods for maintaining low errors when solving elliptic partial differential equations (PDEs) in two dimensions using integral equation methods

[9, 7, 13]. However, there exists a set of problems for which this scheme can fail completely. This includes the following PDEs in :


where is real number. For brevity, we refer to them as the modified PDEs. Note that they are not consistently named in the literature. For example, the modified Helmholtz equation is also known as the screened Poisson equation, the Yukawa equation, the linearized Poisson-Boltzmann equation, and the Debye-Hückel equation. Meanwhile, the modified Stokes equations are also known as the Brinkman equations. These PDEs appear in many different applications: electrostatic interactions in protein and related biological functions, macroscopic electrostatics and fluid flow on the micro scale, to mention a few [15, 14, 10, 11, 5, 8]. They also appear as a result of applying implicit or semi-implicit time marching schemes to the heat equation and the time-dependent Stokes and Navier-Stokes equations [4, 6].

A common trait of the modified PDEs is that their associated layer potentials have kernels that either decay exponentially, or have components that decay exponentially, with a rate that is proportional to . This decay presents a problem for the abovementioned kernel-split quadrature. In short, the quadrature method is based on writing the kernel on a form with smooth functions multiplying explicit singularities, and then integrating each term separately. In order to be accurate, these smooth functions have to be well approximated by piecewise polynomials, which is no longer the case when is large.

We have developed a robust quadrature scheme, based on adaptive refinement, that maintains high accuracy for any

, without sacrificing effficiency. It applies for target points both on and close to the boundary, where regular quadrature is insufficient. In this context, refinement refers simply to an interpolation of known quantities to a locally refined discretization, as opposed to increasing the number of degrees of freedom in the discretized integral equation.

This report presents a work in progress, as the guidelines for setting parameters are entirely heuristic at the moment, even though they provide very satisfying result. We will at a later date update this report with analysis that supports these parameter choices.

The remainder of this report is organized as follows. In section 2 we give an outline of the kernel-split quadrature. Section 3 describes the problem that we are trying to solve, using the modified Helmholtz equation as an example. The new algorithm we propose is presented in section 4, followed by numerical results in section 5.

2 Background

Our goal is to evaluate layer potentials of the form


The layer density is assumed to be smooth, while the kernel is singular at . It can be expressed with explicit singularities as


where , and are smooth functions. We refer to this decomposition as kernel-split.

The boundary is discretized using a composite Gauss-Legendre quadrature. It is subdivided into intervals, denoted panels,


Each panel is described by a parametrization ,


Associated with the parametrization is a speed function

, a normal vector

and the curvature . Introducing the convenience notation , the layer potential from a panel becomes


Each panel is discretized in the parametrization variable using the nodes and weights of an -point Gauss-Legendre quadrature rule, which is of order , such that on each panel we have the discrete quantities


Omitting the panel index , the layer potential contribution from a panel is then computed using the approximation


Due to the singularities in , the above formula requires to be well-separated from (this notion will be clarified further). Otherwise, the scheme of [7] is used, known as product integration. With that, target-specific quadrature weights of order are computed for the known singularities where needed, such that


Substituting (5) into (4) and applying the above product integration gives the so called kernel-split quadrature scheme. By formulating it as a correction, we only need to explicitly evaluate at . Depending on the location of the target point relative to the source panel , the evaluation can be divided up into three different cases:

  1. Singular, with self-interaction. If is one of the quadrature nodes, , then the term multiplying is smooth, with the limit


    where is the curvature of at . Applying product integration to the term, we get

  2. Singular, without self-interaction. If is either on the source panel , or on a neighboring panel, then the term is still smooth, but we do not have to take the limit at into account. This lets us simplify the above to

  3. Nearly singular case. If is close to , but not on a neighboring panel, then we need to deal with both the singularities in (5). This case occurs when either is close to , or when is on a section of the boundary that is distant in arc length or disjoint. The layer potential is then evaluated as


Given a target point , the panels on are partitioned into two sets: far panels, that can be evaluated directly using (11), and near panels, that must be evaluated (or corrected), using either (15), (16), or (17).

3 Problem statement

The kernel-split quadrature scheme outlined above is both efficient and accurate when applied to the single and double layer potentials of several PDEs, such as the Laplace, Helmholtz and Stokes equations [9, 7, 13]. However, for the layer potentials of eqs. 3, 2 and 1, the scheme can fail completely. All of these equations have layer potential kernels including second-kind modified Bessel functions, in the forms and/or . As we shall see, this is problematic for the kernel-split quadrature.

Figure 1: Relative error , when evaluating the modified Helmholtz single layer potential (23) over a flat panel of length , using kernel-split quadrature. The maximum is taken over 100 values of randomly drawn from the box . The white lines are contours of the quantity , providing a heuristic motivation for why a criterion of the form (30) is suitable.

As an illustrating example, we study the single layer kernel of the modified Helmholtz equation (1),


The split of this kernel is based on using a standard decomposition [12, §10.31] to explicitly write out the singularities in ,


where is a modified Bessel function of the first kind, and is the smooth remainder. Inserting this into the kernel (18), after also splitting the logarithm, we identify the terms in (5) as


Even though goes to zero as , the functions and actually grow exponentially as , following the asymptotic as [12, §10.30]. As gets larger, this makes an increasingly bad candidate for polynomial interpolation, which is essential for product integration to be accurate. To illustrate this, we consider the single layer potential with unit density, evaluated using kernel-split quadrature from a flat panel of length ,


To evaluate this using the kernel-split correction (17), we first write


The remaining integral is evaluated by writing the integrand as a polynomial plus remainder,


and then computing that using product quadrature,




are integrals that can be computed recursively, starting from exact formulas. The error in the kernel-split quadrature is dominated by the integral in the right-hand side of (24), which integrates the remainder of the polynomial interpolation. To help quantify this remainder, we consider the case . From the power series expansion of , available at [12, §10.25], we then have that


such that, on the interval , the remainder is bounded by


This clearly grows monotonically with the quantity . Numerical results indicate that this is the case also for general values of near , see fig. 1. To ensure that the kernel-split quadrature error remains below some tolerance , we therefore suggest that and must satisfy a criterion on the form


for some constant , which can be determined using numerical experiments. Once determined, we find that this criterion holds well also for general geometries.

The above criterion can also be reformulated as follows: In order to achieve a tolerance , panel lengths must satisfy


The obvious way of achieving this, for a given , is to discretize using sufficiently short panels. However, this can result in a discretization with orders of magnitude more points than necessary to resolve the geometry and the layer density.

4 Quadrature by recursive subdivision

To circumvent the problem described above, we here introduce an algorithm for local refinement, based on panel subdivision. Given a single source panel , we assume that it is sufficiently short, relative to the quadrature order , for both the geometry and the layer density to be well-represented by a polynomial, interpolated at the quadrature points. We say that it is well-resolved. For a given target point , we can then subdivide into a set of subpanels , interpolate our known quantities from the quadrature nodes on to the quadrature nodes on those subpanels, and then evaluate the layer potential at using the subpanels. To ensure accuracy, this subdivision is formed in a way that guarantees that all subpanels either are short enough to satisfy (30), or are sufficiently far away from , relative to their own length, to not need kernel-split quadrature.

Before we can state our algorithm, a number of preliminaries are needed:

Preimage of target

Let be the mapping from the standard interval to the panel . Then , such that , is the preimage of the target point . The preimage is real-valued if , and complex-valued otherwise. We here assume that we know the value of , but need not be a known function; see [2] for a discussion on how construct a numerical representation.

Subpanels and subintervals

A subdivision of is defined by a set of edges in the parametrization, , such that a subpanel is given by the mapping of the subinterval under . We can, by a linear scaling, define the local mapping that maps the standard interval to as


where . Given the preimage , the local preimage , such that , is given by


Near evaluation criterion

Given the preimage of a point close to a panel (or subpanel)

, it is possible to compute an accurate estimate of the quadrature error when evaluating the layer potential using

-point Gauss-Legendre quadrature. Detailed discussions can be found in [1, 2]. To leading order, the error is proportional , where is the elliptical radius of the Bernstein ellipse on which lies,


where is defined as with . For a given kernel and error tolerance , it is then possible to introduce cutoff radius , such that kernel-split quadrature must be used for


and Gauss-Legendre quadrature is sufficiently accurate otherwise. Later we will use that the inverse of has a particularly simple form in the special case when lies on the imaginary axis,


Interpolation and upsampling

To interpolate data from the original Gauss-Legendre nodes on , to new Gauss-Legendre nodes on a subinterval , we use barycentric Lagrange interpolation [3]. By upsampling we refer to the special case of interpolating from to Gauss-Legedre nodes, both on .

Subinterval length criterion

When a new subpanel is formed on , we need to check if it satisfies the kernel-split accuracy criterion (30), which requires knowledge of the arc length of the subpanel, denoted . Assuming that does not vary rapidly on , a good approximation to is , where is the arc length of . We can now combine this approximation with (30), to get a an accuracy criterion formulated in subinterval size,


In particular, this allows us to write down the maximum length of subintervals on which product integration can be used,


4.1 Algorithm

Our algorithm proceeds with creating a division of into subintervals, which corresponds to a division of into subpanels.

For a target point with preimage such that , the first step is to create a subinterval centered on , with length set to satisfy both of the conditions (35) and (38). The centering ensures that the subpanels will not introduce new edges or quadrature nodes that are close enough to to degrade precision. If this initial subinterval has length , then the preimage of in that local frame will be , with and (36) is applicable.

If we wish the local preimage to be just beyond the limit where kernel-split quadrature is needed, then we must set such that . From (36), we can derive that this is satisfied when ,


For the subinterval to be contained within , it may not be bigger than twice the distance between and the closest interval edge,


Now we set the initial subpanel as large as possible, while still ensuring that the quadrature from it is accurate, and that it falls within ,


This gives us the initial subdivision . The center subinterval is now acceptable, and we proceed by recursively bisecting each remaining subinterval until either its length satisfies (37), or the local preimage of satisfies (35).

For target points such that , we can skip the process of carefully selecting the length of the nearest subinterval, and proceed immediately with recursive bisection of . This completes the algorithm, which we list in its entirety in algorithm 1.

function create_subdivision()
     if  then Preimage outside interval, recursively bisect all of it.
         return recursive_bisection().
          Center interval is now acceptable, recursively bisect remainder subintervals.
     end if
end function
function recursive_bisection()
     if  then
          From (33).
         if  and  then Using (34).
               Kernel-split must be used, but interval still too large. Continue bisection.
         end if
     end if
     return Subinterval passed.
end function
Algorithm 1 Given a panel of length and a nearby target point with preimage , create a subdivision of that allows the the layer potential to be accurately evaluated, using either direct Gauss-Legendre quadrature or kernel-split quadrature.

5 Numerical results

(a) Error when evaluating the layer potential close to the boundary, relative error in max norm computed as .
(b) Time to assemble the matrix blocks that correspond to neighboring panels, which is where the -correction need to be added.
Figure 4: Comparison of the original kernel-split algorithm, denoted “Original”, and our adaptive algortihm, denoted “With subdivision”, when solving our test problem for a large range of . We test the solution up to ; for larger values of the solution is identically zero in the entire domain.

To test the robustness of our method across a range of values, we setup the following test problem: We solve the modified Helmholtz equation (1) inside the annulus defined by a circle of radius 0.3, and a circle of radius 0.6, with Dirichlet boundary conditions given by a fundamental solution (18) located in the inner annulus,


The exact solution to this problem is equal to the expression for the Dirichlet boundary condition, evaluated in . To solve this problem numerically, we represent the solution using the double layer potential


where is the double layer kernel,


Enforcing the boundary condition gives us a second kind integral equation in ,


We solve this using the Nyström method, discretizing the boundary using -point Gauss-Legendre panels, with panels on the inner circle, and panels on the outer. For the bounds (37) and (39) we use use and (empirically determined). Following the notation of (5), the kernel-split of (46) is given by


For a large range of values, we solve the integral equation, and then evaluate the solution at 15 random points on a circle of radius 0.301 (very close to the inner boundary). The results, shown in fig. 4 demonstrate that our subdivision algorithm is capable of avoiding the catastrophic loss of accuracy otherwise present above a threshold , and that the additional cost incurred from it is proportional to . Around one digit of accuracy appears to be lost after the adaptive quadrature is activated, presumably due to the additional interpolation steps involved.

The present test problem is somewhat difficult to work with, since the solution to the modified Helmholtz equation decays exponentially fast. This limits the largest value of that we can test, as noted in the caption to fig. 4. To test larger ranges we can solve problems whose solutions do not decay as fast, such as the modified Stokes equations or the modified biharmonic equation. Such tests will be reported at a later date.

6 Conclusions

We present a robust recursive algorithm that allows the method of Helsing et al. to be applied, for any , to the modified Helmholtz equation, modified biharmonic equation and modified Stokes equations. Before, this was not possible for large , which corresponds to small timesteps with semi-implicit marching schemes for the heat equation and the time-dependent Stokes and Navier-Stokes equations. Our algorithm is fully adaptive, and the additional computational time it requires scales as . Our choice of the parameters and is based on numerical observations and provide excellent results. We will at a later date expand this report by adding an error analysis, together with error estimates that can be used for parameter selection.


The authors gratefully acknowledge the support from the Knut and Alice Wallenberg Foundation under grant no. 2016.0410 (L.a.K.), and by the Swedish Research Council under Grant No. 2015-04998 (F.F. and A.K.T.)


  • af Klinteberg and Tornberg [2017] L. af Klinteberg and A.-K. Tornberg. Error estimation for quadrature by expansion in layer potential evaluation. Adv. Comput. Math., 43(1):195–234, 2017. doi: 10.1007/s10444-016-9484-x.
  • af Klinteberg and Tornberg [2018] L. af Klinteberg and A.-K. Tornberg. Adaptive quadrature by expansion for layer potential evaluation in two dimensions. SIAM J. Sci. Comput., 40(3):A1225–A1249, 2018. doi: 10.1137/17M1121615.
  • Berrut and Trefethen [2004] J.-P. Berrut and L. N. Trefethen. Barycentric Lagrange interpolation. SIAM Rev., 46(3):501–517, 2004. doi: 10.1137/S0036144502417715.
  • Catherine A. Kropinski and Quaife [2011] M. Catherine A. Kropinski and B. Quaife. Fast integral equation methods for Rothe’s method applied to the isotropic heat equation. Comput. Math. with Appl., 61:2436–2446, 2011. doi: 10.1016/j.camwa.2011.02.024.
  • Chen et al. [2015] C. S. Chen, X. Jiang, W. Chen, and G. Yao. Fast solution for solving the modified Helmholtz equation with the method of fundamental solutions. Commun. Comput. Phys., 17(3):867–886, 2015. doi: 10.4208/cicp.181113.241014a.
  • Greengard and Kropinski [1998] L. Greengard and M. C. Kropinski. An integral equation approach to the incompressible Navier–Stokes equations in two dimensions. SIAM J. Sci. Comput., 20(1):318–336, 1998. doi: 10.1137/S1064827597317648.
  • Helsing and Holst [2015] J. Helsing and A. Holst. Variants of an explicit kernel-split panel-based Nyström discretization scheme for Helmholtz boundary value problems. Adv. Comput. Math., 41(3):691–708, 2015. doi: 10.1007/s10444-014-9383-y.
  • Helsing and Jiang [2018] J. Helsing and S. Jiang. On integral equation methods for the first Dirichlet problem of the biharmonic and modified biharmonic equations in nonsmooth domains. SIAM J. Sci. Comput., 40(4):A2609–A2630, 2018. doi: 10.1137/17M1162238.
  • Helsing and Ojala [2008] J. Helsing and R. Ojala. On the evaluation of layer potentials close to their sources. J. Comput. Phys., 227(5):2899–2921, 2008. doi: 10.1016/
  • Jiang et al. [2013] S. Jiang, M. C. A. Kropinski, and B. D. Quaife. Second kind integral equation formulation for the modified biharmonic equation and its applications. J. Comput. Phys., 249:113–126, 2013. doi: 10.1016/
  • Kropinski and Quaife [2011] M. C. A. Kropinski and B. D. Quaife. Fast integral equation methods for the modified Helmholtz equation. J. Comput. Phys., 230(2):425–434, 2011. doi: 10.1016/
  • [12] NIST. Digital library of mathematical functions. Release 1.0.16 of 2017-09-18. URL
  • Ojala and Tornberg [2015] R. Ojala and A.-K. Tornberg. An accurate integral equation method for simulating multi-phase Stokes flow. J. Comput. Phys., 298:145–160, 2015. doi: 10.1016/
  • Vorobjev [2019] Y. N. Vorobjev. Modeling of electrostatic effects in macromolecules. In L. A., editor, Comput. Methods to Study Struct. Dyn. Biomol. Biomol. Process. From Bioinforma. to Mol. Quantum Mech., pages 163–202. Springer International Publishing, Cham, 2019. ISBN 978-3-319-95843-9. doi: 10.1007/978-3-319-95843-9˙6.
  • Zhou and Pang [2018] H.-X. Zhou and X. Pang. Electrostatic interactions in protein structure, folding, binding, and condensation. Chem. Rev., 118(4):1691–1741, 2018. doi: 10.1021/acs.chemrev.7b00305.