Abstract
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 and Sidi [4–7] are relevant. Levin and Sidi in 1981 [4] developed very general class of extrapolation methods for various types of infinite integrals, whether oscillatory or not. Sidi in [5–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 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 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 , 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 . 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 , 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.
2. Hat Functions and Their Associated Properties
Hat functions are defined on the domain . These are continuous functions with shape of hats, when plotted on two-dimensional planes. The interval is divided into subintervals , , of equal lengths where . The hat function’s family of first hat functions is defined as follows [16]: From the definition of hat functions it is obvious that The hat functions are continuous, linearly independent and are in .
A function may be approximated as The important aspect of using extended hat functions in the approximation of function lies in the fact that the coefficients in (5) are given by
3. Method of Solution
To derive the algorithm, we first assume that the domain space of input signal extends over a limited region . 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 there exists a positive real such that , which is the case if , where 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). Writing in (8), we get Using (5) and (6), we may approximate as From (10), (9) may be written as Using (3), (11) may be expressed as It is noteworthy here that the integrals involved in (12) are evaluated by using the following formulae: (see [18, p. 333]), (see [18, p. 333]), (see [19, p. 480]).
4. Error Analysis
Let the R.H.S. of (10) be denoted by ; that is, Now replacing by in (9) we define an approximate of the FHT as follows.
Definition 1. An th approximate finite Hankel transform of , denoted by , is defined as Let denote the absolute error between the FHT and its th approximate ; then we have the following.
Theorem 2. If is approximated by the family of first hat functions as given in (10), then(i), for .(ii), for , ;(iii), where ;
Proof. (i) From (16) and (4), the value of at th nodal point , is given by .
So, , for .
(ii) If lies between two consecutive integer multiples of , that is, , , then from (16), we have
As , from (18), we obtain
By expanding in the form of Taylor’s series, in the powers of , we have
where denotes the th order derivative of . Using (19) and (20), the error between exact and approximate values of is given by
Since and , from (21), we get
(iii) The absolute error between exact FHT and th approximate FHT is given by
As , we have
If , then it follows that
5. Stability Analysis
In this section, the stability of the proposed algorithm is analyzed under the influence of noise. In what follows, the exact data function is denoted by and the noisy data function is obtained by adding a random noise to such that , where , , , and is a uniform random variable with values in such that , . If is denoted by and the approximate HT of the perturbed function is denoted by , then from (17) From (26) and (17), we have Let ; then substituting from (3), we get 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 in the output data .
6. Numerical Illustrations
The test problems included in this paper are solved with and without random perturbations (noises) to illustrate the efficiency and stability of proposed algorithm by choosing three different values of noise as , , and (as discussed in the beginning of Section 5). In the following examples, the errors , are computed and their graphs are sketched for , where = (Exact FHT − approximate FHT obtained from (26) with random noise ).
Further the parameter ranges between and in steps of . Figures 3 and 6 depict the graph of , , for different test functions in Examples 1 and 2, respectively. These graphs are in conformity with Theorem 3. For all the illustrations, the computations are done in MATLAB 7.0.1. For different examples, the least square errors involved in computations of approximate FHT with noises , , are given in Table 1. These are calculated using the formula where is taken in steps of in the range .
Example 1. Consider the function , , given in [14, 20], for which Numerical evaluation of has been achieved by Barakat and Sandler [20], by using Filon quadrature philosophy but again the associated error is appreciable for , whereas our method gives almost zero error in that range. Equation (12) is used to obtain the th approximate HT . The comparison between exact HT and approximate HT is shown in Figure 1. In Figure 2, the errors , , and for are plotted. In Figure 3, is shown for noises , . This figure is in conformity with Theorem 3.
Example 2 (Sombrero function). A very important and often used function is the function that can be defined as The zeroth order HT of is the Sombrero function [13, 21], that will be written as with the following analytical expression: We use (12) to obtain the approximation for the FHT . In Figure 4, we sketch the error between exact HT and approximate FHT of function (without noise) for and . It is evident that, even for such a small value of , the error is appreciably small. Figure 5 compares the errors and for . Figure 6 shows the plot of for and with .
Example 3. Consider the function , , 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 is the corresponding spread function. Barakat and Parshall [22] evaluated using Filon quadrature philosophy but the associated error is again appreciable for , whereas our method gives significantly small error in that range. In Figure 7, the error is shown. Figure 8 is the comparison between and for .
Example 4. Consider the function [24], , .
, , 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 has been calculated for , , and . The comparison between and exact FHT is shown in Figure 9. Further in Figure 10, a comparison between the different errors , , and is shown, for .
Example 5. In this example we choose as a test function the generalized version of the top-hat function, given as , , and is the unit step function defined by (34). Then, In [9], authors took and for numerical calculations. We take , and observe that the associated errors with and without random noises are quite small. The error (without noise) is shown in Figure 11 for . In Figure 12, errors and for are shown.
7. 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:(i)as , , where and is unit step function given by (34),(ii)as , for each fixed and .When denotes the temperature within the cylinder, means that heat is being radiated away from the surface of the cylinder.
Let denote the Bessel differential operator . Then differential equation (37) can be written as Applying the th order finite Hankel transform operator to both sides of the differential equation (38) and using (8) we get where is the th positive root of equation [24] The L.H.S. of (39) may be written as Since is a root of , we have , and thus Using the fact that is a solution of , (42) may be simplified as From (43) and (39), we have If FHT of is denoted by , then (44) may be written as whose solution is given by The constant of integration in (46) is determined by initial condition. Thus Hence from (47), (46) becomes Therefore, by inversion theorem [25], we have So We want to prove that , given by (50), is truly a solution of (37) that satisfies the given initial and boundary conditions. To achieve this, we need the following well known estimates [17]: Using the above estimates, we see that series (50) and the series obtained by applying and separately under the summation sign of (50) converge uniformly on and . Hence by applying on (50) we get Hence we see that (50) satisfies the differential equation (37).
Let us verify the boundary condition (ii); we have and since the convergence is uniform, we can take the inside the summation sign and arrive at the conclusion, as ’s are the roots of the equation .
The initial condition (i) is already taken care of as we evaluated the constant by using it. Through Figures 13–18, we establish the accuracy of the proposed method. For numerical solution of differential equation (37), we have taken and . We have used MATLAB routine “fzero” to find zeros of (40) and we have drawn Figures 13–18, by truncating series (50) at . While evaluating the solution from (50), we have evaluated first from its analytical expression given by (36) and denoted the solution thus obtained by in Figures 13–18 and then evaluating approximates FHT for , using our proposed algorithm for evaluation of the FHT as given by (17). This solution is denoted by in the above mentioned figures.
Figure 13 compares the given initial condition with as and Figure 14 shows the corresponding error . Figure 15 depicts the various profiles of at fixed times , , , , where the various profiles are denoted by , , , and . The absolute errors between and , at fixed times , , , , are shown in Figure 16. As the maximum possible error occurs in the neighborhood of and , we have restricted in in Figures 17 and 18, representing and , respectively, and note that they are in good agreement in the range.
8. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.