Stable Numerical Evaluation of Finite Hankel Transforms and Their Application

A new stable algorithm, based on hat functions for numerical evaluation of Hankel transform of order , is proposed in this paper. The hat basis functions are used as a basis to expand a part of the integrand, , appearing in the Hankel transform integral. This leads to a very simple, efficient, and stable algorithm for the numerical evaluation of Hankel transform. The novelty of our paper is that we give error and stability analysis of the algorithm and corroborate our theoretical findings by various numerical experiments. Finally, an application of the proposed algorithm is given for solving the heat equation in an infinite cylinder with a radiation condition. 1. Introduction Classically, the Hankel transform of order of a function is defined by As the Hankel transform is self-reciprocal, its inverse is given by where is the th order Bessel function of first kind. This form of Hankel transform (HT) has the advantage of reducing to the Fourier sine or cosine transform when . The Hankel transform arises naturally in the discussion of problems posed in cylindrical coordinates and hence, as a result of separation of variables, involving Bessel functions. So, it has found wide range of applications related to the problems in mathematical physics possessing axial symmetry [1]. But analytical evaluations of are rare, so numerical methods are important. The usual classical methods like Trapezoidal rule, Cotes rule, and so forth connected with replacing the integrand by sequence of polynomials have high accuracy if integrand is smooth. But and are rapidly oscillating functions for large and , respectively. To overcome these difficulties, two different techniques are available in the literature. The first is the fast Hankel transform as proposed by Siegman [2]. Here, by substitution and scaling, the problem is transformed in the space of the logarithmic coordinates and the fast Fourier transform in that space. But it involves the conventional errors arising when a nonperiodic function is replaced by its periodic extension. Moreover, it is sensitive to the smoothness of function in that space. The second method is based on the use of Filon quadrature philosophy [3]. In Filon quadrature philosophy, the integrand is separated into the product of an (assumed) slowly varying component and a rapidly oscillating component. In the context of Hankel transform, the former is and the latter is . But the error associated with Filon quadrature philosophy is appreciable for . There are several extrapolation methods developed in the eighties. In particular, the papers by Levin