PaperlayerView exact revision on arXiv

Source-linked AI summary

Quadrature by Matched Asymptotic Expansion (QBMAX): Accelerating Evaluation of Layer Potentials with Boundary Layers

Chaoqi Lin, Xiaoyu Wei, Andreas Klöckner

arXiv:2609.10857v2math.NAcs.NA

TL;DR

Near-boundary layer-potential evaluation is difficult when rapidly decaying kernels vary sharply near their singularities. QBMAX extends QBX by factoring the kernel’s leading asymptotic decay into the Taylor expansion; in 2D and 3D tests, it improves accuracy by up to four digits in some large-parameter cases, with similar performance in some small-parameter or curved-geometry cases.

  • Problem

    Near-boundary layer-potential evaluation is difficult because the kernel is singular or nearly singular and can vary rapidly near the singularity, while a fixed-order QBX Taylor polynomial may poorly capture rapid decay at large k.

    Source: §1 · p. 2

  • Method

    QBMAX extends QBX by removing the leading exponential decay along the expansion line before Taylor approximation.

    Source: §1 · p. 3

  • Results

    QBMAX improves accuracy over QBX by up to four digits in some large-parameter cases, while performing similarly in some small-parameter or curved-geometry cases.

    Source: §1 · p. 4

  • Takeaways & Limitations

    The flat-boundary analysis shows reduced truncation-error constants, and numerical tests generally show less convergence degradation as k increases, with improvement depending on density, geometry, and discretization.

    Source: §5 · p. 30

    The authors do not yet have a theoretical error analysis for curved geometries, and QBMAX gives only marginal improvement over QBX at small k on the starfish domain, possibly due to local curvature.

    Source: §4.3.2 · p. 22

Abstract

from arXiv · show

Evaluating layer potentials with rapidly decaying kernels becomes computationally challenging in the presence of boundary layers, as occurs for Helmholtz layer potentials with large imaginary wavenumbers. We introduce Quadrature by Matched Asymptotic Expansion (QBMAX), a high-order method that extends the Quadrature by Expansion (QBX) framework by incorporating the kernel's asymptotic decay behavior directly into the Taylor expansion of the potential. Numerical results in both 2D and 3D show improvements of up to four digits over QBX for the tested large-parameter cases, while the methods can perform similarly for smaller parameters or more strongly curved geometries. A symbolic weighted operation-count model indicates comparable costs for forming the two expansion kernels. The Taylor expansions have the same formal order, while QBMAX reduces the truncation-error constants in the flat-boundary model and generally exhibits less convergence degradation in the reported experiments.

1. Introduction

QBMAX extends QBX by incorporating the kernel’s leading asymptotic decay into the Taylor expansion. In the reported 2D and 3D tests, it improves accuracy by up to four digits for some large-parameter cases, while performing similarly in some small-parameter or curved-geometry cases.

1. Introduction

  • QBMAX incorporates the kernel’s leading asymptotic decay into the Taylor expansion to reduce truncation error for rapidly decaying potentials.It removes the leading exponential decay along the expansion line before Taylor approximation; the resulting expansion is target-specific.Source: §1 · p. 3
  • The authors describe applying modified Taylor expansions to rapidly decaying kernels and integral equations as novel to the best of their knowledge.They distinguish this application from related uses in singularity removal and function approximation [14].Source: §1 · p. 3
  • A symbolic weighted operation-count model indicates comparable costs for forming the QBMAX and QBX expansion kernels.Across reported refinement studies, the same meshes are used for all tested k values, and QBMAX is substantially less sensitive to k in the rapidly decaying regime.Source: §1 · p. 4
  • For an infinite flat boundary in 2D, the analysis makes explicit how QBMAX reduces truncation-error constants while retaining the same asymptotic convergence rate as QBX.The analysis is presented to clarify error behavior in more general settings.Source: §1 · p. 4
  • QBMAX improves accuracy by up to four digits over QBX in some large-parameter cases, but performs similarly in some small-parameter or curved-geometry cases.Its observed convergence rates generally degrade less as k increases, although the rates vary with geometry.Source: §1 · p. 4

2. Fundamental Solutions and Layer Potentials

The modified Helmholtz fundamental solutions use modified Bessel functions whose rapid decay concentrates the kernel near its singularity at large parameter values. The section introduces single- and double-layer potentials and explains how QBX uses off-surface Taylor expansions for close evaluation, with accuracy and efficiency shaped by truncation and quadrature choices.

2.1. Fundamental Solutions

  • For large k, the modified Bessel kernel decays rapidly away from x = y and is concentrated in a thin region around the singularity.Near z = 0, K0(z) and its derivative have leading singular behavior −ln(z) and −1/z, respectively.Source: §2.1 · pp. 4–5

2.2. Layer Potentials

  • The section defines single- and double-layer potentials on smooth boundaries, with the Green’s function’s leading singularity matching the Laplace kernel.The single-layer potential is continuous, while the double-layer potential has jump relations.Source: §2.2 · p. 5

2.3. Quadrature by Expansion and Taylor Expansions of Layer Potential on a Line

  • QBX evaluates near-boundary targets by approximating the potential along a line from an off-surface expansion center with a Taylor expansion, then integrating its expansion kernel.This regularizes the kernel singularity and allows standard quadrature; the approach relies on analyticity of layer potentials up to the boundary under appropriate regularity assumptions.Source: §2.3 · p. 6
  • QBX has truncation and quadrature errors; practical implementations refine geometry and quadrature resolution, balancing expansion radius and order against accuracy and efficiency.The reported quadrature-error estimates include smooth closed curves in 2D and extensions to 3D; the original QBX uses volumetric rather than line expansions.Source: §2.3 · pp. 6–7
  • QBX relative error and pre-asymptotic convergence can deteriorate as k increases, reducing efficiency for prescribed relative accuracy in the rapidly decaying regime.Expansion radius and order affect truncation error, while reducing the radius or increasing the order requires higher quadrature resolution.Source: §2.3 · p. 7

3. Quadrature by Matched Asymptotic Expansion

The flat-boundary analysis motivates QBMAX by showing how rapidly decaying kernels can produce large relative truncation errors in fixed-order QBX. QBMAX removes the leading decay before Taylor approximation, preserving formal order while reducing truncation-error constants in the analyzed setting.

3.1. Infinite Flat Boundary: Truncation Error Analysis

  • For fixed p, the relative QBX error approaches one as the decay parameter grows, while maintaining accuracy by shrinking the expansion radius requires finer quadrature.The flat-boundary model indicates that this can increase computational cost, particularly in 3D.Source: §3.1 · p. 8
  • For constant density (ω = 0), the QBMAX expansion has zero truncation error and is exact in the infinite-flat-boundary analysis.The analysis considers a Fourier-mode density on an infinite flat boundary in 2D.Source: §3 · p. 7, §3.1 · p. 9
  • For fixed p and ω′, QBMAX’s relative truncation error tends to zero as k′ grows, while QBX’s remains O(1) in the flat-boundary model.Figure 4 plots the error ratio for p = 4, 6, 8; its log10(R) = 1 contour marks a tenfold reduction.Source: §3.1 · p. 9, Figure 4 · p. 10

3.2. QBMAX Expansion

  • QBMAX forms a target-specific Taylor expansion after scaling the potential to compensate for its exponential decay along the expansion line.The scaling function is included in the expansion, then divided out to recover the original layer potential.Source: §3.2 · p. 10

3.3. QBMAX Quadrature Error

  • QBMAX retains the QBX expansion’s formal asymptotic order, (rτ)^(p+1), but removes the leading exp(−kr(1 −τ)) decay factor before Taylor approximation.The quadrature-error argument retains the same order in h for fixed parameters, without asserting a constant uniform in k or kr.Source: §3.2 · p. 11, §3.3 · p. 12

4. Numerical Results

The numerical-results section first compares the cost of forming QBX and QBMAX expansion kernels. Under the stated weighted symbolic operation-count model, their costs are comparable; the fitted exponents describe only the tested range.

4.1. Cost Analysis

  • The cost estimate counts exponentials as 10 flops, other basic function evaluations as 1 flop, and uses SymPy derivatives with common-subexpression elimination.The exponent is fitted by linear least squares to log(weighted operation count) against log(p) for p = 5, …, 8.Source: §4.1 · p. 12
  • Under the weighted symbolic operation-count model, forming QBMAX expansion kernels costs about as much as forming QBX kernels.Figure 5 reports counts for 2D and 3D kernels, with an empirical power-law fit over 5 ≤ p ≤ 8.Source: §4.1 · p. 13, Figure 5 · p. 13
  • The fitted exponent is lower in 3D than in 2D, indicating more effective SymPy common-subexpression elimination for the tested 3D expressions.The authors caution that this short-range fit is not an asymptotic complexity result.Source: §4.1 · p. 13

4.2. Problem Setup for Layer Potential Evaluation

The setup evaluates regularized layer potentials on discretized boundaries using composite quadrature and upsampling. Tests cover close evaluation and jump relations across 2D and 3D geometries, with convergence assessed under mesh refinement.

4.2. Problem Setup for Layer Potential Evaluation

  • Boundaries are divided into panels and integrated with q-point Gauss–Legendre rules in 2D or q²-point tensor-product rules in 3D.The authors note that the method also applies to quadrature rules on triangles.Source: §4.2 · p. 13
  • Quadrature is upsampled by κ to resolve the kernel, distinguishing discretization nodes from the denser source nodes used in the regularized-potential evaluation.The Nyström discretization uses panelwise interpolation and source-node quadrature weights.Source: §4.2 · p. 13
  • The experiments use the unit disk and starfish domain in 2D and the cruller surface in 3D, with analytic or self-convergence references as specified.The listed test cases include close evaluation and numerical verification of jump relations.Source: §4.2 · pp. 14–15
  • Close evaluation samples targets from expansion-center to on-surface positions along interior and exterior expansion lines; successive mesh refinements support h-convergence assessment.At τ = 1, the values are interpreted as interior and exterior limits; h is the maximum panel arc length in 2D and maximum square root of panel area in 3D.Source: §4.2 · pp. 14–15

4.3. 2D Numerical Results – Forward Evaluation

The disk tests show stronger accuracy and convergence stability for QBMAX, especially as k increases. On the starfish, gains are smaller for low k and curved geometry, but become clear for rapidly evanescent Helmholtz potentials.

4.3.1. Disk

  • 13-digit accuracy at expansion centers and about 8-digit accuracy on the boundary were achieved by QBMAX for k=80 in disk close evaluation.QBMAX consistently outperformed QBX across the tested k values.Source: Figure 7 · p. 16, §4.3.1 · p. 16
  • For on-surface single-layer evaluation, QBMAX convergence rates stayed near 6 (interior) and 5.6 (exterior), while QBX rates fell to about 3.5 at k=80.The QBX rates at k=10 were 5.61 (interior) and 5.46 (exterior).Source: §4.3.1 · p. 16, Figure 8 · p. 17
  • For on-surface double-layer evaluation, QBMAX maintained a global convergence rate near 4.8, while QBX declined from 4.8 at k=10 to 3.6 at k=80.QBMAX's rate varied minimally across the tested k values.Source: Figure 10 · p. 18, §4.3.1 · p. 18
  • Before errors plateaued around 12 digits, QBMAX jump-relation convergence rates remained above 6; the plateau was likely due to floating-point limits in computing small distances.The reported jump relations were for Dkσ and S′kσ.Source: §4.3.1 · p. 18, Figure 11 · p. 19
  • For oscillatory disk densities, the QBX-to-QBMAX error gap decreased as frequency m increased; for m≥40, the flat-panel gamma ratio approximated the error ratio within roughly one order of magnitude.In the low-frequency regime m≪k, the gamma ratio diverges because QBMAX has zero truncation error on flat boundaries with constant density.Source: §4.3.1 · p. 19, Figure 12 · p. 20

4.3.2. Starfish

  • For starfish close evaluation at k=10, the methods showed no significant difference for the double layer and at most one order of magnitude difference for the single layer.Both single-layer methods had slopes of approximately 7; at small k, QBMAX's improvement may be marginal on this curved geometry.Source: §4.3.2 · pp. 21–22, Figure 14 · p. 22
  • For starfish jump relations, both methods converged at approximately order 6; the reported analysis has no theoretical error analysis for curved geometries, where curvature may limit improvement.The authors suggest local curvature correction terms as a direction for future work.Source: Figure 15 · p. 22, §4.3.2 · p. 22
  • When Im(k)>Re(k), QBMAX consistently outperformed QBX for starfish Helmholtz close evaluation, often gaining at least one digit, with the gap growing as Im(k) increased.For Im(k)≤Re(k), the methods showed no significant performance difference in the reported cases.Source: Figure 16 · p. 25, §4.3.2 · p. 25

4.4. 2D Numerical Results – Interior Neumann Problem

The starfish interior Neumann problem is formulated as a second-kind integral equation and solved with Nyström discretization and GMRES. The reported refinement tests show the expected high-order convergence and low iteration counts outside underresolved cases.

4.4. 2D Numerical Results – Interior Neumann Problem

  • The Neumann problem was represented by a second-kind integral equation and discretized with the Nyström method; GMRES solved the resulting system for the boundary density.The solution was represented as a single-layer potential with an unknown density.Source: §4.4 · p. 26
  • Convergence orders ranged from about 6.9 to 7.6 for the interior Neumann problem, and GMRES required about 7–8 iterations outside underresolved cases.The iteration count used a relative residual tolerance of 10^-12.Source: §4.4 · p. 26, Figure 18 · p. 27
  • Relative error was evaluated at τ=1 and τ=0.5 across panel-size refinements.The figure reports relative error versus element size and GMRES iteration count versus boundary-element count.Source: §4.4 · p. 26, Figure 18 · p. 27

4.5. 3D Numerical Results

The 3D cruller tests examine close evaluation and jump relations for QBMAX and QBX. QBMAX showed markedly stronger jump-relation convergence at large k, while errors on the surface varied with local curvature.

4.5.1. Cruller

  • For close evaluation of the single-layer potential on the cruller, both methods had slopes of approximately 6.The test used 5000 panels, with expansion radii from 0.047 to 0.11 depending on local geometry.Source: Figure 19 · p. 28, §4.5.1 · p. 28
  • For cruller jump relations, QBMAX achieved global convergence rates from 3.5 to 5.9, while QBX stagnated near 1.6 at k=80.The reported convergence rates are summarized in Table 3.Source: §4.5.1 · p. 28, Table 3 · p. 29
  • Mean jump-relation error was around 10^-6 for QBMAX and 10^-2 for QBX, with QBMAX's highest error in the high-curvature inner region.QBX's maximum error occurred in the outer region; the authors identify local curvature correction as future work.Source: §4.5.1 · p. 28, Figure 21 · p. 30

5. Conclusions

The conclusion presents QBMAX as a high-order method for close evaluation with rapidly decaying kernels, particularly Yukawa potentials in 2D and 3D. It reports reduced flat-boundary truncation error and generally less k-dependent convergence degradation in the numerical tests, while identifying theoretical and computational extensions as future work.

5. Conclusions

  • QBMAX targets close evaluation of layer potentials with rapidly decaying kernels, particularly Yukawa potentials in 2D and 3D.The method factors out the leading exponential decay; its improvement depends on density, geometry, and discretization parameters.Source: §5 · p. 30
  • In the infinite-flat-boundary model, QBMAX reduces relative truncation error and generally reduces the convergence degradation observed for QBX as k increases in the numerical tests.The conclusion does not claim these effects are uniform across densities, geometries, or discretizations.Source: §5 · p. 30
  • Future work includes extending rigorous asymptotic analysis to curved surfaces, integrating fast algorithms such as the Fast Multipole Method, and applying the principles to other parameter-dependent PDEs.The conclusion names modified biharmonic equations and modified Stokes flows as examples.Source: §5 · p. 30

Appendix A. Integral identities

Appendix A derives closed-form integrals involving modified Bessel functions of the second kind for Fourier-mode layer densities on the unit circle. It uses Graf’s addition theorem and continuity to evaluate the single-layer potential on the circle.

Appendix A. Integral identities

  • Graf’s addition theorem represents the modified Bessel kernel for points inside or outside the unit circle, supporting the closed-form integral identities.The appendix treats both |x| < 1 and |x| > 1 cases for the Fourier-mode density.Source: Appendix A · pp. 31–32
  • For the Fourier-mode density σ_m(θ) = e^{imθ}, the single-layer potential on |x| = 1 is e^{imθ_x}I_m(k)K_m(k).The boundary value follows from continuity of S_kσ.Source: Appendix A · pp. 31–32
  • For the radial derivative calculation, the directional derivative of e^{ilθ} along a radial ray vanishes.This makes the corresponding inner-product term zero.Source: Appendix A · p. 32

Appendix B. Derivative bounds for K0

Appendix B establishes derivative bounds for K0 using its series and integral representations. Near zero, derivatives have singularity no stronger than the stated inverse-power or logarithmic bound.

Appendix B. Derivative bounds for K0

  • The appendix invokes the series expansion and integral representation of K0 to derive bounds on its derivatives.It also uses the expansion of cosh(t) − 1; the resulting constants depend only on derivative order.Source: Appendix B · pp. 32–33
  • For 0 < z < 1, the m-th derivative of K0(z) is bounded by C_m max{z^-m, |log z|}.The singularity is of order O(max{z^-m, |log z|}), with C_m depending on m.Source: Appendix B · p. 32

Appendix C. Quadrature error analysis for the QBX expansion

Appendix C bounds quadrature error for the QBX line expansion under smoothness, panelization, and clearance assumptions. It combines Gauss–Legendre quadrature estimates with derivative bounds for K0 to bound individual terms and the summed pth-order expansion error.

Appendix C. Quadrature error analysis for the QBX expansion

  • The analysis assumes a C∞ density and boundary, equal-size non-overlapping panels, and panel parametrizations inherited from a fixed smooth global parametrization.These assumptions keep panel-parameterization derivatives uniformly bounded under refinement.Source: Appendix C · p. 33
  • The QBX line expansion is evaluated at a boundary target with an expansion center satisfying the standard clearance condition.The condition requires the radius-r ball about the center to intersect the boundary only at the target, so |c_x − y| ≥ r for every boundary point y.Source: Appendix C · p. 33
  • Cartesian derivatives of K0(k|z − y|) contribute at most one additional inverse power of |z − y| each.For larger separation, the exponential-decay estimate is bounded and absorbed into the same constant.Source: Appendix C · p. 34
  • Each m-th QBX expansion term receives a quadrature error bound using a q-point Gauss–Legendre rule and a factor proportional to h^{2q+1}.The bound follows from a quadrature error estimate and the Leibniz rule; its constant depends on k, m, q, and ∂Ω.Source: Appendix C · p. 34
  • The summed quadrature error for the pth-order QBX expansion is bounded by a constant depending on k, p, q, and ∂Ω.The derivation sums panelwise error bounds.Source: Appendix C · p. 35
Loading 2609.10857v2…