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

The numerical method of Helsing and co-workers evaluates Laplace and related layer potentials generated by a panel (composite) quadrature on a curve, efficiently and with high-order accuracy for arbitrarily close targets. Since it exploits complex analysis, its use has been restricted to two dimensions (2D). We first explain its loss of accuracy as panels become curved, using a classical complex approximation result of Walsh that can be interpreted as “electrostatic shielding” of a Schwarz singularity. We then introduce a variant that swaps the target singularity for one at its complexified parameter preimage; in the latter space the panel is flat, hence the convergence rate can be much higher. The preimage is found robustly by Newton iteration. This idea also enables, for the first time, a near-singular quadrature for potentials generated by smooth curves in 3D, building on recurrences of Tornberg–Gustavsson. We apply this to accurate evaluation of the Stokes flow near to a curved filament in the slender body approximation. Our 3D method is several times more efficient (both in terms of kernel evaluations, and in speed in a C implementation) than the only existing alternative, namely, adaptive integration.


Introduction
Integral equation methods enable efficient numerical solutions to piecewise-constant coefficient elliptic boundary value problems, by representing the solution as a layer potential generated by a so-called "density" function defined on a lower-dimensional source geometry [4,31]. This work is concerned with accurate evaluation of such layer potentials close to their source, when this source is a one-dimensional curve embedded in 2D or 3D space. In 2D, this is a common occurrence: the curve is the boundary of the computational domain or an interface between different materials. In 3D, such line integrals represent the fluid velocity in non-local slender-body theory (SBT) for filaments in a viscous (Stokes) flow [17,26,28], and also may represent the solution fields in the electrostatic [42] or elecromagnetic [11] response of thin wire conductors. Applications in the Stokes case include simulation of straight [45], flexible [35,46] or closed-loop [33] fiber suspensions in fluids. SBT has recently been placed on a more rigorous footing as the small-radius asymptotic limit of a surface boundary integral formulation [30,33,34]. In this work we focus on non-oscillatory kernels arising from Laplace (electrostatic) and Stokes applications, although we expect that by singularity splitting (e.g. [20]) the methods we present could be adapted for oscillatory or Yukawa kernels.
For numerical solutions, the density function must first be discretized [31,44]. A variety of methods are used in practice. If the geometry is smooth enough, and the interacting objects are far enough away that the density is also smooth, then a nonadaptive density representation is sufficient. Such representations include low-order splines on uniform grids [46], global spectral discretization by the periodic trapezoid rule for closed curves [6,18,31], and, for open curves, global Chebychev [11,35] or Legendre [45] expansions. However, if the density has regions where refinement is needed (such as corners or close interactions with other bodies), a composite ("panel") scheme is better, as is common with Galerkin boundary-element methods [44]. On each panel a high-order representation (such as the Lagrange basis for a set of Gauss-Legendre nodes) is used; panels may then be split in either a graded [24,36,44] or adaptive [37,41,50] fashion.
Our goal is accurate and efficient evaluation of the potential due to a density given on a panel, at arbitrarily close target points. Not only is this crucial for accurate solution evaluation once the density is known, but it is also key to constructing matrix elements for a Nyström or Galerkin solution [31,44] in the case when curves are close to each other [6,23,36]. Specifically, let the panel Γ be an open smooth curve in R d , where d = 2 or d = 3. We seek to evaluate the layer potential at the target point x / ∈ Γ . Here the density f is a smooth function defined on Γ , and K is a kernel that has a singularity as x → y. In two dimensions the dominant singularity can be either logarithmic, or of power-law form where m is generally even. Here k denotes smooth (possibly tensor-valued) functions on R d . In 3D only the latter form arises, for m generally odd. We will assume that Γ is described by the parametrization g : R → R d that maps the standard interval [−1, 1] to Γ , Thus the parametric form of the desired integral (1) is wheref (t) := f (g(t)) is the pullback of the density to parameter space. The difficulty of evaluation of a layer potential with given density is essentially controlled by the distance of the target point x from the curve. In particular, if x is far from Γ , the kernel K (x, y) varies, as a function of y, no more rapidly than the curve geometry itself. Thus a fixed quadrature rule, once it integrates accurately the product of density f and any geometry functions (such as the "speed" g ), is accurate for all such targets. Here, "far" can be quantified relative to the local quadrature node spacing h and desired accuracy: e.g. in the 2D periodic trapezoid rule setting the distance should be around (log 10)h/2π ≈ 0.37 h times the desired number of correct digits (a slight extension of [7, Remark 2.6]). However, as x approaches Γ , the kernel K (x, y) becomes increasingly non-smooth in the region where y is near x, and any fixed quadrature scheme loses accuracy. This is often referred to as the layer potential integral being nearly singular.
There are several approaches to handling the nearly singular case. The most crude is to increase the number of nodes used to discretize the density; this has the obvious disadvantage of growing the linear system size beyond that needed to represent the density, wasting computation time and memory. Slightly less naive is to preserve the original discretization nodes, but interpolate from these nodes onto a refined quadrature scheme (with higher order or subdivided panels [46]) used only for potential evaluations. However, to handle arbitrarily close targets one has to refine the panels in an adaptive fashion, afresh for each target point. As a target approaches Γ an increasing level of adaptive refinement is demanded. The large number of conditionals and interpolations needed in such a scheme can make it slow.
This has led to more sophisticated "special" rules which do not suffer as the target approaches the curve. They exploit properties of the PDE, specifically, in 2D, the relation of harmonic to analytic functions. In the 2D high-order setting, several important contributions have been made in the last decade, including panel-based kernelsplit quadrature [19][20][21]23,36], globally compensated spectral quadrature [6,23], quadrature by expansion (QBX) [7,29,48], harmonic density interpolation [39], and asymptotic expansion [12]. In 3D, little work has been done on high-order accurate evaluation of nearly singular line integrals. Many 3D applications are in the context of active or passive filaments in viscous fluids, but current quadrature methods are limited to either using analytical formulae for straight segments [45], or regularizing the integrand using an artificial length scale [13,25,46].
In this paper, we propose a new panel-based quadrature method for nearly singular line integrals in both 2D and 3D. The method incurs no loss of efficiency for targets arbitrarily close to the curve. (Note that, in the 2D case it is usual for one-sided limits on Γ to exist, but in 3D limits on Γ do not generally exist.) Our method builds on the panel-based monomial recurrences introduced by Helsing-Ojala [23] and extended by Helsing and co-workers [19,20,22,36] [21, Sec. 6], with a key difference that, instead of interpolating in the physical plane (which is associated with the complex plane), we interpolate in the (real) parametrization variable of the panel. Rather than exploiting Cauchy's theorem in the physical plane, the singularity is "swapped" by cancellation for one in the parameter plane; this requires a nearby root-finding search for the (analytic continuation of the) function describing the distance between the panel and the target. It builds on recent methods for computing quadrature error estimates in two dimensions [3]. We include a robust and efficient method for such root searches.
This new formulation brings two major advantages: 1. The rather severe loss of accuracy of the Helsing-Ojala method for curved panelswhich we quantify for the first time via classical complex approximation theory-is much ameliorated. As we demonstrate, more highly curved panels may be handled with the same number of nodes, leading to higher efficiency. 2. The method extends to arbitrary smooth curves in 3D. Remarkably, the search for a root in the complex plane is unchanged from the 2D case.
As a side benefit, our method automatically generates the information required to determine whether it is needed, or whether direct quadrature using the existing density nodes is adequate, using results in [2,3]. The structure of this paper is as follows. Section 2 is concerned with the 2D problem; we give the necessary background and introduce our new "singularity swap" method. Section 3 then extends the method to deal with problems in 3D. Section 4 summarizes the entire proposed algorithm, unifying the 2D and 3D cases. Section 5 presents further numerical tests of accuracy in 2D, and accuracy and efficiency in a 3D Stokes slender body application, with a comparison against standard adaptive quadrature. We conclude in Sect. 6.

Two dimensions
This section gives the necessary background and introduces our method for problems in 2D. Section 2.1 reviews the convergence rate for direct evaluation of the potential, i.e. using a direct quadrature scheme. This also serves to introduce complex notation and motivate special quadratures. Section 2.2 summarizes the complex monomial interpolatory quadrature of Helsing-Ojala [19,23] for evaluation of a given density close to its 2D source panel. In Sect. 2.3 we explain and quantify the loss of accuracy in the aforementioned method due to panel bending. Finally in Sect. 2.4 we use the "singularity swap" to make a real monomial version that is less sensitive to bending, and which will form the basis of the 3D method.

Summary of the direct evaluation error near a panel
Here we review known results about the approximation of the integral (5) directly using a fixed high-order quadrature rule. In the analytic panel case this is best understood by analytic continuation in the parameter. Let t j be the nodes and w j the weights for an n-node Gauss-Legendre scheme on [−1, 1], at which we assume the density is known. 1 Substituting the rule into (5) gives the direct approximation where the nodes are y j = g(t j ), the density samples f j = f ( y j ), and the modified weights vanishes exponentially with n, with a rate which grows with the size of this neighborhood. Specifically, define E ρ , the Bernstein ρ-ellipse, as the image of the disc |z| < ρ under the Joukowski map (For example, Fig. 1a shows this ellipse for ρ ≈ 1.397.) Then, if F is analytic and bounded in E ρ , there is a constant C depending only on ρ and F such that An explicit bound [47,Thm. 19.3] Returning to (5), the integrand of interest is F(t) = K (x, g(t))f (t)|g (t)|, for which we seek the largest ρ such that F is analytic in E ρ . Bringing x near Γ induces a singularity near [−1, 1] in the analytic continuation in t of K (x, g(t)), a claim justified as follows. Each of the kernels K (x, g(t)) under consideration has a singularity when the squared distance goes to zero, which gives x 1 − g 1 (t) = ± i(x 2 − g 2 (t)). Taking the negative case, the singularity (denoted by t 0 ) thus obeys and one may check by Schwarz reflection that the positive case gives its conjugate t 0 . This suggests identifying R 2 with C and introducing Thus, in this complex notation, the Eq. (11) for t 0 is written For now we assume that γ : C → C is analytic in a sufficiently large complex neighborhood of [−1, 1], so the nearest singularity is t 0 = γ −1 (ζ ), the preimage of ζ under γ ; see upper and lower panels of Fig. 1a. Thus our claim is proved. Inverting (8) gives the elliptical radius ρ of the Bernstein ellipse on which any t ∈ C lies, with sign chosen such that ρ > 1.
We then have, for any of the kernels under study, that (9) holds for any ρ < ρ(γ −1 (ζ )). Thus images of the Bernstein ellipses, γ (E ρ ), control contours of equal exponential convergence rate, and thus (assuming that the prefactor C does not vary much with the target point), at fixed n, also the contours of equal error magnitude. A more detailed analysis for the Laplace double-layer case in [2] [3, Sec. 3.1] improves 2 the rate to ρ = ρ(γ −1 (ζ )), predicts the constant, and verifies that the error contours very closely follow this prediction. Indeed, the lower panel of our Fig. 1a illustrates, for an analytic choice of density f , that the contour of constant error (modulo oscillatory "fingering" [3]) well matches the curve γ (E ρ ) passing through the target ζ . The top-left plots in Figs. 2 and 3 show the shrinkage of these error contours (images of Bernstein ellipses) as n, the number of nodes on the panel, grows.
For kernels other than the Laplace double-layer, error results are similar [2,29]. This quantified breakdown in accuracy as the target approaches Γ motivates the special quadratures to which we now turn. ) their error in potential evaluation using the direct scheme; an example target ζ and its preimage t 0 = γ −1 (ζ ) (large green dots); Bernstein ellipse for t 0 and its image (dotted red); Schwarz singularity τ * = i/4k and its preimage t * = i/2k (red * 's); and Cartesian grid lines (grey) in the t-plane and their images (white). b Equipotential curves of g(z) for the same panel Γ , showing "shielding" of the Schwarz singularity (red * ) thus small g ≈ 0.2 and ρ ≈ 1.22. Contrast ρ ≈ 2.14 for the Schwarz preimage in the t-plane. c Convergence in n of polynomial approximation error of f over Γ (solid lines) for several curvatures k, compared to Cρ −n (dotted) predicted by Sect. 2.3 (color figure online)

Interpolatory quadrature using a complex monomial basis
We now review the special quadrature of Helsing and coworkers [19,23,36] that approximates (1) in the 2D case. In the complex notation of (14), given a smooth density defined on Γ , for all standard PDEs of interest (Laplace, Helmholtz, Stokes, etc), the integrals that we need to evaluate can be written as the contour integrals "Appendix A" reviews how to rewrite common boundary integral kernels as such contour integrals; note that f may include other smooth factors than merely the density. Other 2nd-order elliptic 2D PDE kernels may be split into terms of the above form, as pioneered by Helsing and others in the Helmholtz [20], axisymmetric Helmholtz [22], elastostatics [19,Sec. 9], Stokes [37], and (modified) biharmonic [21,Sec. 6] cases.
An innovation of the method was to interpolate the density in the complex coordinate τ rather than the parameter t. Thus, This is most numerically stable if for now we assume that Γ has been transformed by rotation and scaling to connect the points ± 1. Here the (complex) images of the quadrature nodes are {τ j } n j=1 := {γ (t j )} n j=1 . The monomial coefficients c k can be found by collocation at these nodes, as follows. Let c be the column vector with entries {c k }, f be the vector of values { f (τ j )}, and A be the Vandermonde matrix with entries A i j = τ j−1 i , for i, j = 1, . . . , n. Then enforcing (19) at the nodes results in the n × n linear system Ac = f (20) which can be solved for the vector c. Note that working in this monomial basis, even in the case of τ j on the real axis, is successful even up to quite high orders, although this is rarely stated in textbooks.
Remark 1 (Backward stability) It is known for {τ j } being any set of real or complex nodes that is not close to the roots of unity that the condition number of the Vandermonde matrix A in (20) must grow at least exponentially in n [16,38]. However, as pointed out in [23, App. A], extreme ill-conditioning in itself introduces no loss of interpolation accuracy, at least for n < 50. Specifically, as long as a vector c has small relative residual Ac− f / f , then the resulting monomial is close to the data at the nodes. Assuming weak growth in the Lebesgue constant with n, the interpolant is thus uniformly accurate. Such a small residual norm, if it exists for the given right-hand side, is found by a solver that is backward stable. We recommend either partially pivoted Gaussian elimination (e.g. as implemented by mldivide in MATLAB) which is O(n 3 ), or the Björck-Pereyra algorithm [10] which is O(n 2 ).
As with any interpolatory quadrature, once the coefficients {c k } have been computed, the integrals (17) and (18) can be approximated as where the values {q k } and { p m k } are the exact integrals of the monomial densities {τ k−1 } multiplied with the logarithmic and power-law kernels respectively, The benefit of this (somewhat unusual) complex monomial basis is that these integrals can be efficiently evaluated through recursion formulae, as follows.

Exact evaluation of p m k and q k by recurrence
Recall that Γ has been transformed by rotation and scaling to connect the points ±1.
We first review the Cauchy kernel m = 1 case [23, Sec. 5.1]. The contour integral for k = 1 exploits independence of the path Γ : Here N ζ ∈ Z is a winding number that, for the standard branch cut of the logarithm, is zero if ζ is outside the domain enclosed by the oriented closed curve comprising Γ traversed forwards plus [−1, 1] traversed backwards, while N ζ = +1 (−1) when ζ is inside a part of this domain that is enclosed counterclockwise (clockwise) [23,Sec. 7]. The following 2-term recurrence is derived by adding and subtracting ζ τ k−1 from the numerator of the formula (24): One can then recur upwards to get p 1 m for all k. For m > 1 we get analogously, We emphasize that these simple exact path-independent formulae are enabled by the choice of the complex monomial basis in the coordinate τ .
Remark 2 (Transforming for general panel endpoints) We now review the changes that are introduced to the answers I L and I m when rescaling and rotating a general panel Γ parameterized by γ (t), t ∈ [−1, 1], into the required one which connects ± 1, and similarly transforming the target. Define the complex scale factor s 0 := (γ (1) − γ (−1))/2 and origin τ 0 = (γ (1) + γ (−1))/2. Then the affine map from a given location τ to the transformed locationτ , given bỹ can be applied to all points on the given panel and to the target to give a transformed panelΓ and targetζ . From this, applying the formulae above one can evaluateĨ L and I m . By inserting the change of variables one then gets the desired potentials due to the density on Γ at target ζ , Here a new integral is needed (the total "charge"), easily approximated by

The adjoint method for weights applying to arbitrary densities
The above procedure computes I L and I m of (17)- (18) for a specific function f (such as the density) on Γ , using its samples f on the nodes to solve first for a coefficient vector c. In practice it is more useful instead to compute the weight vectors λ L := {λ L j } n j=1 and λ m := {λ m j } n j=1 which, for any vector f , approximate I L and I m when their inner products are taken against f . That is, indicates non-conjugate transpose. One application is for filling near-diagonal blocks of the Nyström matrix arising in a boundary integral formulation: (λ L ) T and (λ m ) T form row blocks of the matrix.
We review a result presented in [23, Eq. (51)] (where a faster numerical method was also given in the m = 1 case). Let λ be either λ L or λ m , let I be the corresponding integral (17) or (18), and let p be the corresponding vector {q k } or { p m k } as computed in Sect. 2.2.1. Writing our requirement for λ, we insert (20) to get where the last form simply expresses (21) or (22). But this last equality must hold for all c ∈ C n , thus which is a linear system whose solution is the weight vector λ. Notice that this is an adjoint method [27] to compute λ.

The effect of panel curvature on Helsing-Ojala convergence rate
The Helsing-Ojala method reviewed above is easy to implement, computationally cheap, and can achieve close to full numerical precision for target points very close to Γ . However, there is a severe loss of accuracy as Γ becomes curved. This is striking in the central column of plots in Fig. 2. This observation led Helsing-Ojala [23] and later users to recommend upsampling from a baseline n = 16 to a larger n, such as n = 32, even though the density and geometry were already accurately resolved at n = 16. The "HO" k = 0.4 plot in Fig. 3 shows that even using n = 32 has loss of accuracy near a moderately curved panel. This forces one to split panels beyond what is needed for the density representation, just to achieve good evaluation accuracy, increasing the cost. The error in this method is essentially entirely due to the error of the polynomial approximation (19) of f on Γ , because each monomial is integrated exactly (up to rounding error) as in Sect. 2.2.1. For this error there is a classical result that the best polynomial approximation converges geometrically in the degree, with rate controlled by the largest equipotential curve of Γ inside of which f can be analytically continued. Namely, given the set Γ , let g be the unique solution to the potential problem in R 2 (which we identify with C), For each ρ > 1, define 3 the level curve C ρ := {τ ∈ C : g(τ ) = log ρ}. An example g and its equipotential curves are shown in Fig. 1b; note that they are very different from the Bernstein ellipse images in the lower plot of Fig. 1a.
Let Γ ⊂ C be a set whose complement (including the point at infinity) is simply connected. Let f : Γ → C extend analytically to a function analytic inside and on C ρ , for some ρ > 1. Then there is a sequence of polynomials p n (z) of degree n = 0, 1, . . . such that

where C is a constant independent of n and z.
To apply this to (19), we need to know the largest region in which f may be continued as an analytic function, which is in general difficult. In practical settings the panels will have been refined enough so that on each panel preimagef is interpolated in the parameter t ∈ [−1, 1] to high accuracy, thus we may assume that any singularities off are distant from [−1, 1]. Yet, in moving from the t to τ plane, the shape of the panel Γ itself can introduce new singularities in f which will turn out to explain well the observed loss of accuracy, as follows.

Proposition 4 (geometry-induced singularity in analytic continuation of the density) Let the pullbackf be a generic function analytic at t
Proof Since the derivative vanishes, the Taylor expansion of the parameterization at t * takes the form which must match, term by term, the Taylor expansion of the pullbackf (t) =f 0 + f 1 (t − t * ) + . . . . But genericallyf 1 = 0, in which case there is a contradiction. Similar (more elaborate) arguments are known for singularities in the extension of Helmholtz solutions [32]. Here τ * is an example of a singularity of the Schwarz function [14,43] for the curve Γ . Recall that the Schwarz function G is defined by G(τ ) = γ (γ −1 (τ )), which is the analytic reflection of τ through the arc Γ . At points where γ = 0, the inverse map γ −1 , hence the Schwarz function, has a square-root type branch singularity. Figure 1a shows that for a moderately bent parabolic panel, the Schwarz singularity is quite close to the concave side. Note that the above argument relies onf having a generic expansion, which we believe is typical; we show below that its convergence rate prediction holds well in practice.
In summary, smooth but curved panels often have a nearby Schwarz singularity, which generically induces a singularity at this same location in the analytic continuation of the density f ; then the equipotential for Γ at this location controls the best rate of polynomial approximation, placing a fundamental limit on the Helsing-Ojala convergence rate in n. Figure 1b shows that, for a moderately bent panel, the close singularity is effectively shielded electrostatically by the concave panel, thus the convergence rate ρ = e g(τ ) is very small (close to 1). In Fig. 1c we test this quantitatively: for a parabolic panel and simple analytic pullback densityf we compare the maximum error in n-term polynomial approximation on Γ with the prediction from the above Schwarz-Walsh analysis. The agreement in rate is excellent over a range of curvatures. This also matches the loss of digits seen in the central column of Figs. 2 and 3: e.g. for even the mild curvature k = 0.25, n = 16 achieves only 5 digits, and n = 32 only 10 digits.
Conjecture 5 (Schwarz singularities at focal points) Note that in the parabola example of Fig. 1, a square-root type Schwarz singularity appears near Γ 's highest curvature point, at distance R/2 on the concave side, where R is the radius of curvature. This is also approximately true for ellipses, whose foci induce singularities [43, §9.4]. We also always observe this numerically for other analytic curves. Hence we conjecture that it is a general result for analytic curves, in an approximate sense that becomes exact as R → 0.

Interpolatory quadrature using a real monomial basis
We now propose a "singularity swap" method which pushes the interpolation problem to the panel preimage t ∈ [−1, 1], thus bypassing the rather severe effects of curvature just explained. Recalling that the target point is ζ , we first write the integrals (17) and (18) in parametric form where we have introduced the displacement function Q and the smooth function h, Now let t 0 be the root of Q nearest to [−1, 1]. (Unless the panel is very curved, there is only one such nearby root and it is simple.) Then Q(t)/(t −t 0 ) has a removable singularity at t 0 , and hence it and its reciprocal are analytic in a larger neighborhood of [−1, 1] than the above integrands. By multiplying and dividing by (t − t 0 ), which, as we detail below, can be handled as special cases of (17) and (18) for a new flat panel Γ = [−1, 1] written in the t plane. We have "swapped" the singularity from the τ to the t plane. As Fig. 1 illustrates, the preimage of the Schwarz singularity has a much higher conformal distance from [−1, 1] than the actual singularity in the τ -plane has from Γ , indicating much more rapid convergence.
Thus we now apply the Helsing-Ojala methods of Sect. 2.2, but to [−1, 1] in the t plane. Namely, as before, let t j and w j , j = 1, . . . , n be the Gauss-Legendre quadrature for [−1, 1]. The first integral in (37) is smooth, so direct quadrature is used. The second integral in (37) is evaluated via (21) and (29), setting Γ = [−1, 1], with the coefficients c k for the function h(t) found by solving the Vandermonde system as in Sect. 2.2. The first term (38) is smooth and plays the role of f (τ ) in (18), so its coefficients are found similarly. Equation (38) is then approximated by (22) with (25)- (28).

Remark 6
Since in this scheme the new flat "panel" [−1, 1] is fixed, one may LU factorize the Vandermonde matrix once and for all, then use forward and back-substitution to solve for the coefficients {c k } given the n samples of each of the two smooth functions,

Remark 7
With a very curved panel it is possible that more than one root of Q(t) is relevant; see Fig. 4. In such a case where poles z 1 and z 2 need to be cancelled, we have derived recursion formulae for p m k (z 1 , z 2 ) and q k (z 1 , z 2 ) which generalize those in Sect. 2.2.1. However, we have not found these to be needed in practice, so omit them.

Finding the roots of Q(t)
The only missing ingredient in the scheme just described is to find the nearest complex root t 0 satisfying Fig. 4 Illustration of the relationship between the roots of the displacement Q(t), to the left, and target points in 2D physical space, to the right. The point ξ corresponds to only one nearby root, w 1 , while the point ζ corresponds to two nearby roots, z 1 and z 2 . Note that only z 1 lies inside the Bernstein ellipse of radius ρ , marked green, such that it requires singularity swapping. The point γ (t * ) is the Schwarz singularity of the curve, marking the point where γ −1 stops being single-valued i.e., the preimage of the (complex) target point ζ under the complexification of the panel map γ ; see Fig. 4. We base our method on that of [3] (where roots were needed only for the purposes of error estimation). Let P n [γ ](t) denote a degree n − 1 polynomial approximant to γ (t) on [−1, 1]. Using the data at the Legendre nodes {t j }, forming P n [γ ](t) is both well-conditioned and stable if we use a Legendre (or Chebyshev) expansion, and has an O(n 2 ) cost. See [3] for details on how to use a Legendre expansion. Continuing each basis function to C gives a complex approximatioñ for which we seek roots nearest [−1, 1]. Given a nearby starting guess, a single root ofQ can be found using Newton's method, which converges rapidly and costs O(n) per iteration. A good starting guess is t init ≈ s(ζ ), recalling that s defined by Eq. (30) maps the panel endpoints to ± 1. With this initialization, Newton iterations converge to the nearest root in most practical cases. The iterations can, however, be sensitive to the initial guess when two roots are both relatively close to [−1, 1], and sometimes converge to the second-nearest root. As illustrated in Fig. 4, this typically happens only when the target point is on the concave side of a very curved panel.
A more robust, but also more expensive, way of finding the nearest root is to first find all n − 1 roots ofQ at the same time, and then pick the nearest one. This can be done using a matrix-based method, which finds the roots as the eigenvalues to a generalized companion (or "comrade") matrix [8]. This has a cost that is O(n 3 ) if done naively, and O(n 2 ) if done using methods that exploit the matrix structure [5].
In order to determine whether Newton iterations or a matrix-based method should be used for finding the root t 0 , we suggest a criterion based on the distance to the nearest Schwarz singularity preimage of the panel, t * , as this is a measure of the size of the region where γ −1 is single-valued. Let ρ be a Bernstein radius beyond which direct Gauss-Legendre quadrature is expected to have relative error . Then, ignoring the prefactor in (9) and solving for ρ, is a useful estimate. Recalling (16), then ρ(t * ) ≤ Cρ , where we choose the constant C = 1.1, indicates that there may be more than one nearby root, and that it is safer to use a matrix-based root finding method. For each panel, t * can be found at the time of discretization, by applying Newton's method to the polynomial P n [γ ](t) = 0, with the initial guess t * ≈ 0.

Remark 9
The root finding process can be made efficient by computing the interpolants P n [γ ](t) and P n [γ ](t) for each segment at the time of discretization, such that they can be reused for every target point ζ .

Tests of improved convergence rate for curved panels
We now compare the three methods: (i) direct Gauss-Legendre quadrature, (ii) Helsing-Ojala complex interpolatory quadrature, and (iii) the proposed real singularity swap quadrature. We evaluate the Laplace double layer potential from a single panel, The panel is the same parabola as in Fig. 1, i.e. g(t) = (t, kt 2 ), or in complex form, γ (t) = t + ikt 2 , and we explore the dependence on curvature k. In complex form, which is necessary to apply the quadratures, the layer potential is (see "Appendix A") Figures 2 and 3 compare the three schemes for various curvatures, for n = 16 and n = 32 respectively. For n = 32 all data is upsampled by Lagrange interpolation from n = 16 (which is already adequate to represent it to machine precision). The error is in all cases measured against a reference solution computed using adaptive quadrature (integral in MATLAB with abstol and reltol set to eps), and the relative error E rel computed against the maximum value of u on the grid. To see how they behave when pushed to their limits, we deliberately apply the quadratures much further away than necessary, i.e., also in the region where direct quadrature achieves full numerical accuracy.
The direct errors (left column of each figure) have magnitudes controlled by the images of Bernstein ellipses under the parameterization, as explained in Sect. 2.1. The Helsing-Ojala scheme (middle column) is much improved over direct quadrature, but, because of its convergence rates dependence (Sect. 2.3) the reduction in accuracy is severe as curvature increases, although it can to some can extent be ameliorated by using n = 32. In addition, there is an error in the quadrature that grows with the distance from the panel, particularly for n = 32, due to amplification of roundoff error in the upward recurrence for the integrals { p 1 , . . . , p n }. This was also noted in [23], where the suggested fix was to instead run the recurrence backwards for distant points, starting from a value of p n computed using the direct quadrature. This avoids the problem of catastrophic accuracy loss, but only because specialized quadrature is used where direct quadrature would be sufficient. The bottom middle plot of Fig. 3 shows at best 10 accurate digits, and that only over quite a narrow range of distances, making the use of the HO algorithm to give a uniform accuracy for all target locations rather brittle.
For the proposed singularity swap quadrature (right columns), there is still a loss of accuracy with increased curvature, but the improvement over Helsing-Ojala is striking. The instability of upwards recurrence is also much more mild.

Remark 10
The loss of accuracy with increased curvature can for the singularity swap quadrature be explained by considering the roots of the displacement function Q(t). We have "swapped" out the nearest root t 0 , so the region of analyticity of the regularized integrand h(t)(t − t 0 )/Q(t) is bounded by the second nearest root of Q(t), illustrated by z 2 in Fig. 4. This limits the convergence rate of the polynomial coefficients {c k }, such that we may need to upsample the panel in order for the coefficients to fully decay. In the case of our parabolic panel, the second nearest root gets closer with increased curvature, which explains why upsampling to 32 points is necessary to achieve full accuracy in some locations. Note that for the upsampling to be beneficial, the coefficients for k > n must be nonzero. This is only the case if the components of the integrand (ρ, γ , γ ) are upsampled separately, before the integrand is evaluated at the new nodes.
where Γ = g ([−1, 1]) is an open curve in R 3 , andf incorporates the density and possibly other smooth denominator factors in the kernel K . Fixing the target x, then introducing the 3D squared distance function R 2 [analogous to (10)] and the smooth function h, we can write (44) as This integrand has singularities at the roots of R(t) 2 , which, since the function is real for real t, come in complex conjugate pairs {t 0 , t 0 }. How to find these roots is discussed shortly in Sect. 3.2; for now we assume that they are known. As in 2D, if Γ is not very curved and the target is nearby there is only one nearby pair of roots.
We construct a quadrature for (47) as in the 2D case (38), except that now there is a conjugate pair of near singularities to cancel. Thus we write (47) as H (t) can be expected to be analytic in a much larger neighborhood of [−1, 1] than the integrand. As in Sect. 2.4, we now represent H (t) in a real monomial basis, which one can evaluate to machine accuracy using recurrence relations as described in the next section. Combining (50) and (49) gives the final quadrature approximation to (44), The adjoint method of Sect. 2.2.2 can again be used to solve for a weight vector λ m whose inner product with h approximates I m . As in Sect. 2.4, the error in (51) is essentially due to the convergence rate of the best polynomial representation (49) on [−1, 1], which is likely to be rapid because of the absence of curvature effects (Sect. 2.3).

Recurrence relations for 3D kernels on a straight line
The above singularity swap has transformed the integral on a curved line Γ in 3D to (48), whose denominator corresponds to that of a straight line in 3D. Such upward recursions are available in [45, App. B], where they were used for Stokes quadratures near straight segments. We will here present them with some improvements which increase their stability for t 0 in certain regions of C. We present the case of a single root pair {t 0 , t 0 } = t r ± it i , for orders m = 1, 3, 5, noting that, while we have derived formulae for double root pairs, we have found that they are never needed in practice.
To simplify notation (matching notation for I and for the distances to the parameter endpoints we write Beginning with m = 1 and k = 1, the integral (50) is This expression suffers from cancellation for t 0 close to [−1, 1], and must be carefully evaluated. To begin with, it is more accurate in the left half plane, even though the integral is symmetric in t r , so the accuracy can be increased by evaluating it after the substitution t r → −|t r |. In addition, the argument of the second logarithm of (54) suffers from cancellation when |t r | < 1 and t 2 i (1 − |t r |) 2 , even after the substitution. In this case we evaluate it using a Taylor series in t i , We achieve sufficient accuracy by applying this series evaluation to points inside the rhombus described by 4|t i | < 1−|t r |, evaluating the series to n = 11 terms. For k > 1 the upward recursions of [45] are stable, such that For m = 3, the formula from [45] for k = 1 contains a conditional statement for points on the real axis (d = 0), where the pair of two conjugate roots t 0 , t 0 merge into a double root. In finite precision, this property causes a region around the real axis where the formula is inaccurate, namely two cones extending outwards from the endpoints ±1. To get high accuracy there, we consider the integral with shifted limits, where the antiderivative S 3 is evaluated by forming a Maclaurin series of the integrand in t i , and then integrating that series exactly in s, thus We find that we get sufficiently high accuracy in the region of interest by truncating this series to 30 terms, denoted byS 3 (s), and finally evaluating the integrals for m = 3 as The series evaluation can be optimized by defining a succession of narrower cones, using fewer terms in each cone. For m = 5 and k = 1 we need to repeat the process of finding a power series that we can use in cones around the real axis, extending from the endpoints. We now consider and the antiderivative of the Maclaurin series of the integrand is Truncating this series to a maximum of 50 terms, denoted byS 5 (s), we now evaluate the integrals for m = 5 as Note that our expression for P 5 2 is different from the one found in [45]; ours avoids a conditional.

Finding the roots of R(t) 2
In three dimensions, the nearby roots of the squared distance function (45) do not correspond to a single point in space, as was the case in 2D (see Fig. 4). Instead, each complex conjugate pair of roots {t 0 , t 0 } near [−1, 1] corresponds to a circle in space looped around Γ . To see this, let g r and g i be the real and imaginary parts of the Then, the roots of (45) satisfy For this to hold in both real and imaginary components for a given t, the point x must satisfy The above equations describe the intersection of a plane with a sphere, and are satisfied by the circle of radius g i (t) that is centered at g r (t), and lies in the plane normal to g i (t). All points in R 3 corresponding to the root t lie on this circle. For an illustration of this, see Fig. 5.
In order to construct a quadrature for a given target point x, we need to find the roots t of (45), which are the points satisfying the complex equation We form a polynomial approximation P n [g](t) of g(t) using the existing real Gauss-Legendre nodes {g(t j )}. We then need to find the roots of the polynomial These can be found using either of the root-finding methods discussed in Sect. 2.4.1. We here limit ourselves to finding the root pair {t 0 , t 0 } closest to [−1, 1], using Newton's method. We find that a stable initial guess can be obtained using a linear mapping in the plane containing x and the two closest discretization points on Γ , which we denote by g j := g(t j ) and g k := g(t k ). The initial guess t init is then set to one of the two points satisfying This is the exact root for the case of Γ a straight line. Compared to the 2D case, we find that the accuracy of our method in 3D is more sensitive to the implementation details of the root-finding algorithm, especially for roots very close to the interval [−1, 1]. Below are a number of noteworthy observations: -The convergence of Newton's method deteriorates for roots that are close to the real axis (or on it, when | Re t 0 | > 1). This is because the root t 0 and its conjugate t 0 then lie very close to each other, and in the case of Im t 0 = 0 form a double root. To ameliorate this problem, we let our root-finding routine switch to Muller's method [40] if Newton has not converged within a certain number of iterations (we use 20). In our tests, this appears to be a stable root-finding process for all points requiring special quadrature, such that we can apply direct quadrature to points where neither Newton's nor Muller's methods converge. For additional robustness, one could also switch to the comrade matrix method in such cases. -Just as in the 2D case, once the root t 0 is known we can determine if special quadrature is needed using the Bernstein radius ρ(t 0 ) criterion at the end of Sect. 2.4.1. -We find that it is important for accuracy to find the roots using the approximation (74) of R 2 . If instead R 2 is approximated as P n [R 2 ](t), then 1-2 digits of accuracy are lost close to [−1, 1]. If roots are computed from this form, for example if the comrade matrix method is used, then they can be polished by taking one or two Newton steps of (74), thereby recovering the lost accuracy. -For panels discretized with more than n = 16 points, we find that the best accuracy is achieved if the roots are found using Legendre expansions truncated to 16 terms. In general, our recommendation is to not use an underlying discretization with n larger than 16, and then upsample to more points when evaluating the SSQ, if necessary (see Remark 10).

Summary of algorithm
Below is a step by step summary of our algorithm for evaluating the potential generated by a single open curve panel Γ , in 2D or 3D, discretized using n Gauss-Legendre points. Specifically, we describe how to get the weight vector λ as defined in Sect. 2.2.2, from which one of the line integrals (17), (18) or (44) can be approximated by I = λ T f for any sample vector f of densities (or density times smooth geometric factors).
Recall that potentials are defined by such line integrals, possibly in the 2D case via complexification as shown in "Appendix A". (We omit the optional finding of the panel's Schwarz singularity, and use of the robust comrade-matrix method given in Sect. 2.4.1, since we find that it is not needed in practice.) There is an overall algorithm option (flag) which may take the values "no upsampling", "upsampling" or "upsampling with upsampled direct". Either variety of upsampling is more expensive, but can increase accuracy; the second variety avoids some special weight computations so is the cheaper of the two. Precomputations for panel Γ : 1. Choose a suitable polynomial basis (e.g. Legendre or Chebyshev; we use the former). The analytic continuation of the chosen basis functions must be simple to evaluate (as is the case for Legendre). Form the degree n − 1 polynomial approximation of the panel parameterization P n [g](t) in this basis; this can be done by solving a n × n linear system with right-hand side the coordinates of the nodes y j = g(t j ), j = 1, . . . , n. 2. Given the desired tolerance , set the critical Bernstein radius ρ via (41).

For each given target point x:
1. Use a cheap criterion to check if x is a candidate for near evaluation, based on the minimum distance between the target point and the panel nodes: x candidate if min where D is a multiple of the panel length h (at n = 16 we use D = h for full accuracy). If not a candidate, use the direct rule (6) then exit; more specifically the weights are 2. Find the complex target preimage t 0 from the basis approximation to g(t): In 2D, use Newton's method applied to (40) to find the root t 0 . In 3D, use Newton (complemented by Muller's method) applied to (74) to find the root pair {t 0 , t 0 }. , and the right-hand side is the p just computed. We find that the Björck-Pereyra algorithm [10] is even faster than using a precomputed LU decomposition of A. 6. Finally, one must apply corrections to turn λ into a set of weights that act on the samples of f , rather than the singularity-swapped smooth functions involving h(t) or H (t). First multiply each weight λ j by the Jacobian factor |g (t j )|, according to (36) or (46). Then: -For the log case in 2D, according to (37), add w j log(Q(t j )/(t j − t 0 )) to each weight λ j . -For the power-law case in 2D, according to (38), multiply each weight by -In 3D, according to (48), multiply each weight by (t j −t 0 )(t j − t 0 )/R(t j ) 2 m/2 .

Numerical tests
In Sect. 2 we already performed basic comparisons of prior methods and the proposed singularity swap quadrature in the 2D case. Now we give results of more extensive 2D tests and the 3D case. The 2D and 3D quadrature methods of this paper have been implemented in a set of MATLAB routines. These are available online [1], together with the programs used for all the numerical experiments reported here. The computationally intensive steps of the 3D quadrature have been implemented in C and interfaced to MATLAB using MEX. The code is called using a list of target points, and has rudimentary parallelization in the form of OpenMP worksharing, so that each thread works on a subset of the target points. This allows very rapid computation of quadrature weights: on the computer used for the numerical results of this section, which is a 3.6 GHz quad-core Intel i7-7700, we can find the roots (preimages) of 4.7 × 10 6 target points per second, and compute 3.2 × 10 6 sets of target-specific quadrature weights per second, for n = 16 (without upsampling).

Evaluating a layer potential in two dimensions
To study quadrature performance when evaluating an actual layer potential, we solve a simple test problem using an integral equation method, then evaluate the solution using either the Helsing-Ojala method or our proposed method. Following [23] we test the Laplace equation Δu = 0 in the interior of the starfish-shaped domain Ω given by γ (t) = (1 + 0.3 cos(5t))e it , t ∈ [0, 2π), and we choose Dirichlet boundary conditions u = u e on ∂Ω with data u e (ζ ) = log 3 + 3i − ζ . We have picked this data to be very smooth, since the errors then are due to quadrature, rather than density resolution, making it easier to compare methods.
We discretize ∂Ω using composite Gauss-Legendre quadrature, with n = 16 points per panel. The panels are formed using the following adaptive scheme. The interval [0, 2π) is recursively bisected into segments that map to panels, until all panels satisfy a resolution condition, and all pairs of neighboring segments differ in length by a ratio this is no more than 2 (i.e. level-restriction). The resolution condition is based on the While somewhat ad-hoc, this produces a solution with a relative error that is comparable to , for a smooth problem (see Figs. 6 and 7).
The solution to the Laplace equation is represented using the double layer potential (43), which results in a second kind integral equation in ρ. This is solved using the Nyström method and the above discretization, which is straightforward due to the layer potential being smooth on ∂Ω. For further details on this solution procedure, see, for example, [20] or [23]. Once we have ρ, we evaluate the layer potential using quadrature at sets of points in Ω, and compute the error by comparing with the exact solution, u e . We compare two methods: 1. The Helsing-Ojala (HO) method summarized in Sect. 2.2. 2. The proposed singularity swap quadrature (SSQ), as outlined in Sect. 4.
In our first test, shown in Fig. 6, we solve the problem using a quite coarse discretization (adaptive with = 10 −6 ), choose ρ = 1.8 (corresponding via (41) to a more conservative error of 10 −8 ), and evaluate the solution on a 300×300 uniform grid covering the domain. Far from the boundary the solution shows around the expected 6 digit accuracy. Notice that, since only eight panels are used, they are highly curved. HO gets only 2 digits atñ = 16 and scarcely more atñ = 32, whereas SSQ gets 3 digits atñ = 16 and the "full" 6 digits atñ = 32. This shows that upsampled SSQ achieves the full expected accuracy using the efficient discretization chosen for this accuracy, whereas HO would require further panel subdivision beyond that needed to achieve the requested Nyström accuracy.
As a second test, we solve the same problem, this time using a fine discretization (adaptive with = 10 −14 ), and pick ρ = 3, corresponding via (41) to around ≈ 10 −15 at n = 16. The result is 32 panels (similar in number to the 35 panels used in [23]). First we evaluate the solution on a uniform 300 × 300 grid covering all of Ω. Then we evaluate the solution on the slice of Ω defined by the mapping γ (t) of the square in C defined by Re t ∈ [1.66π, 1.76π ] and Im t ∈ [d min , 0.15]. Here we use two grids in t: one uniform 300 × 300 with d min = 10 −3 , and one 300 × 300 that is uniform in Re t and logarithmic in Im t, with d min = 10 −8 . Here we only show the results when using upsampling, since both methods require that to achieve maximum accuracy (without upsampling, HO gets 8 digits and SSQ 9 digits). The results of this test are shown in Fig. 7. Both schemes are able to achieve high accuracy here: both get 13 digits close to the boundary, and both drop to 11-12 digits when tested extremely close to the boundary. However, contrary to the previous test, the proposed SSQ does lose around 1 digit relative to HO, in the distance range 10 −3 to 10 −8 . We do not have an explanation for this loss.
To summarize this experiment, SSQ is much better than HO for panels that are resolved to lower accuracies, because these panels may be quite curved. When panels are resolved to close to machine precision the methods are about the same, with HO attaining a slightly smaller error.

Evaluating the field due to a slender body in three dimensions
To test our method in 3D, we evaluate the Stokes flow around a thin filament in the slender body approximation. In this approximation, the flow due to a fiber with centerline Γ and radius ε 1 is given by the line integral (see e.g. [33]) Here f is given force data, defined on the centerline, S is the Stokeslet kernel, defined as and D is the doublet kernel, defined as Note that S and D are tensors, and I is the 3 × 3 identity. To apply our quadrature, we first write (78) in kernel-split form, where and R := x − y. Once in this form, one can compute modified weights for I 1 , I 3 , and I 5 using the proposed SSQ quadrature. For our tests, we use the force data f ( y) = y and radius ε = 10 −3 . We let Γ be the closed curve shown in Fig. 8, which is described by a Fourier series with decaying random coefficients, where the components of c(k) are complex and normally distributed random numbers. 4 We subdivide Γ into panels by recursively bisecting the curve until each panel is resolved to a tolerance , using a criterion similar to that used in Sect. 5.1: a panel is deemed resolved if the expansion coefficients {ŝ k } of the speed function s(t) = γ (t) satisfy max( ŝ n−1 , ŝ n ) < ŝ ∞ . Each panel is then discretized using n = 16 Gauss-Legendre nodes.

Comparison to adaptive quadrature
As a comparison for our method, and also to compute reference values, we have implemented a scheme for nearly singular 3D line quadrature based on per-target adaptive refinement: For a given target point x, nearby panels are recursively subdivided until each panel Γ i satisfies min y∈Γ i |x − y| < h i , where h i is the arc length of Γ i . Empirically, this is sufficient for accurate evaluation of the slender body kernel on a 16-point panel. Each newly formed panel is discretized using n = 16 Gauss-Legendre points, and the required data (γ , |γ |, f ) is interpolated to these points from the parent panel using barycentric Lagrange interpolation [9]. After this new set of panels is formed, the field (78) can be accurately evaluated at x using the Gauss-Legendre quadrature weights, since all source panels then are far away from x, relative to their own length. This scheme is relatively simple to implement, but has the drawback of requiring an increasingly large amount of interpolations and kernel evaluations as x approaches Γ .
In an attempt to obtain a fair comparison of runtimes between the adaptive quadrature and the proposed SSQ, we will only report on the time spent in the essential core tasks for each method, and ensure that they have been implemented with similar high efficiency. For the adaptive quadrature, we will therefore report the time spent on interpolations (reported as t interp ), since those are implemented using a combination of C and BLAS, and omit the time spent on the recursive subdivisions, since that is presently implemented in a relatively slow prototype code. For the SSQ, we will report on the total time spent on finding roots and computing the new quadrature weights (reported t weights ), since all the expensive steps of that algorithm have been implemented in C. For both methods we will also report on the time spent on near evaluations of the Stokeslet and doublet kernels (reported as t eval ), since those are computed using the same code. Far field evaluation times are omitted, since they are identical. Thus t eval + t interp is somewhat less than the total time for an adaptive scheme, whereas t eval + t weights is the total time for our proposed SSQ.

Results
As a first test, we evaluate the u(x) velocity field (78) (41) to full accuracy]. Unlike in the 2D case, we find that this upsampling is necessary for achieving accurate results. The error is measured against a reference computed using the adaptive quadrature and a discretization with 18-point panels, created using = 5 × 10 −14 (this is to ensure that the grids are different, and that the reference grid is better resolved). The resulting error field, shown in Fig. 8, indicates that the SSQ can recover at least 13 digits of accuracy on all of this grid of target points, and that the largest errors are at the points closest to the curve.
To further study the behavior of SSQ close to the curve, and to compare it to the adaptive quadrature, we use both methods to compute u at a set of random target points all located at the same distance d away from Γ . The test setup and results are listed in Table 1. As expected, the computational costs (kernel evaluations and runtime) of the proposed SSQ quadrature does not vary with the distance between the target points and the curve, as long as the target points are within the near evaluation threshold. The adaptive quadrature, on the other hand, has costs that grow slowly as the distance decreases, since more and more source points are required to evaluate the integral. The SSQ also appears to be a factor of several times cheaper, both in number of evaluations (4-7 times less) and in runtime (2.5-5 times faster), even though our timings are underestimates for adaptive quadrature. However, the error for the SSQ in this application starts to grow for close distances: for instance, at d = 10 −4 the table shows that it loses around 2 digits relative to the at which the panels were discretized.

Remark 11
This loss of accuracy in this application for targets at small distance d is not entirely surprising. The SSQ method gives weights λ at the given nodes for integrating the kernels |R| −1 , |R| −3 , and |R| −5 . The resulting integrals diverge, respectively as log d, d −2 , and d −4 . We have checked that, for these kernels with generic densities, the relative errors of SSQ are close to machine precision (around 13 digits), even though as d → 0 the λ j are growing and oscillatory. However, the Stokes kernels (82)-(84) involve near-vanishing numerator factors (such as R R T ) which partially cancel  the singularities, leading to much smaller values for the integrals. In our method such factors are incorporated into a modified densityf in (44), thus catastrophic cancellation as d → 0 is inevitable for this application (this is also true for the straight fiber method of [45]). This could only be avoided by a more specialized set of recursions for the full Stokes kernels (83) and (84).

Conclusions
We have presented an improved version of the monomial approximation method of Helsing and Ojala [19,23,36] for nearby evaluation of 2D layer potentials discretized with panels. At the cost of a single Newton search for the target preimage, our method uses singularity cancellation to "swap" the singular problem on the general curved (complex) panel for one on the standard (real) panel [−1, 1], much improving the error. Hence we call the method "singularity swap quadrature" (SSQ). We emphasize that it is not the conditioning of the monomial interpolation problem that improvesthe Vandermonde matrix is exponentially ill-conditioned in either case-rather, the exponential convergence rate improves markedly in the case of curved panels. In the case of a Nyström discretization adapted to a requested tolerance, this allows SSQ to achieve this tolerance when HO would demand further subdivision. We gave a quantitative explanation for both the HO and SSQ convergence rates using classical polynomial approximation theory: the HO rate is limited by electrostatic "shielding" of a nearby Schwarz singularity by the panel itself. We then showed that the singularity swap idea gives a close-evaluation quadrature for line integrals on general curves in 3D, which has so far been missing from the literature. The (to these authors) delightful fact that complex analysis can help in a 3D application relies on analytically continuing the squared-distance function (45). Note that QBX [7,29] cannot apply to line integrals in 3D, since the potential becomes singular (in the transverse directions) when approaching the curve. We demonstrated 3D SSQ in a slender body theory Stokes application. We found that, at least when targets are only moderately close, comparable accuracies to a standard adaptive quadrature are achieved with several times fewer kernel evaluations, and, in a simple C implementation, speeds several times faster.
As Remark 11 discusses, in the Stokes application, in order to retain full relative accuracy at arbitrarily small target distances a new set of recurrences for these particular kernels would be essential; we leave this for future work. Of course, the slender body approximation breaks down for distances comparable to or smaller than the body radius, so that the need for full accuracy at such small distances is debatable. and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. f · ( y − x)( y − x)( y − x) ·n y − x 4 ds( y) As an example of how these can be used, if we let f =n in (91) and use that νν = 1, we directly get the complex form of the double layer kernel of the Laplace equation,