Quadrature algorithms for high dimensional singular integrands on simplices

Abstract

Galerkin discretizations of integral operators in $\mathbb {R}^{d}$ require an accurate numerical evaluation of integrals $I={\int }_{\!\!S^{(1)}}{\int }_{\!\!S^{(2)}}f(x,y)dydx$ where S ⁽¹⁾, S ⁽²⁾ are d-simplices and the integrand function f has a possibly nonintegrable singularity at x = y. In a previous paper (2011) we constructed several families of quadrature rules ${Q}_{\mathcal {N}}$ for a class of functions f which allow algebraic singularities at x = y, including hypersingular kernels, and are Gevrey smooth for x ≠ y. This holds for kernel functions commonly occurring in integral equations. In this paper we address the implementation aspects for computing $Q_{\mathcal {N}}$ and give a detailed computation algorithm for arbitrary space dimension d and arbitrary mutual location of simplices S ⁽¹⁾ and S ⁽²⁾. The algorithm consists of a “desingularizing” coordinate transformation, which reduces the singular support of the integrand to one variable while preserving Gevrey regularity in all 2d−1 remaining variables and simultaneously simplifies the integration domain to a unit cube. Due to the simple singularity structure after transformation, one can optionally use various combinations of quadrature rules in singular and regular directions. We report on comprehensive convergence studies for the full tensor product and Smolyak-type quadrature rules in the 2d−1 Gevrey-regular variables combined with either composite Gauss-Legendre or Gauss-Jacobi quadrature rules in the singular direction. A Matlab software implementing the algorithm complements this paper and is available from Netlib.