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

This paper presents a quadrature method for evaluating layer potentials in two dimensions close to periodic boundaries, discretized using the trapezoidal rule. It is an extension of the method of singularity swap quadrature, which recently was introduced for boundaries discretized using composite Gauss-Legendre quadrature. The original method builds on swapping the target singularity for its preimage in the complexified space of the curve parametrization, where the source panel is flat. This allows the integral to be efficiently evaluated using an interpolatory quadrature with a monomial basis. In this extension, we use the target preimage to swap the singularity to a point close to the unit circle. This allows us to evaluate the integral using an interpolatory quadrature with complex exponential basis functions. This is well-conditioned, and can be efficiently evaluated using the fast Fourier transform. The resulting method has exponential convergence, and can be used to accurately evaluate layer potentials close to the source geometry. We report experimental results on a simple test geometry, and provide a baseline Julia implementation that can be used for further experimentation.


Introduction
In integral equation methods, one of the long-standing challenges is the evaluation of layer potentials for targets points close to the source geometry.The integrals by which these layer potentials are computed are referred to as nearly singular, since the integrand is evaluated close to a singularity.This near singularity reduces the smoothness of the integrand, thereby severely reducing the accuracy of common quadrature rules such as Gauss-Legendre and the trapezoidal rule.To tackle this, many specialized quadrature methods have been proposed for the case of the domain being a one-dimensional line.For panel-based Gauss-Legendre quadrature in the plane, the high-order method of [7] (often referred to as the Helsing-Ojala method) is efficient and well-established.For the trapezoidal rule, which is exponentially convergent on a closed curve, several fixed-order methods have been proposed that reduce the impact of the near singularity by means of a local correction, see e.g.[5,6,8].In addition, exponentially convergenent near evaluation is available through the "globally compensated" quadrature method introduced in [7] and extended in [3].
Recently, the method of singularity swap quadrature (SSQ) was proposed in [1].The method evaluates nearly singular line quadratures in two and three dimensions by first "swapping" the singularity from the target point in the plane to the preimage of the target point in the curve's parameterization.After the swap, the integral can be evaluated using interpolatory quadrature with a monomial basis, following the method of [7].This was shown to improve the accuracy for curved panels compared to [7], and allowed the method to be extended to line integrals in three dimensions.
The SSQ method was in [1] described for curves discretized using panel-based Gauss-Legendre quadrature.In this brief follow-up, we show how the method can be extended to closed curves in the plane that have been discretized using the trapezoidal rule, resulting in a global correction with exponential convergence.In doing so, we assume that the reader has at least some familiarity with the original method.

Method
Since we are considering the problem in two dimensions, let us identify R 2 with C, and let Γ be a closed curve parameterized by a smooth function γ(t) = g 1 (t) + ig 2 (t), with t ∈ [0, 2π).Furthermore, let σ be a smooth function defined on Γ, which we will refer to as the density.For a point z ∈ C, typically lying close to Γ, our integrals of interest are here Being periodic, these integrals are suitable for discretization using an N-point trapezoidal rule, due to its exponential convergence [9].The quadrature nodes are then γ(t j ), j = 0, . . ., N − 1, with t j = 2πj/N and the quadrature weights w j = 2π/N .We will now outline how these discretized integrals can be evaluated using SSQ.Throughout, we will assume that we only have access to the discrete values of the functions σ, γ and γ at the nodes, and not to their analytic definitions.

The Cauchy integral
We begin with the simplest (and perhaps most common) case, the Cauchy integral, Let us now assume that z is close to Γ, which corresponds to the integrand in (3) having a singularity close to the real line at the preimage t * ∈ C such that 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 | Im t * | .In order to reduce the error for small 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 γ.We apply a Fast Fourier Transform (FFT) to the node values γ(t 0 ), . . ., γ(t N −1 ) to get the truncated Fourier series where K = N −1 2 , assuming N odd, and γk = F[γ](k) are the Fourier coefficients of γ.This has a one-time cost of O(N log N ), for a given discretization.Substituting this expansion into (4), we can find t * at an O(N ) cost using Newton's method (assuming O(1) iterations to convergence).A good initial guess for t * is the corresponding t-value of the grid point on Γ that is closest to z, We assume that Γ is parametrized in a counter-clockwise direction, so that Im t * > 0 implies z interior to Γ, 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 γ(t) − γ(t * ) −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, 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 , . . ., t N −1 , we use an FFT to compute its truncated Fourier series (at an O(N log N ) cost), such that 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 fk decay sufficiently fast.The integrals p k can be evaluated analytically (at an O(N ) cost) by rewriting them as contour integrals on the unit circle S. Let Then e ikt = ξ(t) k and We can evaluate this integral using the residue theorem.If |ζ| < 1 (or equivalently Im t * > 0, corresponding to z being an interior point), then the integrand has a simple pole at ζ, with residue In addition, if k ≤ 0, then the integrand has a pole of order 1 − k at the origin, with residue The two residues cancel for k ≤ 0 and |ζ| < 1 (i.e.Im t * > 0), so we get This completes the method, which has a total cost of 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 Im t * | .In constrast, the coefficients fk (t * ) decay rapidly and with a rate that stays nearly constant as z → Γ.
A weakness of the method is that in order to find t * , one must evaluate the analytic continuation of γ using a truncated Fourier series, which naturally is most accurate close to Γ.For z far away from Γ 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.

Remark 1
The reader may note that one could in fact apply a variant of the Helsing-Ojala quadrature directly to (3).First solve interpolation problem where τ i = γ(t i ) and σ i = γ(t i ).Then evaluate the layer potential as where the integrals can be exactly evaluated using residue calculus, just as is done above for the unit circle.This will in fact work perfectly fine for some geometries (the unit circle in particular), but the interpolation problem can be severely ill-conditioned in a way that is hard to control, since it depends on the geometry itself.This could potentially be investigated further, but that is deemed beyond the scope of this work.

Higher order integrals
For higher order integrals (m > 1), we have 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 We then need expressions for the integrals Evaluating the residues in the same way as for the m = 1 case, we get the general expression for integer m ≥ 1, (24)

Log kernel
Consider now the log integral, also known as the Laplace single-layer potential, for a real-valued density σ, where f is assumed to be a smooth function.Proceeding in a similar fashion as before, we find t * and separate the integral as 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| = Re log r, such that we can evaluate it using SSQ, Just as with p k , the integrals q k can be transformed into integrals on the unit circle, 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 Im t * .

Im t * < 0
Beginning with the case |ζ| > 1 (or Im t * < 0), we note that ξ − ζ = 0 on the entire unit disc.We can therefore choose the branch cut of the complex logarithm such that it never intersects the unit circle.This choice allows us to evaluate the integral in (29) as a contour integral on the unit circle, with a possible pole of order 1 − k at the origin, Since f0 ∈ R, we only need the real part of q k at k = 0. We then have that q is for Im t * < 0 computed as (31) We are now left with the case |ζ| < 1 (or Im t * > 0).Here, ξ − ζ = 0 at at point on the unit disk, so the integral of log(ξ − ζ) on the unit circle must intersect the branch cut in the ξ plane going from ζ to infinity.We choose the branch cut such that it passes through the point ξ = 1.We define 1 + and 1 − to be the limits approaching 1 from above and below in the plane, such that Then the integral on the unit circle S is, In order to evaluate this, we let C be the closed contour following the unit circle from 1 + to 1 − , then going from 1 − to ζ on the negative side of the branch cut, and then going from ζ to 1 + on the positive side of the branch cut, The integrand is analytic inside C, except for a pole of order 1 − k at the origin when k ≤ 0, with residue given by (30).The integrals along the branch cut are Combining the results, and just as before keeping only the real part of the k = 0 integral, we arrive at the final results for Im t * > 0,

Experiments
In order to demonstrate the method outlined above, we here report experimental results for the by now well-known "starfish" geometry γ(t) = (1 + 0.3 cos 5t)e it .The method has been implemented in Julia, and is available online as open source together with the scripts that produce the numerical in this manuscript 1 .

Laplace double layer potential
As a first demonstration, we reuse the Laplace test problem from [1].That is, we evaluate the solution to the Laplace equation ∆u = 0 in the domain bounded by Γ, with boundary condition u = u e (z) = log|3 + 3i − z|.The solution to this problem can be represented using the Laplace double layer potential (DLP) u(z) = Im I 1 (z), which leads to a second-kind integral equation in  the real-valued density σ.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 Γ 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 × 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 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 O(10 −13 ).

Convergence
For a more systematic investigation of the method, we create the following setup.Let Γ be the starfish domain from the previous section, and discretize it using N points.We then create 100 evaluation points as z = γ(t * ), where Re t * ∈ [0, 2π) and Im t * = d.In order to test points both interior and exterior to Γ, we use d = {±0.01,±0.02, ±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 σ and reference solution used are: • For I L we use the σ(t) = g 1 (t)g 2 (t) and compute the reference solution using adaptive Gauss-Kronrod quadrature.
• For I m and interior points (Im t * > 0) we use σ(τ ) = τ 3 + τ , with reference solution • For I m and exterior points (Im t * < 0) we use σ(τ ) = τ −1 , with reference solution 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.

Conclusions
We have extended the singularity swap quadrature (SSQ, [1]) method to closed 2D curves discretized using the trapezoidal rule.This extension builds on the simple observation that interpolatory quadrature can be used on a periodic integrand if the problem is first swapped to the unit circle.The method relies on representing quantities as Fourier expansions, and as a consquence the domintaing cost of the method is that of applying the fast Fourier transform (FFT).
The method is accurate, exponentially convergent, and relatively simple to implement (and we provide source code).The fact that the cost is O(N log N ) per target point makes the method unfeasible for large problems, but that is a consequence of the underlying trapezoidal rule being global (and exponentially convergent).For maximum efficiency, composite Gauss-Legendre quadrature with the original SSQ likely remains the better option.
In the original SSQ paper, it was shown that the method can be applied to line integrals in 3D by finding the singularity preimage through analytic continuation of the squared-distance function.The same technique can be applied to the present method, for closed curves in three dimensions.What remains is analytical evaluation of the integrals of the interpolatory quadrature.This is straightforward for the logarithmic kernel, but the integrals corresponding to power-law kernels are more challenging.Deriving expressions for these is currently work in progress.

Figure 1 :
Figure 1: Example results for the Cauchy integral on the starfish geometry and density from sec-3.2, for two different t * .The Fourier coefficients of the integrand in (3) decay approximately as ρ −|k| , where ρ = e | Im t0| , which gets very slow as z → Γ.On the other hand, the coefficients of the regularized function f in (9) decay with a rate that changes little as z → Γ.

Figure 2 :
Figure 2: Error in the Laplace double layer potential on a domain with N = 400 points on the boundary, comparison between the trapezoidal rule and SSQ.

Figure 3 :
Figure 3: Convergence of the trapezoidal rule and SSQ for the log kernel and m = 1, 2, 3, measured as the maximum error on a set of target points such that Im t * = d.Left column shows interior points, right column shows exterior points.