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 ] > −1, is proposed in this paper. The hat basis functions are used as a basis to expand a part of the integrand, rf(r), 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.


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 ] = ±1/2.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  < 1.There are several extrapolation methods developed in the eighties.In particular, the papers by Levin and Sidi [4][5][6][7] are relevant.Levin and Sidi in 1981 [4] developed very 2 International Journal of Analysis general class of extrapolation methods for various types of infinite integrals, whether oscillatory or not.Sidi in [5][6][7] emphasized oscillatory integral involving Fourier and Hankel Transforms, among other more general cases.In addition the subject is dealt with in great detail in a book by Sidi [8,Chapter 5].
Later in 2004, Guizar-Sicairos and Gutierrez-Vega [9] obtained a powerful scheme to calculate the HT of order ] ≥ 0 by extending the zero order HT algorithm of Yu et al. [10] to higher orders.Their algorithm is based on the orthogonality properties of Bessel functions.Postnikov [11] proposed, for the first time, a novel and powerful method for computing zero and first order HT by using Haar wavelets.Refining the idea of Postnikov [11], Singh et al. [12,13] obtained three efficient algorithms for numerical evaluation of HT of order ] > −1 using linear Legendre multiwavelets, Legendre wavelets, and rationalized Haar wavelets which were shown to be superior to the other mentioned algorithms.Recently Singh et al. [14] and Pandey et al. [15] tested their proposed algorithms on some examples for their stability with respect to noise   where   is a uniform random variable in [−1, 1], added to the input signal ().But none of the above mentioned algorithms were analyzed theoretically for stability and no error analysis was provided.
The lack of error and stability analysis in the papers cited above motivated us for the present work.In the present paper, we use hat basis functions for approximation of () and replace it by its approximation in (1), thereby getting an efficient and stable algorithm for the numerical evaluation of the HT of order ] > −1.Error and stability analysis of the algorithm are given in Sections 4 and 5, respectively, and further corroborated by the numerical experiments performed on the test functions in Section 6. Test functions with known analytic HT are used with random noise term   added to the input field (), where   is a uniform random variable with values in [−1, 1], to illustrate the stability and efficiency of the proposed algorithm.In Section 7, the algorithm is applied for solving the heat equation in an infinite cylinder with a radiation condition.

Method of Solution
To derive the algorithm, we first assume that the domain space of input signal () extends over a limited region 0 ≤  ≤ .From physical point of view, this assumption is reasonable due to the fact that the input signal () which represents the physical field either is zero or has an infinitely long decaying tail outside a disc of finite radius .Therefore, in many practical applications either the input signal () has a compact support or for a given  > 0 there exists a positive real  such that | ∫ ∞  () ] ()| < , which is the case if () = (  ), where  < −3/2 as  → ∞.Hence in either case, from (1), we have By scaling (7) may be written as which is known as finite Hankel transform (FHT) [14,15,17].
The FHT is a good approximation of the HT given by (1).

Error Analysis
Let the R.H.S. of (10) be denoted by   (); that is, Now replacing () by   () in ( 9) we define an th approximate F], () of the FHT F] () as follows.
Definition 1.An th approximate finite Hankel transform of (), denoted by F], (), is defined as Let   () denote the absolute error between the FHT F] () and its th approximate F], (); then we have the following.

Stability Analysis
In this section, the stability of the proposed algorithm is analyzed under the influence of noise , (As   ≤ 1, ℎ = 1) . ( Thus we have proved that. Theorem 3. When the input signal () is corrupted with noise , the proposed algorithm reduces the noise at least by a factor of 1/2 in the output data F], ().
where   is taken in steps of 0.2 in the range [0, 30].
Example 1.Consider the function () = √(1 −  2 ), 0 ≤  ≤ 1, given in [14,20], for which Numerical evaluation of  1 () has been achieved by Barakat and Sandler [20], by using Filon quadrature philosophy but again the associated error is appreciable for  < 1, whereas our method gives almost zero error in that range.Equation ( 12) is used to obtain the th approximate HT F1, ().The comparison between exact HT  1 () and approximate HT F1, () is shown in Figure 1.In Figure 2, the errors  0 (),  The zeroth order HT of Circ(/) is the Sombrero function [13,21], that will be written as  0 () with the following analytical expression: We use (12) to obtain the approximation for the FHT F0, ().
In Figure 4, we sketch the error  0 () between exact HT  0 () and approximate FHT F0, () of Circ function (without noise) for  = 1 and  = 10.It is evident that, even for such a small value of , the error is appreciably small.Figure 5 compares the errors  1 () and  2 () for  = 100.)], 0 ≤  ≤ 1, given in [14,22] for which zeroth order exact HT is It is a well-known pair which arises in optical diffraction theory [23].The function () is known as optical transfer function of an aberration-free optical system with a circular aperture, and  0 () is the corresponding spread function.Barakat and Parshall [22] evaluated  0 () using Filon quadrature philosophy but the associated error is again appreciable for  < 1, whereas our method gives significantly small error in that range.In Figure 7, the error  0 () is shown.Figure 8 is the comparison between  1 () and  2 () for  = 100.
Example 4. Consider the function [24], Re() > −1, 0 <  < , where () is a unit step function defined by The ]th order HT of () is given by Using the method of solution, developed in Section 3, the approximate HT F], () has been calculated for ] = 3/2, In [9], authors took  = 1 and ] = 4 for numerical calculations.We take  = 1, ] = 5 and observe that the associated errors with and without random noises are quite small.The error  0 () (without noise) is shown in Figure 11   for  = 100.In Figure 12, errors  1 () and  2 () for  = 100 are shown.

Applications
In this section we solve heat equation in cylindrical coordinates inside an infinitely long cylinder of radius unity, by using the theory of finite Hankel transform (FHT) developed in the preceding pages.The initial temperature is given by () and radiation takes place at the surface into the surrounding medium maintained at zero temperature.We seek a function (, ), where  is radius and  is time, ( does not depend upon  and ) in the following application.The mathematical model of this problem is the diffusion equation in polar coordinates given by with the following initial and boundary conditions: where 0 ≤  ≤ 1 and () is unit step function given by (34), (ii) as  → 1 − , (/) +  → 0 for each fixed  > 0 and  > 0.
When (, ) denotes the temperature within the cylinder,  > 0 means that heat is being radiated away from the surface of the cylinder.

Conclusions
A new stable and efficient algorithm based on hat functions for the numerical evaluation of Hankel transform is proposed and analyzed.Choosing hat basis functions to expand the input signal () makes our algorithm attractive in their application in the applied physical problems as they eliminate the problems connected with the Gibbs phenomenon taking place in [9,11].We have given for the first time (as per our information) error and stability analysis which was one of the main drawbacks of the earlier algorithms and by various numerical experiments we have corroborated our theoretical findings.Stability with respect to the data is restored and excellent accuracy is obtained even for small sample interval and high noise levels in the data.From the various figures it is obvious that the algorithm is consistent and does not depend on the particular choice of the input signal.The accuracy and simplicity of the algorithm provide an edge over the many other algorithms.Finally, an application of the proposed algorithm is given for solving the heat equation in an infinite cylinder with a radiation condition.