Is the fast Hankel transform faster than quadrature?

Abstract

The fast Hankel transform (FHT) implemented with digital filters has been the algorithm of choice in EM geophysics for a few decades. However, other disciplines have predominantly relied on methods that break up the Hankel transform integral into a sum of partial integrals that are each evaluated with quadrature. The convergence of the partial sums is then accelerated through a nonlinear sequence transformation. While such a method was proposed for geophysics nearly three decades ago, it was demonstrated to be much slower than the FHT. This work revisits this problem by presenting a new algorithm named quadrature-with-extrapolation (QWE). The QWE method recasts the quadrature sum into a form conceptually similar to the FHT approach by using a fixed-point quadrature rule. The sum of partial integrals is efficiently accelerated using the Shanks transformation computed with Wynn’s [Formula: see text] algorithm. A Matlab implementation of the QWE algorithm is compared with the FHT method for accuracy and speed on a suite of relevant modeling problems including frequency-domain controlled-source EM, time-domain EM, and a large-loop magnetic source problem. Surprisingly, the QWE method is faster than the FHT for all three problems. However, when the integral needs to be evaluated at many offsets and the lagged convolution variant of the FHT is applicable, the FHT is significantly faster than the QWE method. For divergent integrals such as those encountered in the large loop problem, the QWE method can provide an accurate answer when the FHT method fails.