IMA Journal of Numerical Analysis
2014, October
A set of accurate quadrature rules applicable to a class of integrable func- tions with isolated singularities is designed and analyzed theoretically in one and two dimensions. These quadrature rules are based on the trapezoidal rule with corrected quadrature weights for points in the vicinity of the singularity. To compute the cor- rection weights small size ill-conditioned systems have to be solved. The convergence of the correction weights is accelerated by the use of compactly supported functions that annihilate boundary errors. Convergence proofs with error estimates for the resulting quadrature rules are given in both one and two dimensions. The tabulated weights are specific for the singularities under consideration, but the methodology extends to a large class of functions with integrable isolated singularities. Further- more in one dimension we have obtained a closed form expression based on which the modified weights can be computed directly.
We have constructed and also proved the accuracy of high-order quadrature rules that handle singular functions of the type $s(x)=|x|^{\gamma}$ with $\gamma>-1$ in one dimension and $\gamma=-1$ in two dimensions. The quadrature rules are based on the well-known trapezoidal rule, but with modified weights close to the singularity. They can be applied in a straight-forward manner by using the modified weights which have been computed and tabulated in this paper. In the construction of the quadrature rules we used compactly supported functions to annihilate boundary errors and accelerate the convergence of the correction weights. If the quadrature rule is to be applied to functions that are {neither} periodic nor compactly supported within the domain, it should be combined with boundary corrections for the trapezoidal rule of sufficiently high order to exhibit the full convergence order it has been designed to have. The correction operator associated with the singularity remains the same whether or not boundary corrections are necessary. The modified weights are obtained as the solution of an ill-conditioned linear system of equations. We have described how this system is set-up and discussed its ill-conditioning. The convergence rate of the weights depends on the \textit{flatness} of the compactly supported function used to annihilate boundary errors, i.e. on the number of derivatives that vanish at the singular point. We illustrate this property in numerical examples and prove mathematically its validity in the one-dimensional case. With an improved convergence rate for the weights, a larger $h$ can be used, which to some extent alleviates the ill-conditioning problem. To be able to compute the weights with $16$ correct digits, a multi-precision library is used. In one-dimension, we were able to find an analytical expression for the right hand side of the system such that converged weights can be computed directly (see the Appendix). The methodology for setting up the system for the weights can be extended also to other singularities, e.g. the fundamental solution of the Stokes equations, \cite{marin_2011}, which is not radially symmetric in all its tensorial components as the functions $|x|^{\gamma}$, and more complicated symmetries must be taken into account. The integration of $1/|\vect x|$ in $\mathbb R^2$ as considered here, can also be viewed as an integration of the fundamental solution of Laplace over a flat surface in $\mathbb R^3$ . This relates to the discretization of boundary integral formulations, where integrals are to be evaluated over the boundaries of the domain, may it be an outer boundary, the surface of a scattering object in an electromagnetic application or the surface of an immersed particle in Stokes flow.