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

TL;DR

This paper introduces an adaptive quadrature algorithm that maintains high accuracy for parameter-dependent layer potentials across a wide range of parameters, overcoming limitations of existing methods.

Contribution

It presents a novel adaptive sampling approach with recursive bisection to improve kernel-split quadrature for parameter-dependent equations.

Findings

Achieves accurate evaluation for a broader range of the parameter 

Maintains accuracy with a computational cost scaling as log

Extends the applicability of kernel-split quadrature to complex problems

Abstract

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 $\alpha$, and kernel-split quadrature loses its accuracy rapidly when this parameter grows beyond a certain threshold. This paper 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 the parameter $\alpha$, at an increased cost that scales as $\log\alpha$. Using this algorithm allows kernel-split quadrature to be both accurate and efficient for a much wider range of problems than previously possible.

Equations

$$(\Delta -{\alpha }^2)u=0,$$

$$\Delta (\Delta -{\alpha }^2)u=0,$$

$$(\Delta -{\alpha }^2)u-\nabla p=0,$$

$$u(x)={\int }_{\partial \Omega }G(x,y)\sigma (y)dS(y).$$

$$G(x,y)=G^S(x,y)+G^L(x,y)\log \left|y-x\right|+G^C(x,y)\frac{(y-x)\cdot \hat{n}(y)}{{\left|y-x\right|}^2},$$

$$\partial \Omega =\bigcup _i{\Gamma }_i.$$

$${\Gamma }_i=\left\{{\gamma }_i(t)\in {\mathbb{R}}^2|t\in [-1,1]\right\}.$$

$${\int }_{{\Gamma }_i}G(x,y)\sigma (y)dS(y)$$

$$(9)y_{ij}$$

$$(12)s_{ij}$$

$${\int }_{\Gamma }G(x,y)\sigma (y)dS(y)$$

$${\int }_{\Gamma }f(x,y)\log \left|y-x\right|dS(y)$$

$${\int }_{\Gamma }f(x,y)\frac{(y-x)\cdot \hat{n}(y)}{{\left|y-x\right|}^2}dS(y)$$

$$\underset{\begin{array}{c}x\to y\\ x,y\in \partial \Omega \end{array}}{\lim }\frac{(y-x)\cdot \hat{n}(y)}{{\left|y-x\right|}^2}=-\frac{\kappa (y)}{2},$$

$$\begin{array}{rl}{\int }_{\Gamma }G(y_i,y)\sigma (y)dS(y)\approx & \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n\left[G(y_i,y_j){\sigma }_j{s}_j{\lambda }_j^G+G^L(y_i,y_j){\sigma }_j\left(w_j^L(y_i)-s_j\log \left|y_j-y_i\right|{\lambda }_j^G\right)\right]\\ & +G^S(y_i,y_i)+G^L(y_i,y_i){\sigma }_i{w}_i^L(y_i)-G^C(y_i,y_i)\frac{\kappa (y_i)}{2}.\end{array}$$

$$\begin{array}{rl}{\int }_{\Gamma }G(x,y)\sigma (y)dS(y)\approx & \sum _{j=1}^n\left[G(x,y_j){\sigma }_j{s}_j{\lambda }_j^G+G^L(x,y_j){\sigma }_j\left(w_j^L(x)-s_j\log \left|y_j-x\right|{\lambda }_j^G\right)\right].\end{array}$$

$$\displaystyle\begin{split}\int_{\Gamma}G(x,y)\sigma(y)\operatorname{d\!}S(y)\approx\sum_{j=1}^{n}&\Bigg{[}G(x,y_{j})\sigma_{j}s_{j}\lambda_{j}^{G}\\ +&G^{L}(x,y_{j})\sigma_{j}{\left(w^{L}_{j}(x)-s_{j}\log\mathinner{\!\left\lvert y_{j}-x\right\rvert}\lambda_{j}^{G}\right)}\\ +&G^{C}(x,y_{j})\sigma_{j}{\left(w^{C}_{j}(x)-s_{j}\frac{(y_{j}-x)\cdot{\hat{n}}(y_{j})}{\mathinner{\!\left\lvert y_{j}-x\right\rvert}^{2}}\lambda_{j}^{G}\right)}\Bigg{]}.\end{split}$$

$$G(x,y)=K_0(\alpha |y-x|),$$

$$K_0(\rho )=K_0^S(\rho )-I_0(\rho )\log \rho ,\rho \in {\mathbb{R}}^{+},$$

$$I_{\nu }(\rho )={\left(\frac{1}{2}\rho \right)}^{\nu }\sum _{k=0}^{\infty }\frac{{\left(\frac{1}{4}{\rho }^2\right)}^k}{k!(\nu +k)!},\nu \in \mathbb{N},\rho \in \mathbb{R},$$

$$I_{\nu }^{(n)}(\rho )=2^{-n}\sum _{m=0}^n\left(\genfrac{}{}{0ex}{}{n}{m}\right)I_{2m+(\nu -n)}(\rho ),\nu ,n\in \mathbb{N},\rho \in \mathbb{R},$$

$$G^S(x,y)$$

$$G^L(x,y)$$

$$G^C(x,y)$$

$$\alpha h\le C_{\epsilon },$$

$$h\le h_{\text{max}}\text{mathchar}58=C_{\epsilon }/\alpha .$$

$$u(x)={\int }_{-h/2}^{h/2}{K}_0(x,y)dy.$$

$$u(x)\approx \sum _{j=1}^n\left(G(x,y_j)-G^L(x,y_j)\log \left|y_j-x\right|\right){\lambda }_j-{\int }_{-h/2}^{h/2}{I}_0(\alpha |y-x|)\log |y-x|dy,$$

$$I_0(\alpha |y-x|)$$

$$p_{n-1}(x,y)$$

$$p_{n-1}(x,y_j)$$

$${\int }_{-h/2}^{h/2}{I}_0(\alpha |y-x|)\log |y-x|dy=\mathrm{Re}\left({\int }_{-h/2}^{h/2}{I}_0(\alpha |y-x|)\log (y-x)dy\right)$$

$$=\mathrm{Re}\left(\sum _{k=0}^{n-1}{c}_k(x)q_k(x)\right)+\mathrm{Re}\left({\int }_{-h/2}^{h/2}{r}_n(x,y)\log (y-x)dy\right).$$

$$q_k(x)={\int }_{-h/2}^{h/2}{y}^k\log (y-x)dy,k=0,\dots ,n-1,$$

$$R_n(\alpha ,x,h)={\int }_{-h/2}^{h/2}{r}_n(x,y)\log |y-x|dy=\text{Re}\left({\int }_{-h/2}^{h/2}{r}_n(x,y)\log (y-x)dy\right).$$

Full text

Keywords— Integral equations, Partial differential equations, Layer potentials, modified Helmholtz equation, modified Stokes equation

1 Introduction

This paper presents an extension of the panel-based, kernel-split quadrature scheme by Helsing and Ojala [13], for evaluating singular and nearly singular layer potentials in two dimensions. It is one of the current state of the art methods for maintaining low errors when solving homogeneous elliptic partial differential equations (PDEs) in two dimensions using integral equation methods [13, 11, 17]. However, there exists a set of problems for which this scheme can fail completely. This includes the following PDEs in $\mathbb{R}^{2}$:

[TABLE]

where $\alpha$ is a positive 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 [20, 19, 14, 15, 7, 12].

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 $\alpha$. This decay presents a problem for the above mentioned 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 locally well approximated by polynomials, which is increasingly difficult for larger values of $\alpha$.

Large values of $\alpha$ are of interest in the context of elliptic marching. In elliptic marching a semi-implicit temporal discretization is applied to the governing equations. A time-step then involves solving a sequence of elliptic equations, such as the modified PDEs: e.g. the heat equation, the time-dependent Stokes and Navier-Stokes equations correspond to the modified Helmholtz equation, the modified biharmonic equation and the modified Stokes equations, respectively [6, 8, 9, 3]. Regardless of the specifics of the time-discretization scheme the resulting equations involve the parameter $\alpha$, where $\alpha^{2}$ is inversely proportional to the time-step size. The smaller the time-step the lower temporal error, but the spatial problem becomes harder to solve accurately. This is a new problem in the context of elliptic marching, combined with boundary integral methods applied to the resulting modified PDEs. Earlier, high order and accurate schemes were not available, thus small time-steps were redundant since the spatial error dominated, and the low-error regimes were not feasible to explore on standard desktop computers. Today efficient methods are available that give higher precision and allow for more complicated problems. It is evident that the restriction of the time-step is a significant bottleneck, as the kernel-split quadrature may fail for high temporal resolutions.

We have developed a robust quadrature scheme, based on adaptive refinement, that maintains high accuracy for any $\alpha$, without sacrificing efficiency. It applies to 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. The additional cost, in terms of assembly time per panel, scales as $\mathcal{O}(\log\alpha)$.

The remainder of this paper 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. In section 4 we present error analysis, again with the modified Helmholtz equation as template. The new algorithm we propose is presented in section 5, followed by numerical results in section 6.

2 Background

Our goal is to evaluate layer potentials in $\mathbb{R}^{2}$ of the form

[TABLE]

The layer density $\sigma(y)$ is assumed to be smooth, implying that the boundary is assumed to be smooth as well. The kernel $G(x,y)$ is singular at $x=y$ and can be expressed with explicit singularities as

[TABLE]

where $G^{S}$ is a smooth function. The functions $G^{L}$ and $G^{C}$ are smooth functions that multiply a log-type singularity and a Cauchy-type singularity, respectively. We refer to this decomposition as kernel-split.

The boundary ${\partial\Omega}$ is discretized using a composite Gauss–Legendre quadrature. It is subdivided into intervals $\Gamma_{i}$, denoted panels,

[TABLE]

Each panel $\Gamma_{i}$ is described by a parametrization $\gamma_{i}$,

[TABLE]

We refer to $x$ in (4) as a target point and a point $y$ as a source point. A panel to which a source point belongs is referred to as a source panel. Associated with the parametrization is a speed function $s_{i}(t)=\mathinner{\!\left\lvert\gamma_{i}^{\prime}(t)\right\rvert}$, a normal vector ${\hat{n}}_{i}(t)$ and the curvature $\kappa_{i}(t)$. Introducing the convenience notation $\sigma_{i}(t)=\sigma(\gamma_{i}(t))$, the layer potential from a panel $\Gamma_{i}$ becomes

[TABLE]

Each panel is discretized in the parametrization variable $t$ using the nodes and weights $(t_{j}^{G},\lambda_{j}^{G})$ of an $n$-point Gauss–Legendre quadrature rule on the canonical interval $[-1,\,1]$, which is of order $2n$, such that on each panel we have the discrete quantities

[TABLE]

Omitting the panel index $i$, the layer potential contribution from a panel $\Gamma$ is then computed using the approximation

[TABLE]

Due to the singularities in $G$, the above formula requires $x$ to be well-separated from $\Gamma$ in order to be accurate (see section 5, paragraph Near evaluation criterion). Otherwise, the scheme of [11] is used, known as product integration. With that, target-specific quadrature weights $w^{L}$ and $w^{C}$ of order $n$ are computed for the known singularities where needed, such that

[TABLE]

Substituting (5) into (4) and applying the above product integration gives the so called kernel-split quadrature scheme. Depending on the location of the target point $x$ relative to the source panel $\Gamma$, the evaluation can be divided up into three different cases:

Singular, with self-interaction. If $x\in\Gamma$ is one of the quadrature nodes, $x=y_{i}$, then the term multiplying $G^{C}$ is smooth, with the limit

[TABLE]

where $\kappa(y)$ is the curvature of $\partial\Omega$ at $y$. Applying product integration to the $G^{L}$ term, we get

[TABLE]

Note that this is the only case where $G^{S}$ is evaluated explicitly, and only in one point. 2. 2.

Singular, without self-interaction. If $x\in{\partial\Omega}$ is either on the source panel $\Gamma$ but not a quadrature node, or on a neighboring panel, then the $G^{C}$ term is still smooth, and we do not have to take the limit at $x\to y$ into account. This lets us simplify the above to

[TABLE] 3. 3.

Nearly singular case. If $x$ is close to $\Gamma$, but not on a neighboring panel, then we need to address both singularities in (5). This case occurs when either $x\in\Omega$ is close to ${\partial\Omega}$, or when $x\in{\partial\Omega}$ is on a section of the boundary that is distant in arc length or disjoint. The layer potential is then evaluated as

[TABLE]

Given a target point $x$, the panels on ${\partial\Omega}$ are partitioned into two sets: far panels, that can be evaluated directly using (14), and near panels, that must be evaluated (or corrected), using either (18), (19), or (20).

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 [13, 11, 17] on geometries with a smooth boundary. However, for the modified PDEs (1), (2), and (3), the scheme can fail completely. All of these equations have layer potential kernels including second-kind modified Bessel functions, of the forms $K_{0}(\alpha|y-x|)$ and/or $K_{1}(\alpha|y-x|)$. As we shall see, this is problematic for the kernel-split quadrature.

As an illustrating example, we study the single layer kernel of the modified Helmholtz equation (1). The associated Green’s function is

[TABLE]

where the scaling factor $1/2\pi$ has been omitted. The split of this kernel is based on using a standard decomposition [16, §10.31] to explicitly write out the singularities in $K_{0}$,

[TABLE]

where $K_{0}^{S}$ is the smooth remainder, and $I_{0}$ is the modified Bessel function of the first kind. For future reference, we have by [16, §10.25,§10.29] that

[TABLE]

and the derivatives of $I_{\nu}$ can be expressed as

[TABLE]

which makes them easy to evaluate using standard numerical libraries,

Inserting (22) into the kernel (21), after also splitting the logarithm, we identify the terms in (5) as

[TABLE]

Product integration is a semi-analytical method. This means that the target-specific quadrature weights in (15) are found by an analytic treatment of the singularity or near singularity. The accuracy of the method for a given evaluation point $x$ is basically limited by how well a function $f(x,y)$ can be resolved by an ($n-1$)th-degree Legendre polynomial. For a well-resolved function, the integration error is usually only at most a few orders of magnitude larger than round-off, independent of the location of the evaluation point. We now want to evaluate the layer potential (4) with the layer density $\sigma$, using the split above. With this $G^{L}$, $f(x,y)$ in (15) will be $-I_{0}(\alpha|y-x|)\sigma(y)$. Hence, for good accuracy in integrating the logarithmic singularity, it is this product that must have a small interpolation error.

The function $G(x,y)$ goes to zero as $|y-x|$ goes to infinity, since $K_{\nu}(\rho)\sim\sqrt{\pi/(2\rho)}e^{-\rho}$ as $\rho\to\infty$ for $\nu=1,2$, but its components $G^{S}(x,y)$ and $G^{L}(x,y)$ actually grow exponentially as $e^{\alpha|y-x|}$, with opposing signs, following the asymptotic $I_{\nu}(\rho)\sim e^{\rho}/\sqrt{2\pi\rho}$ as $\rho\to\infty$ for $\nu=1,2$ [16, §10.25,§10.30]. As $\alpha$ gets larger, this makes $G^{L}$ an increasingly bad candidate for polynomial interpolation, which will render the product integration inaccurate. In addition, when evaluated in limited precision, this split is prone to numerical cancellation, due to the limiting forms of $K_{0}$ and $I_{1}$. We will argue that the error mainly is a function of the quantity $\alpha h$, which can be heuristically motivated by figure 1. To ensure that this error remains below some tolerance $\epsilon$, we therefore suggest that $\alpha$ and panel length $h$ must satisfy a criterion on the form

[TABLE]

for some constant ${C_{\varepsilon}}$, which can be determined using numerical experiments. The above criterion can also be reformulated as follows: In order to achieve a tolerance $\epsilon$, panel lengths must satisfy

[TABLE]

The naive way of achieving this, for a given $\alpha$, is to discretize ${\partial\Omega}$ 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. Global refinement is also redundant for another reason: as stated above the functions $K_{0}(\rho)$ and $K_{1}(\rho)$ decay as $\sqrt{\pi/(2\rho)}e^{-\rho}$, meaning that they are very localized for large $\alpha$, i.e. when polynomial interpolation might fail, and are almost zero in finite precision away from $\rho=0$. Thus only a small portion of the boundary needs refinement.

4 Error estimates

As was indicated above, there are two main sources of errors in the kernel split quadrature that grow as $\alpha$ increases. The first part arises from an interpolation error, and the second is due to numerical cancellation. In order to better understand these errors, we will perform a limited analysis for the single layer potential of the modified Helmholtz equation on a flat panel in section 4.1.

This understanding of the error structure is useful also for the associated kernels for the other modified PDEs since they have a similar singularity structure to the single layer modified Helmholtz potential, see appendix B. They are all combinations of $K_{0}$ and/or $K_{1}$, thus in the form of explicit singularities they contain $I_{0}$ and/or $I_{1}$ multiplying a log-type singularity.

Based on our simple analysis in section 4.1, we will in section 4.2 formulate an error estimate and illustrate its performance for the single layer modified Helmholtz potential. The error estimates can be used to set $C_{\epsilon}$ in (28), as discussed in section 4.4. Other kernels are discussed in section 4.3.

In this section as well as the next, we will identify vectors $x,y\in\mathbb{R}^{2}$ with the corresponding complex numbers $x,y\in\mathbb{C}$.

4.1 Error analysis for a flat panel

We consider the single layer potential for the modified Helmholtz equation with unit density, evaluated using kernel-split quadrature from a flat panel of length $h$,

[TABLE]

To evaluate (30) using the kernel-split correction (20), we first write

[TABLE]

using $G^{C}(x,y)=0$ (27). In order to evaluate the remaining integral, the function $I_{0}(\alpha|y-x|)$ is approximated by an $(n-1)$th-degree polynomial $p_{n-1}$ that interpolates $I_{0}(\alpha|y-x|)$ at the $n$ nodes $y_{j}$,

[TABLE]

where $r_{n}$ is the polynomial interpolation error. We consider only even values for $n$. The integral on the right-hand side of (31) can be written as a sum: the approximation by integrating the polynomial $p_{n-1}$, and an integral over the polynomial interpolation error

[TABLE]

Here, we have used the definition of the complex logarithm. We have that

[TABLE]

are integrals that can be computed recursively, starting from exact formulas [10].

The error in the kernel-split quadrature has two sources. The first is the integral over the interpolation error from the integral in the right-hand side of (31), namely

[TABLE]

For brevity we also refer to the integral over the interpolation error as the interpolation error. It will be clear from context which one that is referred to.

The second is a cancellation error due to the limiting forms of $G$ and $G^{L}$. Recall that $G^{S}=G-G^{L}$; this function grows exponentially as $-e^{\alpha|y-x|}/\sqrt{2\pi\alpha|y-x|}$, which is the same asymptotics as for $G^{L}$, but with opposite sign. The function $G$ is a decaying function as $\alpha|y-x|$ goes to infinity, thus $G^{S}$ and $G^{L}$ have to cancel. As the quadrature terms in (31) grows in magnitude catastrophic cancellation errors follow.

We will now perform the analysis of the interpolation error for target points $x$ along the real axis. In section 4.2, we will consider the performance of this estimate for target points in the plane.

4.1.1 Interpolation error

Our goal is to estimate the magnitude of the interpolation error $R_{n}(\alpha,x,h)$ in (38), and we will do so for $x\in\mathbb{R}$. The polynomial $p_{n-1}$ (33) is formed by interpolating at the Legendre nodes, $y_{j}=t_{j}^{G}h/2$. According to [18], the remainder is given by

[TABLE]

where the absolute value of the argument of $I_{0}$ can be removed, since $I_{0}$ is an even function. This is convenient, as it simplifies differentiation.

The $n$ Legendre nodes are the roots of the Legendre polynomial of order $n$, denoted $P_{n}$, orthogonal on $[-h/2,\,h/2]$. Hence

[TABLE]

where $\ell_{n}$ is the coefficient of the leading order monomial term in $P_{n}$. From Rodrigues’ formula [16, §18.5], it can be shown that $\ell_{n}=2^{-n}(2n)!/(n!)^{2}$. After rescaling $\xi$ we have

[TABLE]

To present an error bound for (41) we can apply

[TABLE]

However, the resulting error bound greatly overestimates the error, thus is not of practical use. We pursue an alternative approach by estimating $r_{n}$ by choosing a value for $\xi$ for all pairs $(x,y)$. For this purpose, heuristics suggest $\xi=0$, as demonstrated below.

Combining (38) and (41) for $\xi=0$ we get the estimate

[TABLE]

with $\Lambda_{S}$ defined as

[TABLE]

where the subscript $S$ denotes single layer potential. Note that $I_{0}^{(n)}$ is an even function for even values of $n$; by the relation (24) $I_{0}^{(n)}$ is a linear combination of $I_{\nu}$ with even $\nu$, and it follows from the definition (23) that these terms are even. Thus we may take the absolute value of the argument of $I_{0}^{(n)}$ without altering the result. This is important, as later on we will consider complex target points $x$.

In (44) we have

[TABLE]

obtained by

[TABLE]

where the integral is rewritten using $\log(h/2(\rho-2x/h))=\log(h/2)+\log(\rho-2x/h)$ and the fact that each $P_{n}(\rho)$ integrates to zero over the interval.

This “Legendre-log integral” can be evaluated using the recursion formulas in section A.

For $h=2$, both the absolute value as well as the real part of $\Lambda_{S}$ in (44) are shown as a function of $x$ in figure 2. The absolute value of $\Lambda_{S}$ is a smooth function that achieves its maximum for $x$ close to the endpoints of the interval. This is due to the $\log$-factor in the integrand being singular for target points $x$ in $[-1,1]$, combined with the growth of the factor $I_{0}$ towards the edges of the panel. For target points $x$ outside the source panel the Legendre-log integral $g_{n}$ decreases rapidly.

This plot is only illustrating the error estimate derived for this simplifed case. How well it actually estimates the error will be discussed after we have considered also the second part of the error.

4.1.2 Numerical cancellation error

In addition to the interpolation error just discussed, the kernel-split (31) can also suffer from numerical cancellation when $\alpha$ is large. To see this, note that $G(x,y)-G^{L}(x,y)\log(\alpha\mathinner{\!\left\lvert y-x\right\rvert})=G^{S}(x,y)$, which has the asymptotics $-e^{\alpha\mathinner{\!\left\lvert y-x\right\rvert}}/\sqrt{2\pi\alpha\mathinner{\!\left\lvert y-x\right\rvert}}$ as $\alpha\mathinner{\!\left\lvert y-x\right\rvert}$ goes to infinity. This is the same asymptotics as for $I_{0}$, but with opposing sign. These two terms have to cancel, since $K_{0}$, i.e. $G$, is a decaying function. The sum of terms with large magnitude, but opposing signs, is prone to numerical cancellation in limited precision.

The cancellation error is straightforward to estimate by

[TABLE]

where $\epsilon_{\mathrm{mach}}$ is the machine epsilon.

4.2 Error estimate

We now have an error estimate $|\Lambda_{S}(\alpha,n,x,h)|$ (44) for the interpolation error $R_{n}(\alpha,x,h)$ (38), and an estimate $\Xi_{S}(\alpha,x,h)$ (47) for the cancellation error. Both were derived for target points along the real axis, but we now evaluate the estimated maximum error over a set $D$ of discrete target points in the plane, and compare it to the maximum measured actual error. There are two goals: to study the errors’ and error estimates’ dependence on $\alpha h$, and to construct a combination of the error estimates to predict the total error.

Let $D$ be the set $100$ sampled complex points with positive imaginary part uniformly within a Bernstein ellipse with foci $\pm 1$, and denote its elements as $\tilde{x}$. The radius for the Bernstein ellipse is set to be equal to $3^{16/n}$, as points farther away from the panel do not need the kernel-split quadrature scheme to be accurate [4]. This set is visualized in figure 3. By setting $\tilde{x}=2x/h$, the estimate (44) can be rewritten as

[TABLE]

and similarly the cancellation error (47) becomes

[TABLE]

We define

[TABLE]

The maximum of these two quantities gives the error estimate

[TABLE]

for given values of $\alpha$, $n$ and $h$. From (48) and (49) it clear that $E_{S}$ is a function with $\alpha h$ and $h$ as two separate arguments, motivating the formulation

[TABLE]

Hence, we expect the error divided by $h$ to be a function of $\alpha h$ only, and this is confirmed in figure 3 and figure 4. Here, we plot the error scaled with $1/h$ for different values of $h$ versus $\alpha h$, for different values of $n$. We see that the error curves collapse for different values of $h$, just as predicted. Furthermore, from the plots it is clear that the interpolation error and cancellation error dominate in different regimes. As $n$ is increased, the interpolation error decreases, and the cancellation error will continue to be dominant for larger values of $\alpha h$. Considering the left plot of figure 4, we can most clearly see how the measured error follows the cancellation error estimate $\Xi_{S,\max}$ first, and then shifts to follow the interpolation error estimate $\Lambda_{S,\max}$ as it becomes dominant. We can conclude that the error estimate $\tilde{E}_{S}$ predicts the actual error, scaled with $1/h$, quite well.

4.3 Other kernels

We now discuss error estimates for the other kernels of interest for the modified PDEs. The double layer kernel for modified Helmholtz is discussed in detail, which is used to motivate a similar treatment of the error estimates for the modified biharmonic equation and the modified Stokes equations. Their corresponding kernels, and the decompositions thereof into explicit singularities, are shown in appendix B.

The kernels for layer potentials associated with the modified Helmholtz equation (double layer), modified biharmonic and modified Stokes are harder to study analytically than the single layer kernel for the modified Helmholtz equation; $G^{L}$ for all of them contain $I_{0}$ and/or $I_{1}$ multiplying a smooth function, which we collectively denote $F$. This means that the interpolation error (39) includes the $n$th derivative of a product. However, we observe that the derived error estimates of the interpolation error for the single layer kernel can be applied to the these kernels as well, with a few modifications. This is not surprising, as the difficulties lie in properly resolving $I_{0}$ and/or $I_{1}$, not resolving the function $F$. The cancellation error (51) is also straightforward to modify.

Before proceeding with the error estimates we present the split of the double layer kernel for the modified Helmholtz equation into explicit singularities. The associated kernel is, again omitting the scaling $1/2\pi$,

[TABLE]

and by [16, §10.31] one has

[TABLE]

The content of the functions $G^{0}$, $G^{L}$ and $G^{C}$ depend on the location of the target point $x$. If $x$ belongs to the boundary, i.e. $x\in\partial\Omega$, then the kernel-split is

[TABLE]

For off-boundary $x$ the kernel-split is

[TABLE]

The main difference between the two decompositions is that $1/|y-x|$ in (55) goes into $G^{S}$ for $x\in\partial\Omega$, due to the known limit value (17), while for $x\notin\partial\Omega$ it results in a non-zero $G^{C}$. For both kernel-splits the function $G^{L}$ is smooth and consists of the modified Bessel function $I_{1}$ multiplied with $F=(y-x)\cdot\hat{n}(y)/|y-x|$. Here, $F$ does not pose the difficulties with polynomial interpolation as $I_{1}$ does.

We now create an error estimate for the interpolation error. Again identifying the vectors $x$ and $y$ in $\mathbb{R}^{2}$ with complex numbers $x$ and $y$ in $\mathbb{C}$, and following the same steps as for the single layer potential we get analogous to (41)

[TABLE]

since $I_{1}$ is an odd function. The derivation for the single layer kernel that we based this result on was done for $x$ in $\mathbb{R}$. We will now, as before, apply it to target points $x$ in $\mathbb{C}$, for which ${\left(\xi h/2-x\right)}\cdot\hat{n}(\xi h/2)$ is non-zero. To avoid the $n$th derivative of the product we make the following argument. We exclude the value $x=\xi h/2$ since then $r_{n}$ is equal to zero. Away from zero the derivative of the sgn function is zero. Furthermore, it is $I_{1}$ that is difficult to represent accurately with polynomials and it contributes considerably more to the error than $F$ does, thus we keep only the dominant term obtained by the chain rule and write

[TABLE]

As for the single layer we set $\xi$ to zero, which results in the estimate

[TABLE]

where $g_{n}(\tilde{x})$ is given by (45) and $\tilde{x}=2x/h$ as before. Again we have reinserted the absolute value in the argument of the modified Bessel function, by removing the $\mathrm{sgn}$ function, since by (23) and (24) $I_{1}^{(n)}$ is an odd function for even $n$.

For the cancellation error, similarly

[TABLE]

where $(y-x)\cdot\hat{n}(y)/|y-x|$ has been replaced with $\mathrm{Im}(\tilde{x})/|\tilde{x}|$.

Clearly the error estimates depend purely on $\alpha h$, without any multiplying $h$-factor, in difference to the single layer kernel. For the double layer potential the error estimate is

[TABLE]

where $\Lambda_{D,\max}$ and $\Xi_{D,\max}$ are defined analogous to (50) and (51) for (65) and (66). The predicted pure dependence on $\alpha h$ for the measured numerical error can be seen in figure 5. The estimate for the cancellation error is not as tight as for the single layer kernel for $n=32$. Still, at its worst the estimates are off by a digit of the actual error. Also, we see that it is possible to ignore the chain rule and only differentiate $I_{1}$, but not $F$, and obtain good estimates.

To derive error estimates for the kernels associated with the other modified PDEs a similar approach as above can be taken. The error estimate $E$ will have the form

[TABLE]

where simply $H(h)=h^{p}$ for some $p$.

The discussion that we have put forward has given an understanding of the two sources of errors. Practically, only a value for ${C_{\varepsilon}}$ needs to be determined; it sets panel lengths in the subdivision algorithm, introduced in the next section. Instead of deriving an error estimate, one can simply use figures such as figure 3–5 for the measured numerical errors for this purpose. Here, one needs to find the proper scaling $H(h)$ by which to scale the errors before plotting, such that they coincide for different $h$. In figure 6 such results are presented for the double layer kernel for the modified Stokes equations (106), with $p=1$. The error curves do not collapse as neatly as for the modified Helmholtz equation, indicating that there is a scaling factor other than $h$. Still, the results can be used to set ${C_{\varepsilon}}$, as demonstrated in section 6.

4.4 Condition for $\alpha h$

As was shown above, the error when evaluating the double layer modified Helmholtz kernel with an $n$ point Gauss–Legendre quadrature rule over one panel with length $h$ only depends on $n$ and $\alpha h$. That means that given the error tolerance $\varepsilon$, ${C_{\varepsilon}}=\alpha h$ that yield this error can be read off from a plot such as in figure 5 for the given $n$. Fixing that ${C_{\varepsilon}}$, the condition (as written in (28)) is $\alpha h\leq{C_{\varepsilon}}$.

For other kernels, such as the single layer modified Helmholtz kernel, there is not a pure $\alpha h$ dependence on the error, as we have seen. However, to keep the procedure simple, we would suggest the following approach. Pick a typical value of the panel length $h$ for the discretization at hand. Then, divide the given tolerance $\varepsilon$ by this $h$, and use this scaled tolerance when determining ${C_{\varepsilon}}$ from plots of the scaled error for the appropriate $n$, such as in figure 3 and figure 4.

In the subdivision algorithm, the reduction in error through this extra factor of $h$ of refinement will then not be taken into account, and hence the panels might be refined to be somewhat smaller than needed. But this will allow us to keep the same simple structure of the subdivision algorithm for all different kernels.

5 Quadrature by recursive subdivision

To circumvent the problem of interpolatory quadrature failing for large $\alpha$, we here introduce an algorithm for local refinement, based on panel subdivision. Given a single source panel $\Gamma$, we assume that it is sufficiently short, relative to the quadrature order $n$, for both the geometry and the layer density to be well-represented by a polynomial, interpolated at the $n$ quadrature points. We say that it is well-resolved. For a given target point $x$, we can then subdivide $\Gamma$ into a set of $M$ subpanels $\{\Gamma_{i}\}_{i=1}^{M}$, interpolate our known quantities from the quadrature nodes on $\Gamma$ to the quadrature nodes on those subpanels, and then evaluate the layer potential at $x$ using the subpanels. To ensure accuracy for a given tolerance, this subdivision is formed in a way that guarantees that all subpanels either are short enough to satisfy (28), obtained via (52) or similarly, or are sufficiently far away from $x$, 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 $\gamma(t)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\mathbb{C}$ be the mapping from the standard interval $[-1,1]$ to the panel $\Gamma$. Then $z$, such that $\gamma(z)=x$, is the preimage of the target point $x$. The preimage $z$ is real-valued if $x\in{\partial\Omega}$, and complex-valued otherwise. We here assume that we know the value of $z$, but $\gamma(t)$ need not be a known function; see [2] for a discussion on how construct a numerical representation.

Subpanels and subintervals

A subdivision of $\Gamma$ is defined by a set of edges in the parametrization, $\{-1=t_{1},t_{2},\dots,t_{M+1}=1\}$, such that a subpanel $\Gamma_{i}$ is given by the mapping of the subinterval $[t_{i},t_{i+1}]$ under $\gamma$. We can, by a linear scaling, define the local mapping that maps the standard interval $[-1,1]$ to $\Gamma_{i}$ as

[TABLE]

where $\Delta t_{i}=t_{i+1}-t_{i}$. Given the preimage $z$, the local preimage $z_{i}$, such that $\gamma_{i}(z_{i})=x$, is given by

[TABLE]

Near evaluation criterion

Given the preimage $z$ of a point $x$ close to a panel (or subpanel) $\Gamma$, it is possible to compute an accurate estimate of the quadrature error when evaluating the layer potential using $n$-point Gauss–Legendre quadrature. Detailed discussions can be found in [1, 2]. To leading order, the error is proportional $\rho(z)^{-2n}$, where $\rho$ is the elliptical radius of the Bernstein ellipse on which $z$ lies,

[TABLE]

where $\sqrt{z^{2}-1}$ is defined as $\sqrt{z+1}\sqrt{z-1}$ with $-\pi<\arg{\left(z\pm 1\right)}\leq\pi$. For a given kernel $G$ and error tolerance $\varepsilon$, it is then possible to introduce a cutoff radius ${R_{\varepsilon}}$, such that kernel-split quadrature must be used for

[TABLE]

and otherwise Gauss–Legendre quadrature is sufficiently accurate. We refer to this as the target point and the source points being well-separated, otherwise they are considered to be close. See [2] on how to set ${R_{\varepsilon}}$ for a given tolerance; we use ${R_{\varepsilon}}=3.5$ for $n=16$, which corresponds to a tolerance of $10^{-14}$.

Later we will use that the inverse of $\rho$ has a particularly simple form in the special case when $z$ lies on the imaginary axis,

[TABLE]

Interpolation and upsampling

To interpolate data from the $n$ original Gauss–Legendre nodes on $[-1,1]$, to $m$ new Gauss–Legendre nodes on a subinterval $[t_{i},t_{i+1}]\subset[-1,1]$, we use barycentric Lagrange interpolation [5]. By upsampling we refer to the special case of interpolating from $n$ to $2n$ Gauss–Legedre nodes, both on $[-1,1]$.

Subinterval length criterion

When a new subpanel $\Gamma_{i}$ is formed on $\Gamma$, we need to check if it satisfies the kernel-split accuracy criterion (28), which requires knowledge of the arc length of the subpanel, denoted $h_{i}$. Assuming that $\gamma^{\prime}(t)$ does not vary rapidly on $\Gamma$, a good approximation to $h_{i}$ is $h_{i}\approx h\Delta t_{i}/2$, where $h$ is the arc length of $\Gamma$. We can now combine this approximation with (28), to get an accuracy criterion formulated in subinterval size,

[TABLE]

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

[TABLE]

Here, ${C_{\varepsilon}}$ can be obtained via (52), as a precomputation step.

5.1 Algorithm

Our algorithm proceeds with creating a division of $[-1,1]$ into subintervals, which corresponds to a division of $\Gamma$ into subpanels.

For a target point $x$ with preimage $z$ such that $\operatorname{Re}z\in(-1,1)$, the first step is to create a subinterval centered on $\operatorname{Re}z$, with length set to satisfy both of the conditions (72) and (75). The centering ensures that the subpanels will not introduce new edges or quadrature nodes that are close enough to $z$ to degrade precision.

If this initial subinterval has length $\Delta t_{c}$, then the preimage of $x$ in that local frame will be $z_{c}={\mathrm{i}}b_{c}$, with $b_{c}=2\operatorname{Im}z/\Delta t_{c}$ and (73) is applicable.

If we wish the local preimage $z_{c}$ to be just beyond the limit where kernel-split quadrature is needed, then we must set $\Delta t_{c}$ such that $\rho(z_{c})={R_{\varepsilon}}$. From (73), we can derive that this is satisfied when $\Delta t_{c}=\Delta t_{\text{direct}}$,

[TABLE]

For the subinterval to be contained within $[-1,1]$, it may not be bigger than twice the distance between $\operatorname{Re}z$ and the closest interval edge,

[TABLE]

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 $[-1,1]$,

[TABLE]

Here, $\Delta t_{\text{direct}}$ depends on the location of the target point, while ${\Delta t_{\max}}$ does not. Thus, for points sufficiently far away ${\Delta t_{\max}}$ might be larger than $\Delta t_{\text{direct}}$.

By (78) we have the initial subdivision $\{-1,\operatorname{Re}z-\Delta t_{c}/2,\operatorname{Re}z+\Delta t_{c}/2,1\}$. The center subinterval is now acceptable, and we proceed by recursively bisecting each remaining subinterval until either its length satisfies (74), or the local preimage of $x$ satisfies (72).

For target points such that $\operatorname{Re}z\notin(-1,1)$, we can skip the process of carefully selecting the length of the nearest subinterval, and proceed immediately with recursive bisection of $\{-1,1\}$. This completes the algorithm, which we list in its entirety in algorithm 1.

6 Numerical results

To test the robustness of our method across a range of $\alpha$ values, we solve the modified Helmholtz equation and the modified Stokes equations. Based on the discussion in section 4 the parameter ${C_{\varepsilon}}$ is set to confirm that prescribed error tolerances are satisfied.

6.1 The modified Helmholtz equation

We solve the modified Helmholtz equation (1) inside the annulus centered at the origin, defined by a circle of radius 0.3, and a circle of radius 0.6, with Dirichlet boundary conditions given by a fundamental solution (21) centered in the inner circle,

[TABLE]

The exact solution to this problem is equal to the expression for the Dirichlet boundary condition, evaluated in $\Omega$. See figure 7 for a visualization of (80) for $\alpha=10,\,100$. We see that $K_{0}(\alpha|x-x_{0}|)$ decays rapidly for $\alpha=100$.

To solve the above problem numerically, we represent the solution using the double layer potential

[TABLE]

where $G$ is the double layer kernel (54).

Enforcing the boundary condition (80) gives a second kind integral equation in $\sigma$,

[TABLE]

We solve this using the Nyström method, discretizing the boundary using $16$-point Gauss–Legendre panels, with $15$ panels on the inner circle, and $30$ panels on the outer. For the bound (76) we use use ${R_{\varepsilon}}=3.5$, and for (74) we set ${C_{\varepsilon}}=3.7$ by reading off figure 5 to achieve a tolerance of $10^{-14}$. Following the notation of (5), the kernel-split of (54) is given by (56) – (61).

For a large range of $\alpha$ 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). Quadrature by recursive bisection is applied both when solving for the density (83) and for evaluating the solution (82).

The results, shown in figure 8, demonstrate that our subdivision algorithm is capable of avoiding the catastrophic loss of accuracy otherwise present above a threshold $\alpha$, and that the additional cost incurred from it is proportional to $\log\alpha$. We do not satisfy the error tolerance of $10^{-14}$; around one digit of accuracy appears to be lost, presumably due to the additional interpolation steps involved.

We also try other values for ${C_{\varepsilon}}$ to confirm that we stay under a given tolerance. Based on figure 5 we have ${C_{\varepsilon}}=18.7,\,12.8,\,8.2$ for the tolerances $10^{-2},\,10^{-6},\,10^{-10}$. The largest corresponding relative error for each tolerance, taken over the plotted range of values of $\alpha$, are $3.7\cdot 10^{-5},\,6.7\cdot 10^{-8},\,1.8\cdot 10^{-11}$. Clearly the errors are at least one digit below the prescribed tolerance. It is not suprising that the set value for ${C_{\varepsilon}}$ gives a lower error than figure 5 suggests. In the subdivision algorithm introduced in this paper, the first panel in the recursive scheme always has as its center the projection of the target point on the boundary. Therefore target points are never close to panel edges, where errors from the kernel-split quadrature tend to be greater than for target points towards the panel’s center [13].

6.2 The modified Stokes equations

We solve the modified Stokes equations (3) in a setting analogous to (80); the geometry, and the discretization thereof, is the same, and the boundary data and the solution is given by the associated fundamental solution (106). As opposed to the modified Helmholtz equation, no error estimates are presented in this paper for the modified Stokes equation, which can be used to set ${C_{\varepsilon}}$. Instead, one can use the measured numerical errors for a flat panel for different values of $h$, each scaled with $1/h^{p}$, as suggested in section 4.3. Here $p$ is chosen such that the error curves collapse. For the modified Stokes in a double layer formulation, such results are shown in figure 6, with $p=1$.

Following the methodology presented in section 4.4, we set ${C_{\varepsilon}}$ to achieve certain tolerances. For the tolerances $10^{-2},\,10^{-6},\,10^{-10},\,10^{-14}$, we set ${C_{\varepsilon}}=16.7,\,11.5,\,7.4,\,4.4$ based on figure 6. The corresponding maximum errors over all alpha are $3.8\cdot 10^{-5},\,4.8\cdot 10^{-8},\,1.1\cdot 10^{-11},\,5.3\cdot 10^{-14}$, shown in figure 9. The magnitudes of the errors, for a given tolerance, is consistent with the results for the modified Helmholtz equation, presented in the previous section. The errors are well below the prescribed tolerance, meaning more subdivisions are applied than necessary. The hypothesis is the same as for the modified Helmholtz; for the subdivision algorithm the target point, or the projection thereof, on the boundary is always centered on the new panel. Thus its never close to the panel edges, where the error tends to be greater.

7 Conclusions

We present a robust recursive algorithm that allows the method of Helsing et al. to be applied, for any $\alpha$, to the modified Helmholtz equation, modified biharmonic equation and modified Stokes equations on smooth geometries. Before, this was not possible for large $\alpha$, which corresponds to small time-steps 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 $\log\alpha$. Our choice of the parameters ${C_{\varepsilon}}$ and ${R_{\varepsilon}}$ is based on numerical observations and provide excellent results.

Acknowledgments

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.)

Declarations

The authors declare no competing interests.

Appendix A Formulas for recursive computation of the Legendre-log

integrals

A.1 Result

The integrals

[TABLE]

can be computed recursively using the formulas

[TABLE]

The first two formulas, for $g_{0}(x)$ and $g_{1}(x)$, have finite limits as $x\to\pm 1$.

A.2 Proof

The formulas for $g_{0}$ and $g_{1}$ follow from direct integration of (84), with $P_{0}(\rho)=1$ and $P_{1}(\rho)=\rho$. For $n>1$, the following two results are useful:

[TABLE]

Insertion of the recursion formula (86) into (84) gives

[TABLE]

Integration by parts of remaining integral, using (87) and $((\rho-x)\log(\rho-x))^{\prime}=1+\log(\rho-x)$,

[TABLE]

We have $P_{n}(1)=1$ and $P_{n}(-1)=(-1)^{n}$, so it follows that $I_{1}=0$. Furthermore,

[TABLE]

Consequently, since $n\geq 0$,

[TABLE]

which leads to the recursion formulas (85).

Appendix B Kernel-splits

Here we complement the kernel-split for the modified Helmholtz equation with kernel-splits for the modified biharmonic equation and the modified Stokes equations.

B.1 Modified biharmonic equation

The single layer kernel for the modified biharmonic equation is

[TABLE]

Inserting (22) the explicit split becomes

[TABLE]

Here, $G^{L}$ contains the modified Bessel function $I_{0}$, like the single layer kernel for the modified Helmholtz equation. The double layer kernel is

[TABLE]

by (55). The resulting decomposition is

[TABLE]

B.2 The modified Stokes equations

We present only the double layer kernel for the modified Stokes equations, also known as the stresslet. In the Einstein summation convention it has the closed-form

[TABLE]

where $\delta_{ij}$ is the Kronecker delta and $\mathbf{r}=\mathbf{x}-\mathbf{y}$. Here, the functions $\mathcal{G}_{1}$–$\mathcal{G}_{3}$ are

[TABLE]

Since they are expressed in terms of the modified Bessel functions $K_{0}$ and $K_{1}$ we can use the decompositions (22) and (55) to write the expressions in explicit singularities. We have

[TABLE]

where

[TABLE]

With these definitions the kernel-split form of the stresslet is

[TABLE]

[TABLE]

where

[TABLE]