Singularity swapping method for nearly singular integrals based on trapezoidal rule

In this subsection, we are mainly concerned with the integral

where $f$ is an analytic and periodic function on $[0,2\pi]$. Rewrite $f(t)=u(z)$ with $z=e^{it}$. The original integral (4) becomes

Throughout this paper, a complex contour such as $|z|=1$ is always to be understood as traversed once in the counterclockwise direction. Since $f$ is analytic on $[0,2\pi]$, by analytic extension, there exists $r>1$ such that $u$ is analytic in the region $r^{-1}\leq|z|\leq r$. For any positive integer $n$, the approximation to $I(u)$ using trapezoidal rule with $n$ quadrature nodes is simply given by

To analyze the error of $T_{n}(u)-I(u)$, note that the function $\frac{1}{z^{n}-1}\frac{1}{iz}$ has $n$ simple poles at $e^{2\pi ik/n}$, $k=1,2,\cdots,n$, with residues equal to $1/(in)$. By the residue theorem, it holds

Since $u$ is analytic in $r^{-1}\leq|z|\leq r$, we can change $|z|=1$ to the circle $|z|=r^{-1}$ in equation (5) without changing the value. Thus, combining (5) and (7) gives the representation of $T_{n}(u)-I(u)$ as an integral,

Using (8) and the boundedness of $u(z)$, we obtain the exponential convergence of trapezoidal rule.

Based on (8), one also finds that $T_{n}(u)=I(u)$ for $u(z)=z^{-n+1},\dots,z^{n-1}$ and $z^{n}-z^{-n}$. Therefore, the trapezoidal rule can be treated as a Gaussian quadrature on the unit circle $|z|=1$[25].

We now turn to the evaluation of nearly singular integral defined on $|z|=1$, which is given in the form of

where $u$ is analytic in the annulus $R^{-1}\leqslant|z|\leqslant R$ for some $R>1$ and $K$ is analytic in a smaller annulus $r^{-1}\leqslant|z|\leqslant r$ with $1<r<R$.

Presumably, $r$ is very close to $1$, which means $K$ is badly behaved on the unit circle. According to (9), trapezoidal rule gives very poor convergence rate that is on the order of $O(r^{-n})$. To overcome this difficulty, one way is to adaptively place very dense quadrature nodes on the unit circle that are close to the singularity of $K$. Another way is to numerically recalculate the weights of trapezoidal rule based on the variance of Euler-Macularin formula so that the singularity of $K$ is taken into account[2]. Both ways will only lead to an algebraic convergence of (10). Here we propose a modified trapezoidal rule, in which the quadrature weights can be explicitly obtained and the error still keeps exponential decay. The main idea is to replace $K$ with its Laurent expansion $K_{n}$ and then apply the trapezoidal rule to $K_{n}u$.

More precisely, let

be the $n$th order Laurent expansion of $K$ at zero, where

Note that the $n$-th term of the Laurent series (11) is halved, as well as the $-n$-th term, which is not mandatory but we do this for symmetry and consistency with the classical trigonometric interpolation, as shown later. The modified trapezoidal rule for integral (10) is simply

Now let us analyze the error of formula (13). It follows from equation (7) that

Write $K=K_{n}^{-}+K_{n}+K_{n}^{+}$, where

Using the fact that $K_{n}^{-}u$ is analytic in $r^{-1}\leqslant|z|\leqslant R$ and $K_{n}u+K_{n}^{+}u$ is analytic in $R^{-1}\leqslant|z|\leqslant r$, it holds

From the Cauchy’s estimate

we have the following results on $|z|=R$

Same estimates hold for $\left|\frac{K_{n}(z)}{z^{-2n}-1}\right|$ and $|K_{n}^{+}(z)|$ on $|z|=R^{-1}$. Combining (2.2), (17) and (18) gives

Furthermore, equation (2.2) implies $T_{2n}(K_{n}u)=I(Ku)$ for $u(z)=z^{-n+1},\dots,z^{n-1}$ and $z^{n}+z^{-n}$.

To summarize, we have the following result.

Note that the presented method can be regarded as interpolating $u$ at $e^{2\pi ik/(2n)}$ into the form of $a_{-n+1}z^{-n+1}+\cdots+a_{n-1}z^{n-1}+a_{n}(z^{n}+z^{-n})$ and then integrating analytically, so it can be treated as an interpolatory quadrature. Compared to the original trapezoidal rule, it only has half of the polynomial order of accuracy, as we have preassigned $2n$ quadrature nodes at $e^{2\pi ik/(2n)}$. It is possible to construct a generalized Gaussian quadrature that exactly integrates $2n$ functions using only $n$ function evaluations [13], but that would require solving nonlinear equations and the corresponding quadrature nodes will vary with $K$ [6]. From that perspective, the presented method of fixed quadrature nodes is more convenient for practical computation.

To illustrate the application of the modified trapezoidal rule (13), we give the explicit expressions for the nearly singular integrals of some typical kernel functions, which can be used to evaluate the near field of single and double layer potentials for Laplace’s and Helmholtz equations.

When $\Im t_{0}=0$, formula (22) reduces to the well known Kress quadrature [17]:

Note that the vector $\{R_{n}(t_{0},j\pi/n)\}_{j=0}^{2n-1}$ can be calculated by the Fast Fourier Transform (FFT). In the next section, we will reduce the Laplace’s and Helmholtz single layer potential to the evaluation of (21).

Such a kernel can be appeared in the gradient of Laplace’s and Helmholtz single layer potentials. When $|z_{0}|=1$, the left-hand side exists as Cauchy principal value, and a limit is needed to take on the right-hand side when $z_{0}=e^{\pi ik/n}$ for some $k=1,\cdots,2n$. One can also make use of the formula (24) to construct the quadrature rule for kernel (3) with $m\geq 1$ through differentiation of $K(z)=\frac{1}{z-z_{0}}$.

Corollary (2.5) is essentially the same as Corollary 2.4. Later in the paper, we will use it to compute the near field of double layer potential for Laplace’s and Helmholtz equations.

In this part, we mainly focus on the second issue and discuss the kernels given in the form of (2) and (3), more specifically, the single and double layer potentials for Laplace’s equations, due to their wide applications in the potential theory. We expect the same idea can be extended to other kernels in a straightforward way. Identify $\mathbb{R}^{2}$ with $\mathbb{C}$ and assume that $\Gamma$ is parameterized by $g:\mathbb{R}\to\mathbb{C}$ that maps the interval $[0,2\pi)$ to $\Gamma$. Thus the parametric form of the single layer potential with logarithmic kernel (2) is

where $f$ is an analytic function that consists of the density function and the Jacobian $|g^{\prime}(t)|$. Bringing $x$ near $\Gamma$ introduces a singularity $t_{0}$ near real axis in the analytic continuation of $\log|x-g(t)|^{2}$. Specifically, it has a singularity when the complexified squared distance

goes to zero, which gives $x_{1}-g_{1}(t)=\pm i(x_{2}-g_{2}(t))$. Taking the negative case, the singularity $t_{0}$ thus obeys

and one may check by Schwarz reflection that the positive case gives its conjugate $\bar{t}_{0}$. Equation (29) can be rewritten as $x=g(t_{0})$. The method to efficiently find $t_{0}$ will be discussed in the next subsection 3.2. For now we assume that $g:\mathbb{C}\to\mathbb{C}$ is 2$\pi$-periodic, analytic and single-valued in a sufficiently large strip $\{t\in\mathbb{C}:|\Im t|<\log r\}$, so the only solution to $x=g(t_{0})$ is $t_{0}=g^{-1}(x)$, which is called the preimage of $x$ under $g$.

Now we rewrite the integral (27) as

The first integral in (30) is analytic in a larger domain $\{t\in\mathbb{C}:|\Im t|<\log r\}$, so the classic trapezoidal rule (6) can be used. The second integral can be numerically evaluated by the modified trapezoidal rule (22) in Corollary 2.3. Hence, the singularities $t_{0}$ and $\overline{t_{0}}$ from $\log|x-g(t)|^{2}$ is swapped to the function $\log\left(4\sin\frac{t-t_{0}}{2}\sin\frac{t-\overline{t_{0}}}{2}\right)$ through (30).

For the Laplace double layer potential, it is given by

where $n(g(t))$ is unit outward normal vector of $\Gamma$ at $g(t)$. One way to swap its singularity is given by [1]

The second term $\frac{(1-e^{-i(t-t_{0})})g^{\prime}(t)}{x-g(t)}f(t)$ in (32) is smooth and plays the same role as $f(t)$ in (26). We can therefore follow the quadrature rule given by (26) to discretize (32). However, there is another approach to swap the singularity:

which we have found is numerically more accurate than (32).

One issue for the double layer evaluation is that when the target point $x$ is close to one of the sources $g(j\pi/n)$ with $0\leqslant j\leqslant 2n-1$, the accuracy is lost due to the singularity at $t=j\pi/n$. Similar phenomenon has also been observed in [1]. To address this issue, we apply the singularity subtraction technique by rewriting (31) into

where

and $D$ is the bounded domain with boundary $\Gamma$. Then swap the singularity as in (33). The new integrand vanishes at $t=j\pi/n$, so the rounding errors can be ignored. There is no such simple subtraction techniques for single layer potentials [7]. However, it is also not necessary since the single layer potential suffers much less from cancellation errors, as one would see in the numerical experiments.

As a comment on the extension of the techniques mentioned above to many other standard linear PDEs of interest (Helmholtz, Maxwell, Navier, etc), the integrals that we need to evaluate may be written as a combination of (27) and (31). For example, for the Hemholtz single layer potential, the kernel $K(x,y)$ can be split by

and for the Helmholtz double layer potential, the kernel $K(x,y)$ can be split by

The first terms in both (36) and (37) have logarithmic type singularities and the second term in (37) is the kernel of Laplace double layer potential. Therefore one may apply the same technique as Laplace to evaluate them.

Now the only missing ingredient in the scheme is how to find the singularity $t_{0}=g^{-1}(x)$ efficiently. Newton iteration is definitely an available option. However, it requires the evaluation of function values in the complex plane and the efficiency highly depends on the initial guess. Here we propose a quadrature method to find an approximation of $t_{0}$, which does not need any initial guess and in many situations is good enough to be directly applied in the singularity swapping method.

For the convenience of derivation, let $z_{0}=e^{it_{0}}$ and $g(t)=\gamma(e^{it})$, so that the equation $g(t_{0})=x$ becomes $\gamma(z_{0})=x$. Suppose $u$ is analytic in the annulus $r^{-1}\leqslant|z|\leqslant r$ except a simple pole $z_{0}$ with residue equal to $W$, where $z_{0}$ and $W$ are unknown. Similar to the analysis of equation (7), it holds

and

From the two formulas above, we can conclude that

which provides a method to find the pole $z_{0}$ with high accuracy. The advantage of this formula is only function values on the unit circle are used and does not need any initial guess. To find an approximate solution to equation $\gamma(z_{0})=x$, we can apply (40) to $u(z)=\frac{\gamma^{\prime}(z)}{\gamma(z)-x}$, which gives

In order to numerically test the accuracy of $\hat{z}_{0}$, we compare the performance of three different quadrature methods in the evaluation of nearly singular integrals. The first method is the classical trapezoidal rule given by (6), the second method is the interpolation quadrature given by (24) with exact $z_{0}$, and the third method is also the interpolation quadrature (24) but with $z_{0}$ replaced by $\hat{z}_{0}$. Figure 1 tests the convergence of three integrals that are defined on $|z|=1$ and have the same nearly singular kernel $(z-1.1)^{-1}$. On top of the kernel, the three integrands are $z^{-10}+z^{10}$, $e^{z}$ and $(z-2)^{-1}$, respectively. For the first function $z^{-10}+z^{10}$, our analysis indicates that the interpolation quadrature (24) with $z_{0}=1.1$ is exact for $n\geqslant 20$, which is clearly shown in Figure 1(a). Using $\hat{z}_{0}$ only loses one algebraic accuracy. The second function $e^{z}$ is entire in $\mathbb{C}$. The interpolation quadrature achieves super-geometric convergence. Using $\hat{z}_{0}$ only causes the convergence speed to slow down slightly, as shown in 1(b). The third function $(z-2)^{-1}$ is analytic in a neighborhood of $|z|=1$ but not throughout the complex plane. By Theorem 2.2, the exact interpolation quadrature has a geometric convergence that is on the order of $2.2^{-n/2}$, which is confirmed by Figure 1(c). Here using $\hat{z}_{0}$ causes the convergence error to grow by a factor of $n$, which is still geometrically convergent. On the other hand, application of the classic trapezoidal rule converges very slow throughout all these three integrals. The results clearly show the interpolation quadrature can significantly improve the convergence rate, either using exact $z_{0}$ or the approximate value $\hat{z}_{0}$.

To summarize, numerical observation indicates replacing $z_{0}$ with $\hat{z}_{0}$ into the quadrature formula (24) gives a new residue term

which is slightly larger than the approximation error analyzed in (20). We briefly analyze this term by taking $|z_{0}|>1$ as an example:

The last identity holds because $T_{2n}\left(\frac{1+z/z_{0}}{2}z^{n}u\right)=O(|z_{0}|^{-n})$. In a word, the factor $n$ in (42) comes from the $n$-th power in the error term of (20).

If such an error is unacceptable, usually one Newton iteration

is enough to eliminate the error caused by singularity approximation. In order to improve the efficiency further, we can decrease the number of nodes in (41) and increase the number of Newton iterations. It is numerically found that for the same accuracy, one Newton iteration is roughly equivalent to reducing the number of nodes by half. By balancing the operations of the two parts, we can obtain a good approximation $\hat{z}_{0}$ within almost $O(n)$ operations.