Singularity swap quadrature for nearly singular line integrals on closed curves in two dimensions

We begin with the simplest (and perhaps most common) case, the Cauchy integral,

$\displaystyle I_{1}(z)$

$\displaystyle=\oint_{\Gamma}\frac{\sigma\operatorname{d\!}\tau}{\tau-z}=\int_{0}^{2\pi}\frac{\sigma(t)\gamma^{\prime}(t)\operatorname{d\!}t}{\gamma(t)-z}.$

(3)

Let us now assume that $z$ is close to $\Gamma$, which corresponds to the integrand in (3) having a singularity close to the real line at the preimage $t^{*}\in\mathbb{C}$ such that

$\displaystyle\gamma(t^{*})=z.$

(4)

In [4] and [2], it is shown that this significantly reduces the convergence rate of the trapezoidal rule, with the error being proportional to $e^{-N|\operatorname{Im}{t^{*}}|}$. In order to reduce the error for small $\operatorname{Im}t^{*}$, we will now show how to evaluate $I_{1}(z)$ using SSQ.

Step one is to find the preimage $t^{*}$, using Newton’s method applied to a Fourier expansion of $\gamma$. We apply a Fast Fourier Transform (FFT) to the node values $\gamma(t_{0}),\dots,\gamma(t_{N-1})$ to get the truncated Fourier series

	$\displaystyle\gamma(t)$	$\displaystyle\approx\sum_{k=-K}^{K}\hat{\gamma}_{k}e^{ikt},$		(5)
	$\displaystyle\gamma^{\prime}(t)$	$\displaystyle\approx\sum_{k=-K}^{K}ik\hat{\gamma}_{k}e^{ikt},$		(6)

where $K=\frac{N-1}{2}$, assuming $N$ odd, and $\hat{\gamma}_{k}=\mathcal{F}[\gamma](k)$ are the Fourier coefficients of $\gamma$. This has a one-time cost of $\mathcal{O}(N\log N)$, for a given discretization. Substituting this expansion into (4), we can find $t^{*}$ at an $\mathcal{O}(N)$ cost using Newton’s method (assuming $\mathcal{O}(1)$ iterations to convergence). A good initial guess for $t^{*}$ is the corresponding $t$-value of the grid point on $\Gamma$ that is closest to $z$,

$\displaystyle t^{*}\approx\underset{t_{j}}{\operatorname{arg\,min}}\mathinner{\!\left\lvert\gamma(t_{j})-z\right\rvert}.$

(7)

We assume that $\Gamma$ is parametrized in a counter-clockwise direction, so that $\operatorname{Im}t^{*}>0$ implies $z$ interior to $\Gamma$, and vice versa.

Step two is to swap out the singularity using a function that regularizes the integrand, and leads to an integrable singularity. Contrary to the Gauss-Legendre panel case, we can not regularize the denominator $\left(\gamma(t)-\gamma(t^{*})\right)^{-1}$ using the factor $(t-t^{*})$, as the latter is not periodic. Instead, we use $(e^{it}-e^{it^{*}})$, which corresponds to swapping the singularity to a point close to the unit circle (this will be evident shortly). Then,

$\displaystyle I_{1}(z)$

$\displaystyle=\int_{0}^{2\pi}\underbrace{\frac{\sigma(t)\gamma^{\prime}(t)\left(e^{it}-e^{it^{*}}\right)}{\gamma(t)-\gamma(t^{*})}}_{f(t,t^{*})}\frac{\operatorname{d\!}t}{e^{it}-e^{it^{*}}}.$

(8)

The function $f$ is now the original integrand with the singularity at $t^{*}$ removed, and we expect it to be smooth. Since its values are known at the equidistant nodes $t_{0},\dots,t_{N-1}$, we use an FFT to compute its truncated Fourier series (at an $\mathcal{O}(N\log N)$ cost), such that

$\displaystyle I_{1}(z)$

$\displaystyle\approx\sum_{k=-K}^{K}\hat{f}_{k}(t^{*})\underbrace{\int_{0}^{2\pi}\frac{e^{ikt}\operatorname{d\!}t}{e^{it}-e^{it^{*}}}}_{p_{k}(t^{*})}=\bm{\hat{f}}\cdot\bm{p}.$

(9)

This can be interpreted as our integral being approximated by an interpolatory quadrature on the unit circle, which we can evaluate to high accuracy as long as the coefficients $\hat{f}_{k}$ decay sufficiently fast. The integrals $p_{k}$ can be evaluated analytically (at an $\mathcal{O}(N)$ cost) by rewriting them as contour integrals on the unit circle $S$. Let

	$\displaystyle\xi(t)$	$\displaystyle=e^{it},$		(10)
	$\displaystyle\xi^{\prime}(t)$	$\displaystyle=ie^{it},$		(11)
	$\displaystyle\zeta$	$\displaystyle=e^{it^{*}}.$		(12)

Then $e^{ikt}=\xi(t)^{k}$ and

	$\displaystyle p_{k}(t^{*})$	$\displaystyle=\frac{1}{i}\int_{0}^{2\pi}\frac{e^{i(k-1)t}ie^{it}\operatorname{d\!}t}{e^{it}-e^{it^{*}}}$		(13)
		$\displaystyle=\frac{1}{i}\int_{0}^{2\pi}\frac{\xi(t)^{k-1}\xi^{\prime}(t)\operatorname{d\!}t}{\xi(t)-\zeta}$		(14)
		$\displaystyle=\frac{1}{i}\oint_{S}\frac{\xi^{k-1}\operatorname{d\!}\xi}{\xi-\zeta}.$		(15)

We can evaluate this integral using the residue theorem. If $|\zeta|<1$ (or equivalently $\operatorname{Im}t^{*}>0$, corresponding to $z$ being an interior point), then the integrand has a simple pole at $\zeta$, with residue

$\displaystyle\operatorname{Res}\left[\frac{\xi^{k-1}}{\xi-\zeta},\zeta\right]=\zeta^{k-1},\quad\mbox{ if }|\zeta|<1.$

(16)

In addition, if $k\leq 0$, then the integrand has a pole of order $1-k$ at the origin, with residue

$\displaystyle\operatorname{Res}\left[\frac{\xi^{k-1}}{\xi-\zeta},0\right]=-\zeta^{k-1},\quad\mbox{ if }k\leq 0.$

(17)

The two residues cancel for $k\leq 0$ and $|\zeta|<1$ (i.e. $\operatorname{Im}t^{*}>0$), so we get

$\displaystyle p_{k}(t^{*})=2\pi e^{i(k-1)t^{*}}\begin{cases}1,&\text{if }\operatorname{Im}t^{*}>0\text{ and }k\geq 1,\\ -1,&\text{if }\operatorname{Im}t^{*}<0\text{ and }k\leq 0,\\ 0,&\text{otherwise}.\end{cases}$

(18)

This completes the method, which has a total cost of $\mathcal{O}(N+N\log N)$ per target point.

The method can perhaps be understood in terms of Fourier coefficents. The accuracy of the trapezoidal rule depends on the regularity of the integrand, which in term is reflected in the decay of the Fourier coefficients of the integrand. As seen in fig. 1, these appear to decay as $e^{-|k\operatorname{Im}t^{*}|}$. In constrast, the coefficients $\hat{f}_{k}(t^{*})$ decay rapidly and with a rate that stays nearly constant as $z\to\Gamma$.

A weakness of the method is that in order to find $t^{*}$, one must evaluate the analytic continuation of $\gamma$ using a truncated Fourier series, which naturally is most accurate close to $\Gamma$. For $z$ far away from $\Gamma$ the rootfinding process will struggle, and the accuracy of SSQ can deteriorate. However, this happens at ranges where the base trapezoidal rule is sufficiently accurate.

For higher order integrals ($m>1$), we have

$\displaystyle I_{m}(z)$

$\displaystyle=\oint_{\Gamma}\frac{\sigma\operatorname{d\!}\tau}{\left(\tau-z\right)^{m}}=\int_{0}^{2\pi}\frac{\sigma(t)\gamma^{\prime}(t)\operatorname{d\!}t}{(\gamma(t)-z)^{m}},\quad m\in\mathbb{Z}^{+}.$

(21)

Evaluating this integral with SSQ is completely analogous to the case when $m=1$, with the difference that we instead regularize with $(e^{it}-e^{it^{*}})^{m}$, such that

$\displaystyle f(t,t^{*})=\sigma(t)\gamma^{\prime}(t)\left(\frac{e^{it}-e^{it^{*}}}{\gamma(t)-\gamma(t^{*})}\right)^{m}.$

(22)

We then need expressions for the integrals

$\displaystyle p_{k}^{m}(t^{*})=\frac{1}{i}\oint_{S}\frac{\xi^{k-1}\operatorname{d\!}\xi}{\left(\xi-\zeta\right)^{m}}.$

(23)

Evaluating the residues in the same way as for the $m=1$ case, we get the general expression for integer $m\geq 1$,

$\displaystyle p_{k}(t^{*})=2\pi\frac{\prod_{j=1}^{m-1}(k-j)}{(m-1)!}e^{i(k-m)t^{*}}\begin{cases}1,&\text{if }\operatorname{Im}t^{*}>0\text{ and }k\geq m,\\ -1,&\text{if }\operatorname{Im}t^{*}<0\text{ and }k\leq 0,\\ 0,&\text{otherwise}.\end{cases}$

(24)

Consider now the log integral, also known as the Laplace single-layer potential, for a real-valued density $\sigma$,

	$\displaystyle I_{L}(z)$	$\displaystyle=\oint_{\Gamma}\sigma\log|\tau-z|\>|\operatorname{d\!}\tau|$		(25)
		$\displaystyle=\int_{0}^{2\pi}\underbrace{\sigma(t)|\gamma^{\prime}|}_{f(t)}\log|\gamma(t)-z|\operatorname{d\!}t,$		(26)

where $f$ is assumed to be a smooth function. Proceeding in a similar fashion as before, we find $t^{*}$ and separate the integral as

$\displaystyle I_{L}(z)=\int_{0}^{2\pi}f(t)\left(\log|\gamma(t)-z|-\log|e^{it}-e^{it^{*}}|\right)\operatorname{d\!}t+\int_{0}^{2\pi}f(t)\log|e^{it}-e^{it^{*}}|\operatorname{d\!}t.$

(27)

The first integral is now regularized, and can be evaluated using the trapezoidal rule. The second integral we rewrite using a truncated Fourier series and the relation $\log|r|=\operatorname{Re}\log r$, such that we can evaluate it using SSQ,

$\displaystyle\int_{0}^{2\pi}f(t)\log|e^{it}-e^{it^{*}}|\operatorname{d\!}t\approx\operatorname{Re}\sum_{k=-N/2}^{N/2}\hat{f}_{k}\underbrace{\int_{0}^{2\pi}e^{ikt}\log(e^{ikt}-e^{ikt^{*}})\operatorname{d\!}t}_{q_{k}(t^{*})}=\bm{\hat{f}}\cdot\bm{q}.$

(28)

Just as with $p_{k}$, the integrals $q_{k}$ can be transformed into integrals on the unit circle,

$\displaystyle q_{k}(t^{*})=\frac{1}{i}\int_{S}\xi^{k-1}\log(\xi-\zeta)\operatorname{d\!}\xi.$

(29)

Note that this is not necessarily a closed contour integral, as there is a branch cut in the complex logarithm that needs to be taken into account. We will now show to evaluate $q_{k}(t^{*})$ depending on the sign of $\operatorname{Im}t^{*}$.

As a first demonstration, we reuse the Laplace test problem from [1]. That is, we evaluate the solution to the Laplace equation $\Delta u=0$ in the domain bounded by $\Gamma$, with boundary condition $u=u_{e}(z)=\log\mathinner{\!\left\lvert 3+3i-z\right\rvert}$. The solution to this problem can be represented using the Laplace double layer potential (DLP) $u(z)=\operatorname{Im}I_{1}(z)$, which leads to a second-kind integral equation in the real-valued density $\sigma$. Once this integral equation is solved, the layer potential $u(z)$ can be evaluated in the domain using quadrature, and compared to the exact solution $u_{e}(z)$.

In fig. 2 we show results for the Laplace DLP when $\Gamma$ is discretized using $N=400$ points that are equidistant in the parametrization $t$. The layer potential is evaluated using the trapezoidal rule on a $400\times 400$ grid inside the domain, but for points close to the boundary we also evaluate using SSQ and compare the results. As expected, evaluation using the trapezoidal rule leads to $\mathcal{O}(1)$ errors near the boundary, while SSQ is accurate at all points where it is used. However, full machine precision accuracy is not recovered; the maximum error in SSQ close to the boundary is $\mathcal{O}(10^{-13})$.

For a more systematic investigation of the method, we create the following setup. Let $\Gamma$ be the starfish domain from the previous section, and discretize it using $N$ points. We then create 100 evaluation points as $z=\gamma(t^{*})$, where $\operatorname{Re}t^{*}\in[0,2\pi)$ and $\operatorname{Im}t^{*}=d$. In order to test points both interior and exterior to $\Gamma$, we use $d=\{\pm 0.01,\pm 0.02,\pm 0.04\}$. For a range of $N$, we then compute the integrals $I_{L},I_{1},I_{2},I_{3}$ using both trapezoidal and singularity swap quadrature, and for each $N$ save the maximum error across all $z$. The density $\sigma$ and reference solution used are:

• •

  For $I_{L}$ we use the $\sigma(t)=g_{1}(t)g_{2}(t)$ and compute the reference solution using adaptive Gauss-Kronrod quadrature.

• •

  For $I_{m}$ and interior points ($\operatorname{Im}t^{*}>0$) we use $\sigma(\tau)=\tau^{3}+\tau$, with reference solution $I_{m}(z)=2\pi i\sigma^{(m)}(z)/(m-1)!$.

• •

  For $I_{m}$ and exterior points ($\operatorname{Im}t^{*}<0$) we use $\sigma(\tau)=\tau^{-1}$, with reference solution $I_{m}(z)=2\pi iz^{-m}$.

The results of this investigation are shown in fig. 3. As expected, the convergence rate of the trapezoidal rule depends strongly on the distance to the boundary, which is not the case for SSQ. There appears to be a weak dependence on the rate for interior points, but that could be an effect of the chosen test problem. More concerning is that lowest attainable error appears to increase for smaller $d$ when $m>1$. It is unclear what causes this effect, although in experiments there appears to be an upper bound to the error growth.

The key contributions of this work were made with support from the Knut and Alice Wallenberg Foundation, under grant no. 2016.0410.