Abstract
The aim of the paper is to propose an efficient and stable algorithm that is quite accurate and fast for numerical evaluation of the Fourier-Bessel transform of order , using wavelets. The philosophy behind the proposed algorithm is to replace the part of the integral by its wavelet decomposition obtained by using CAS wavelets thus representing as a Fourier-Bessel series with coefficients depending strongly on the input function . The wavelet method indicates that the approach is easy to implement and thus computationally very attractive.
1. Introduction
The Fourier-Bessel transform (also designated as Hankel transform) is a very useful tool of mathematical physics [1]. It is a very useful instrument in a wide range of physical problems which have an axial symmetry. It is particularly important in optics and two-dimensional image processing, it naturally occurs in image reconstruction from projections or from reflected pulses, and it is a useful tool in the analysis and synthesis of three-dimensional wave fields. The present development is essentially motivated by optics application. The influence of the Laplacian on a function in cylindrical coordinates is equal to the product of the squared parameter of the transformation and the transform of the function [2]
There are two types of the Hankel transform. The first one is defined on the semi-infinite interval. In this case the direct and inverse transforms of the th kind are represented as a symmetric pair. When we are dealing with problems that show circular symmetry, Hankel transforms may be very useful [3, 4]. Laplace’s partial differential equation in cylindrical coordinates can be transformed into an ordinary differential equation by using the Hankel transform. Because the Hankel transform is the two-dimensional Fourier transform of a circularly symmetric function, it plays an important role in optical data processing [5–7]. In optics, the Hankel transform appears in many contexts, not the least of which is the propagation of cylindrically symmetric laser beams. Most classical optical systems like mirrors or lenses are axially symmetric devices. Hankel transform also proved to be extremely useful in problems associated with seismology, geophysics [8, 9], electroscattering, acoustics, hydrodynamics, image processing [10], time dependent Schrodinger equation, and so forth.
Mathematical Background. The Fourier-Bessel transform may be defined by the following expression: where is the th-order Bessel function of the first kind.
In the case of the finite Hankel transform only a direct transform has an integral form. Without loss of generality its expression is (see [11]). Practical calculation of direct and inverse Hankel transform is connected with two problems. The first problem is based on the fact that not every transform in the real physical situation has analytical expression for result of inverse Hankel transform. The second one is the determination of functions as a set of their values for numerical calculations. The classical trapezoidal rule, Cotes rule, and other rules connected with the replacement of the integrand by sequence of polynomials have high accuracy if integrand is a smooth function. But (or ) is a quick oscillating function if (or ) is large. There are two general methods of the effective calculation in this area. The first is the fast Hankel transform [12]. The specification of that method is transforming the function to the logarithmical space and fast Fourier transform in that space. This method needs a smoothing of the function in log space. The second method is based on the separation of the integrand into product of slowly varying component and a rapidly oscillating Bessel function [13]. But it needs the smoothness of the slow component for its approximation by lower-order polynomials.
To overcome these difficulties, various different techniques are available in the literature. Several papers have been written on the numerical evaluation of the HT in general and the zeroth order in particular [14–24]. There are two general methods of the effective calculation in this area. The first is the fast Hankel transform [25]. The specification of that method is transforming the function to the logarithmical space and fast Fourier transform in that space. This method needs a smoothing of the function in log space. The second method is based on the separation of the integrand into product of slowly varying component and a rapidly oscillating Bessel function [26]. But it needs the smoothness of the slow component for its approximation by lower-order polynomials. From variety of algorithm, a potential user would probably find it difficult to select any one algorithm that might be best for a particular application. For an overview of these algorithms and their numerical complexity, the reader is referred to [27–31].
The organization of the paper is as follows: Section 2 gives a brief description of the CAS wavelets followed by the derivation of the algorithm in Section 3. The efficiency and stability of the algorithm are shown by applying it to four test functions with known analytical transform in Section 4. At the end, a brief conclusion and future work are given in Section 5.
2. Properties of CAS Wavelets
2.1. Wavelets and CAS Wavelets
Wavelets constitute a family of functions constructed from dilation and translation of a single function called the mother wavelets. When the dilation parameter is 2 and the translation parameter is 1 we have the following family of discrete wavelets [32]:where form a wavelet orthonormal basis for .
CAS wavelets involve four arguments , , , and , where , is assumed to be any nonnegative integer, is any integer, and is normalized time. CAS wavelets are defined as [33]whereAn efficient algorithm has been presented for the Fourier-Bessel transform.
3. Outline of Algorithm
The function representing physical fields either are zero or have an infinitely long decaying tail outside a disk of finite radius . Hence, in most practical applications either the signal has a compact support or for a given there exists such that .
Therefore, in either case, known as the finite Hankel transform (FHT), is a good approximation of the HT as given by (2). Writing in (8), we getWe may expand as follows:where .
By truncating infinite series (10) at levels and , we obtain an approximate representation for aswhere the matrices and are given by Substituting (11) in (9), we getTaking and , (13) reduces to Now, we relabel and write (14) aswhere ’s are the th-place integral in (14).
The integrals arising in (14) are evaluated by using the formulae(see [34]) and are calculated with the help of Simpson’s one-third rule, Simpson’s three-eighth rule, composite Simpson’s one-third rule, and composite Simpson’s three-eighth rule, respectively. In numerical analysis, Simpson’s rule and composite Simpson’s rule are method for numerical integration, the numerical approximation of definite integrals.
4. Numerical Results
In this section, we test the proposed algorithm (15) by evaluating the approximate Hankel transforms of 4 well-known test functions with known analytical Hankel transforms. Note that in all the examples the truncation is done at level , , and in (15). We observed that the accuracy of the method is very high even at such a low level of truncation.
Example 1. Let , ; then(obtained from [34] by putting , ), where is a Lommel function of two variables,The comparison of the approximation (dotted line) with the exact Hankel transform (solid line) is shown in Figures 1, 3, and 5 and the error in Figures 2, 4, and 6.
Simpson’s One-Third Rule. See Figures 1 and 2.
Simpson’s Three-Eighth Rule. See Figures 3 and 4.
Composite Simpson’s One-Third Rule. See Figures 5 and 6.
Composite Simpson’s Three-Eighth Rule. See Figures 7 and 8.
Example 2. The following example was solved numerically by [35].
ForWe solve the above problem by the proposed algorithm and observe that our method gives a result comparable to [35]. Note that and are indicated by (solid line) and (dotted line) in Figures 9, 11, 13, and 15 and the error is shown in Figures 10, 12, 14, and 16.
Simpson’s One-Third Rule. See Figures 9 and 10.
Simpson’s Three-Eighth Rule. See Figures 11 and 12.
Composite Simpson’s One-Third Rule. See Figures 13 and 14.
Composite Simpson’s Three-Eighth Rule. See Figures 15 and 16.
Example 3 (sombrero function). A very important and often used function is the Circ function that can be defined as [22]This function is quite common in optical problems where it is used, for instance, to represent a circular pupil of radius . The Fourier-Bessel transform of (20) is the well-known “sombrero function.”
The zeroth-order Hankel transform of is the sombrero function [29], given byThe exact and numerical transforms differ very slightly but the differences are hardly visible.
Simpson’s One-Third Rule. See Figures 17 and 18.
Simpson’s Three-Eighth Rule. See Figures 19 and 20.
Composite Simpson’s One-Third Rule. See Figures 21 and 22.
Composite Simpson’s Three-Eighth Rule. See Figures 23 and 24.
Example 4. Let , ; then,a well-known result. The pair arises in optical diffraction theory [36]. The function is the optical transfer function of an aberration-free optical system with a circular aperture, and is the corresponding spread function.
Barakat and Sandler [26] evaluated numerically using Filon quadrature philosophy but the associated error is appreciable for , whereas our method gives almost zero error in that range.
Simpson’s One-Third Rule. See Figures 25 and 26.
Simpson’s Three-Eighth Rule. See Figures 27 and 28.
Composite Simpson’s One-Third Rule. See Figures 29 and 30.
Composite Simpson’s Three-Eighth Rule. See Figures 31 and 32.
5. Conclusion
Since the basis functions used to construct the wavelets are orthogonal and have compact support, it makes them more useful and simple in actual computations. Also, since the numbers of mother wavelet’s components are restricted to one, they do not lead to the growth of complexity of calculations. Our choice of wavelets makes them more attractive in their applications in the applied physical problems as they eliminate the problems connected with the Gibbs phenomenon taking place in [30]. A good agreement between the obtained solution and some well-known results has been obtained. Four test examples are provided to show the advantage of using wavelets. This method is capable of greatly reducing the size of calculations while still maintaining high accuracy of the numerical solution.
Proposed wavelet method is very simple and attractive. The implementation of the current approach in analogy to existing methods is more convenient and the accuracy is high. The numerical example and the compared results support our claim. The difference between the exact and approximate solutions for each example was plotted graphically to determine the accuracy of numerical solutions.
5.1. Future Work
Since computational work is fully supportive of compatibility of the proposed algorithm, the same may be extended to other physical problems also. A very high level of accuracy explicitly reflects the reliability of this scheme for such problems. We would like to stress that the approximate solution includes not only time information but also frequency information due to the localization property of wavelet basis; with some change we can apply this method with the help of other wavelet bases.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to express their sincere gratitude to the editor and reviewers for their constructive comments and suggestion.