Error estimate based adaptive quadrature for layer potentials over axisymmetric surfaces

TL;DR

The paper addresses the challenge of evaluating nearly singular layer potentials on axisymmetric surfaces. It introduces the Singularit y Swap Surface Quadrature (S3Q), a two-stage adaptive scheme that first performs azimuthal singularity swapping with a Fourier-based representation and then handles the polar integral via adaptive, semi-analytical quadrature, supported by rigorous complex-analytic error estimates. Key contributions include automatic parameter selection for adaptivity, new quadrature and interpolation error bounds, and an extension of the SSQ framework to 3D axisymmetric geometries, with demonstrated accuracy for harmonic and Stokes potentials on spheroidal geometries. The approach yields controlled accuracy while maintaining manageable computational cost, and it provides a solid foundation for extensions to Helmholtz problems and near-surface evaluations in more complex axisymmetric configurations. Near very close evaluations, the method faces cancellation-driven limitations, suggesting complementary strategies (e.g., asymptotics or hedgehog-type extrapolation) for the thinnest near-surface layer.

Abstract

Layer potentials represent solutions to partial differential equations in an integral equation formulation. When numerically evaluating layer potentials at evaluation points close to the domain boundary, specialized quadrature techniques are required for accuracy because of rapid variations in the integrand. To efficiently achieve a specified error tolerance, we introduce an adaptive quadrature method with automatic parameter adjustment for axisymmetric surfaces, facilitated by error estimation. Notably, while each surface must be axisymmetric, the integrand itself need not be, allowing for applications with complex geometries featuring multiple axisymmetric bodies. The proposed quadrature method utilizes so-called interpolatory semi-analytical quadrature in conjunction with a singularity swap technique in the azimuthal angle. In the polar angle, such a technique is used as needed, depending on the integral kernel, combined with an adaptive subdivision of the integration interval. The method is tied to a regular quadrature method that employs a trapezoidal rule in the azimuthal angle and a Gauss-Legendre quadrature rule in the polar angle, which will be used whenever deemed sufficiently accurate, as determined by a quadrature error estimate [C. Sorgentone and A.-K. Tornberg, Advances in Computational Mathematics, 49 (2023), p. 87]. Error estimates for both numerical integration and interpolation are derived using complex analysis, and are used to determine the adaptive panel subdivision given the evaluation point and desired accuracy. Numerical examples are presented to demonstrate the method's efficacy.

Figures (11)

• Figure 1: Colors display the measured error \ref{['eq:reg_quad_err']} in $\log_{10}$ scale, for the Laplace layer potentials, evaluated with ${n_t}={n_\varphi}=n$, and only one interval in $t$. The solid black contour lines illustrating the error estimate from SORGENTONE2023 at levels $10^{-2},~10^{-4},~10^{-6},~10^{-8}$ correspond well to the measured error. The layer density function is set to $\sigma(\boldsymbol{\gamma}(\theta,\varphi))=\sin(5\theta)e^{-\cos^2(\varphi)}+1.03$.
• Figure 1: Typical behavior of the functions appearing in the integrand of \ref{['eq:layer_potential_abc']}. In panels (a) and (c) the functions are shown and $\theta=\operatorname{Re}(\theta_0^\lambda)$, for this specific target point, is marked. Panels (b) and (d) show the relative decay of their Chebyshev coefficients. The quantity $|\theta-\theta_0|^{2p-1}(R_\lambda^2)^{-(p-1/2)}$ is plotted in green in panel (c), with the decay of its Chebyshev coefficients in panel (d) showing its low frequency content. Here, $\theta_0^\lambda$ is a root of $R_\lambda^2$ defined by \ref{['eq:theta0_spheroid']} and \ref{['eq:Rlambda']}, respectively. In this example, $k({\mathbf{x}},{\mathbf{y}})=1$ and $\sigma(\boldsymbol{\gamma}(\theta,\varphi))=\sin(5\theta)e^{-\cos^2(\varphi)}+1.03$.
• Figure 1: Sketch of the contour $C$, the base interval $E$, interpolation point $\tau$, and the deformations $C_1=C_1^-\cup C_1^+$ and $C_2=C_2^-\cup C_2^+$ circumventing the singularities $t_0$ and $\overline{t_0}$ with branch cuts $B(t_0)$ and $B(\overline{t_0})$, respectively.
• Figure 1: Harmonic single layer potential computed using the singularity swap surface quadrature (S3Q) with different error tolerances $\epsilon$. Panel (a) shows how the measured error for each target point where S3Q was used stay close to $\epsilon$. Panel (b) illustrates the measured error and the contour level where the regular quadrature error estimate equals $\epsilon=10^{-6}$. Panel (c) shows how the distribution of ${n_{\textrm{GL}}}$-point Gauss--Legendre panels, $n_{\textrm{pan}}$, set by the adaptive algorithm in the polar direction varies with $\epsilon$.
• Figure 1: Panel (a) shows the value of $\mu_k^{3/2}(\alpha)$ and panel (b) the associated absolute forward recurrence error. The black lines, illustrating the corresponding estimates in Table \ref{['tab:decay_rate_rec_err']} at levels $10^{-2},~10^{-4},~10^{-6},~10^{-8},~10^{-10},~10^{-12}$, agree well with their two counterparts.

Theorems & Definitions (26)

• Lemma 2.1: Root of $R^2$ in $\varphi$ SORGENTONE2023
• Lemma 2.2
• Lemma 2.3
• Remark 2.4
• Remark 2.5: Accuracy and stability of recurrence formulas for $\mu_k^p$
• Theorem 2.6
• Remark 2.7: The function $\mathcal{F}$
• Lemma 2.8
• Remark 2.9: Stability and accuracy of recurrence formulas $\nu_k^p$
• Remark 2.10: Arbitrary function interpolation and stability