High-order quadratures for singular functions and their applications to the evaluation of Fourier and Hankel transforms

Evaluation of integrals plays an important role in the numerical solution of many problems in physics, engineering, and other areas; hence the interest in the study and design of rapidly convergent quadrature rules. For example, the discrete Fourier transform (DFT) is a popular numerical tool since it provides an excellent approximation to the continuous Fourier transform of a periodic function. The reason for the accuracy of the approximation is that the trapezoidal rule, which the DFT implements, is superalgebraically convergent for periodic functions. Numerical quadrature schemes for smooth functions are a well understood subject (see, (11), (3), (28)). When a function is singular, the need for a high order quadrature scheme is met by Gaussian quadrature. Unfortunately, the nodes at which a function is tabulated are often non-Gaussian (such as equispaced nodes, Chebyshev nodes, etc.), for which existing quadrature schemes are inadequate. In this thesis we introduce a class of quadrature formulae applicable to both non-singular and singular functions, generalizing the classical end-point corrected trapezoidal quadrature rules. While the standard end-point corrected trapezoidal rules are usually derived by means of the Euler Maclaurin formula, their generalizations are obtained as solutions of certain systems of linear algebraic equations. A procedure is developed for the construction of very high-order quadrature rules, applicable to functions with apriori specified singularities, and relaxing the requirements on the distribution of nodes. We also present two applications based on these high-order quadratures. An algorithm is developed for the rapid evaluation of the Fourier transform of functions with singularities. The algorithm is based on a combination of the fast Fourier transform (FFT) with a quadrature scheme tailored to the singularity. A related algorithm for the fast Hankel transform is also presented. The algorithm decomposes the Hankel transform into a product of two integral operators, the first of which is evaluated rapidly by a combination of the fast cosine transform with a quadrature formula of the type developed in this thesis. The second operator is evaluated rapidly by a combination of a version of the fast multipole method with yet another quadrature formula derived in this thesis. All calculations are performed to full double precision accuracy. Numerical experiments are presented demonstrating the practical usefulness and efficiency of all the algorithms developed. Tables of quadrature weights are included for singularities of the form $s({x})=\vert{x}\vert\sp\lambda$ for a variety of values of $\lambda,$ and s(x) = log$(\vert{x}\vert).$