High Order Regularization of Nearly Singular Surface Integrals

Basic Information

• Paper ID: 2510.13639
• Title: High Order Regularization of Nearly Singular Surface Integrals
• Authors: J. Thomas Beale (Duke University), Svetlana Tlupova (Farmingdale State College, SUNY)
• Classification: math.NA, cs.NA (Numerical Analysis)
• Publication Date: October 16, 2025
• Paper Link: https://arxiv.org/abs/2510.13639

Abstract

Solutions to partial differential equations can often be represented as surface integrals with kernel functions related to singular fundamental solutions. Accurate evaluation of integrals at points on or near the surface requires special methods. Building on previous analytical work, this paper derives formulas for high-precision regularized integrals that enable the use of standard quadrature methods without special treatment near singular points. The study covers single-layer and double-layer integrals for harmonic functions and Stokes flow. Nearly singular cases (evaluation at points close to the surface) are necessary when surfaces are in close proximity or when evaluating at grid points near the surface. The paper derives regularized kernel formulas with errors O(δᵖ), where δ is the smoothing radius and p = 3, 5, 7. By choosing δ = κhᵍ (q < 1), discretization errors are controlled as h → 0, with expected convergence rates O(hᵖᵍ) observed in various examples.

Research Background and Motivation

Problem Description

1. Core Problem: In boundary integral methods for partial differential equations, when evaluation points lie on or near the surface, the integral kernel functions become singular or nearly singular, causing numerical computation difficulties.
2. Problem Significance:
   • Boundary integral methods are widely applied in solving elliptic PDEs, Stokes flow, and other problems
   • Nearly singular integrals are unavoidable when handling surfaces in close proximity or computing values at grid points near surfaces
   • Accurate evaluation of these integrals is critical for overall solution accuracy
3. Limitations of Existing Methods:
   • Singularity subtraction techniques require analytical computation of the most singular part
   • Extrapolation methods (QBX) and hedgehog methods are computationally complex
   • Simple regularization methods have limited accuracy
   • Existing high-order methods often require multiple computations and extrapolation
4. Research Motivation: Develop a high-order regularization method that can:
   • Achieve high precision with a single parameter δ
   • Use standard quadrature rules
   • Apply to evaluation points on and near surfaces
   • Provide controllable convergence

Core Contributions

1. High-Order Regularization Formulas: Derives 3rd, 5th, and 7th order regularized kernel functions with errors O(δ³), O(δ⁵), and O(δ⁷) respectively
2. Unified Framework: Provides unified regularization methods for single-layer and double-layer potentials of harmonic functions, as well as Stokeslet and stresslet integrals in Stokes flow
3. Parameter Selection Strategy: Proposes parameter choice δ = κhᵍ where q < 1, achieving convergence rate O(hᵖᵍ) in total error
4. Efficient Extension Method: Combines Mayo's method to efficiently compute solutions on the entire grid from integral values near the surface
5. Practical Application Verification: Validates the method's effectiveness in harmonic function interface problems and Stokes flow

Method Details

Task Definition

Input:

• Density functions f(x) or g(x) on surface Γ
• Evaluation point y (possibly on or near the surface)
• Grid spacing h

Output:

• High-precision numerical values of regularized surface integrals
• Solution function values over the entire computational domain

Constraints: Evaluation point y can be expressed as y = x₀ + bn, where x₀ is the nearest point on Γ, n is the outward normal vector, and b is the signed distance

Core Regularization Strategy

1. Basic Regularization Idea

Replace the singular kernel G(r) = -1/(4π|r|) with a smoothed version:

Gδ(r) = G(r)s₁(|r|/δ)

where s₁(ρ) = erf(ρ) is the error function.

2. High-Order Correction Method

Based on error expansion:

Sδ(y) = S(y) + C₁δI₀(b/δ) + C₂δ³I₂(b/δ) + C₃δ⁵I₄(b/δ) + O(δ⁷)

Eliminate leading error terms by modifying the smoothing factor s₁ to s₁⁽ᵖ⁾:

s₁⁽⁷⁾(ρ) = erf(ρ) + (c₁ρ + c₂ρ³ + c₃ρ⁵)e^(-ρ²)

3. Regularization of Various Integral Types

Single-Layer Potential Integral:

S(y) = ∫_Γ G(x-y)f(x)dS(x)

Regularized using modified s₁⁽ᵖ⁾(ρ).

Double-Layer Potential Integral:

D(y) = ∫_Γ ∂G(x-y)/∂n(x) g(x)dS(x)

Regularized using subtraction form and modified s₂⁽ᵖ⁾(ρ).

Stokeslet Integral:

uᵢ(y) = 1/(8π) ∫_Γ Sᵢⱼ(y,x)fⱼ(x)dS(x)

Regularized using combination of s₁⁽ᵖ⁾ and s₂⁽ᵖ⁾.

Stresslet Integral:

vᵢ(y) = 1/(8π) ∫_Γ Tᵢⱼₖ(y,x)qⱼ(x)nₖ(x)dS(x)

Requires additional s₃⁽ᵖ⁾(ρ) function.

Technical Innovations

1. Systematic Correction: Determines correction coefficients by solving linear systems, avoiding complex analytical calculations
2. Unified Mathematical Framework: Discovers that correction coefficients for different integral types satisfy the same linear system, simplifying implementation
3. Extensible Design: Lower-order versions can be obtained by simple truncation of higher-order formulas
4. Surface Specialization: Provides simplified formulas for evaluation points on the surface, improving computational efficiency

Experimental Setup

Test Problems

1. Harmonic Functions on Spheres: Uses spherical harmonics to construct known solutions
2. Molecular Surfaces: Four-atom molecular model with complex geometry
3. Ellipsoids: Ellipsoids with varying aspect ratios
4. Stokes Flow of Translating Sphere: Standard test problem with analytical solution
5. Two Closely Spaced Spheres: Tests nearly singular cases

Evaluation Metrics

• Maximum Error: max|u_computed - u_exact|
• L₂ Error: (∑|error|²/N)^(1/2)
• Convergence Rate: Determined by comparing errors at different grid spacings h

Parameter Settings

• Regularization Parameter: δ = κh^q, where:
  • p=3: q=2/3, κ₀=1,2
  • p=5: q=4/5, κ₀=2,3
  • p=7: q=5/7, κ₀=3,4
• Grid Spacing: h from 1/32 to 1/1024
• Quadrature Rule: Wilson's spherical partition quadrature

Implementation Details

• Uses GMRES to solve integral equations (tolerance 10⁻¹⁰)
• Ignores regularization effects for points beyond distance 8δ
• Employs treecode acceleration for far-field computations
• Fourth-order discrete Laplacian for grid extension

Experimental Results

Main Results

1. Harmonic Function on Sphere Verification

• δ = 4h: Observes O(h⁷) convergence
• δ = 2h^(5/7): Achieves O(h⁵) convergence
• Results consistent for single-layer and double-layer potentials

2. Convergence on Complex Geometries

Molecular surface test results:

• 3rd order kernel: O(h²) convergence (theoretical expectation O(h^(pq)) = O(h²))
• 5th order kernel: O(h⁴) convergence
• 7th order kernel: O(h⁵) convergence

Ellipsoid tests:

• Verifies stability across different geometric shapes
• Achieves optimal O(h⁵) precision with κ₀ = 4

3. Stokes Flow Application

Stokes flow of translating sphere:

• Pressure and velocity integrals: O(h⁵) precision
• Solutions on grid: O(h⁴) precision
• Velocity gradients: O(h⁴) precision

4. Nearly Singular Cases

Two spheres with separation ε = 1/163:

• 5th order method: Stable O(h⁴) convergence
• 7th order method: Achieves O(h⁵) convergence
• Error at near-contact points comparable to edge point errors

Ablation Studies

Impact of Parameter Choice

• Excessively small κ₀ values lead to unstable convergence
• κ₀ = 2, 3, 4 respectively achieve optimal performance for p = 3, 5, 7
• Choice of q value directly affects overall convergence rate

Comparison of Regularization Orders

• Higher-order methods allow larger δ values
• High-order methods significantly improve precision at same grid resolution
• Computational cost increase is reasonable relative to precision improvement

Case Studies

Grid Extension Method Verification

Computing entire grid values from integral values near surface using Mayo's method:

1. Integral values near surface: O(h⁵) precision
2. Grid function values: O(h⁴) precision
3. First-order differences: O(h⁴) precision

This "precision degradation" is due to:

• Truncation error of discrete Laplacian operator
• Propagation effects of errors near the surface

Related Work

Main Research Directions

1. Singularity Subtraction Methods: Helsing et al.'s high-order subtraction techniques
2. QBX Methods: Klöckner et al.'s quadrature by expansion methods
3. Local Correction Methods: Nitsche et al.'s trapezoidal rule corrections
4. Heat Potential Methods: Asymptotic analysis-based methods
5. Harmonic Density Interpolation: Pérez-Arancibia et al.'s methods

Advantages of This Work

1. Single-Parameter Method: Simpler than extrapolation methods
2. Unified Framework: Applicable to multiple integral types
3. High-Order Accuracy: Achieves up to 7th order regularization
4. Practicality: Easy to implement and apply

Conclusions and Discussion

Main Conclusions

1. Successfully derives 3rd, 5th, and 7th order regularization formulas achieving expected high precision
2. Parameter selection strategy δ = κh^q effectively controls total error
3. Combined with Mayo extension method, efficiently solves entire computational domain
4. Method demonstrates stable performance on complex geometries and nearly singular cases

Limitations

1. Geometric Requirements: Requires sufficiently smooth surfaces
2. Parameter Tuning: κ value selection requires empirical guidance
3. Computational Complexity: High-order formulas increase implementation complexity
4. Theoretical Analysis: Lacks rigorous convergence proofs

Future Directions

1. Develop adaptive parameter selection strategies
2. Extend to more complex PDE problems
3. Combine with fast algorithms to improve efficiency
4. Theoretical convergence analysis

In-Depth Evaluation

Strengths

1. Method Innovation: Systematic high-order regularization avoiding multiple extrapolations
2. Practical Value: Provides complete implementation guidance and parameter selection strategies
3. Comprehensive Verification: Full testing from simple geometries to complex applications
4. Clear Writing: Rigorous mathematical derivations with detailed implementation details

Weaknesses

1. Theoretical Gaps: Lacks rigorous error analysis and convergence proofs
2. Parameter Dependence: κ value selection still requires experience; lacks automation strategies
3. Limited Scope: Primarily restricted to smooth surfaces; may not apply to geometries with sharp corners
4. Computational Cost: Insufficient analysis of computational overhead for high-order formulas

Impact

1. Academic Contribution: Provides new high-precision tools for boundary integral methods
2. Application Prospects: Broad application potential in computational fluid dynamics, electromagnetics, and other fields
3. Reproducibility: Provides detailed implementation details and open-source code

Applicable Scenarios

1. High-precision boundary integral computations
2. Handling multiple surfaces in close proximity
3. Numerical solution of Stokes flow and potential flow problems
4. High-precision solution of interface problems

References

The paper cites 31 related references, primarily including:

• 5 Beale & Tlupova (2024): Theoretical foundation of this work
• 18 Mayo (1985): Grid extension method
• 30 Wang et al. (2020): Treecode acceleration algorithm
• 6 Beale et al. (2016): Wilson quadrature method

This paper makes important contributions to the field of numerical analysis, providing practical and effective solutions for high-precision computation of nearly singular surface integrals. The systematic and comprehensive nature of the method demonstrates excellent application prospects.