Locally Corrected Multidimensional Quadrature Rules for Singular Functions

Accurate numerical integration of singular functions usually requires either adaptivity or product integration. Both interfere with fast summation techniques and thus hamper large-scale computations. This paper presents a method for computing highly accurate quadrature formulas for singular functions which combine well with fast summation methods. Given the singularity and the N nodes, we first construct weights which integrate smooth functions with order-k accuracy. Then we locally correct a small number of weights near the singularity, to achieve order-k accuracy on singular functions as well. The method is highly efficient and runs in O(Nk²ᵈ + Nlog ² N) time and O(k²ᵈ + N) space. We derive precise error bounds and time estimates and confirm them with numerical results which demonstrate the accuracy and efficiency of the method in large-scale computations. As part of our implementation, we also construct a new adaptive multidimensional product Gauss quadrature routine with an effective error estimate, and compare it with a standard package. The approach generalizes to interpolate and differentiate scattered data and to integrate singular functions over curves and surfaces in several dimensions.