# Stable Numerical Evaluation of Finite Hankel Transforms and Their Application

Authors: Manoj P. Tripathi; B. P. Singh; Om P. Singh
Journal: International Journal of Analysis
 (2014)

---

## Abstract

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.


---

## Body

## 1. Introduction

Classically, the Hankel transformFν of order ν of a function f(r) is defined by
(1)Fνp=∫0∞rfrJνprdr,p>0.
As the Hankel transform is self-reciprocal, its inverse is given by
(2)fr=∫0∞pFνpJνprdp,
where Jν 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 Fν(p) 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. Butrf(r)Jν(pr) and pFν(p)Jν(pr) are rapidly oscillating functions for large r and p, 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 rf(r) and the latter is Jν(pr). But the error associated with Filon quadrature philosophy is appreciable for p<1. 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 ν≥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 αθi where θi is a uniform random variable in [-1,1], added to the input signal f(r). 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 ofrf(r) 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 αθi added to the input field f(r), where θi 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.


## 2. Hat Functions and Their Associated Properties

Hat functions are defined on the domain0,1. These are continuous functions with shape of hats, when plotted on two-dimensional planes. The interval 0,1 is divided into n subintervals ih,(i+1)h, i=0,1,2,…,n-1, of equal lengths h where h=1/n. The hat function’s family of first (n+1) hat functions is defined as follows [16]:
(3)ψ0t=h-th,0≤t<h,0,otherwise,ψit=t-i-1hh,(i-1)h≤t<ih,i+1h-th,ih≤t<i+1h,i=1,2,…,n-1,0,otherwise,ψnt=t-(1-h)h,1-h≤t≤1,0,otherwise.
From the definition of hat functions it is obvious that
(4)ψikh=1,i=k,0,i≠k.
The hat functions ψj(t) are continuous, linearly independent and are in L20,1.A functionf∈L2[0,1] may be approximated as
(5)ft≅∑i=0i=nfiψit=f0ψ0t+f1ψ1t+f2ψ2t+⋯+fnψnt.
The important aspect of using extended hat functions in the approximation of function f(t) lies in the fact that the coefficients fi in (5) are given by
(6)fi=fih,fori=0,1,2,…,n,whereh=1n.


## 3. Method of Solution

To derive the algorithm, we first assume that the domain space of input signalf(r) extends over a limited region 0≤r≤R. From physical point of view, this assumption is reasonable due to the fact that the input signal f(r) which represents the physical field either is zero or has an infinitely long decaying tail outside a disc of finite radius R. Therefore, in many practical applications either the input signal f(r) has a compact support or for a given ε>0 there exists a positive real R such that ∫R∞rf(r)Jν(pr)dr<ε, which is the case if f(r)=o(rλ), where λ<-3/2 as r→∞. Hence in either case, from (1), we have
(7)F^νp=∫0RrfrJνprdr.
By scaling (7) may be written as
(8)F^νp=∫01rfrJνprdr,
which is known as finite Hankel transform (FHT) [14, 15, 17]. The FHT is a good approximation of the HT given by (1). Writing rf(r)=g(r) in (8), we get
(9)F^νp=∫01grJνprdr.
Using (5) and (6), we may approximate g(r) as
(10)gr≅∑i=0ngihψir.
From (10), (9) may be written as
(11)F^νp≅∫01∑i=0ngihψirJνprdr=∑i=0ngih∫01ψirJνprdr.
Using (3), (11) may be expressed as
(12)F^νp≅g0∫0hh-rhJν(pr)dr+∑i=1n-1gih∫(i-1)hihr-i-1hhJνprdrmmmiinmm+∫ih(i+1)h(i+1)h-rhJν(pr)dr+gnh∫1-h1r-1-hhJνprdr.
It is noteworthy here that the integrals involved in (12) are evaluated by using the following formulae:
(13)∫0aJνtdt=2lim⁡N→∞∑n=0NJν+2n+1a,Reν>-1
(see [18, p. 333]),
(14)∫0at1-νJνtdt=12ν-1ν-a1-νJν-1a
(see [18, p. 333]),
(15)∫0atμJνtdt=aμν+μ+1/2ν-μ+1/2×lim⁡N→∞∑n=0Nν+2n+1Γν-μ+1/2+nΓν+μ+3/2+nmmmmmmmm×Jν+2n+1a,Reν+μ+1>0
(see [19, p. 480]).


## 4. Error Analysis

Let the R.H.S. of (10) be denoted by gn(r); that is,
(16)gnr=∑i=0ngihψir.
Now replacing g(r) by gn(r) in (9) we define an nth approximate F^ν,n(p) of the FHT F^ν(p) as follows.Definition 1.
Annth approximate finite Hankel transform of f(r), denoted by F^ν,n(p), is defined as
(17)F^ν,np=∫01gnrJνprdr=∑i=0ngih∫01ψirJνprdr.
Let εn(p) denote the absolute error between the FHT F^ν(p) and its nth approximate F^ν,n(p); then we have the following.Theorem 2.
Ifg(r) is approximated by the family of first (n+1) hat functions as given in (10), then(i)
g(jh)-gn(jh)=0, for j=0,1,2,…,n.(ii)
g(r)-gn(r)≤1/2n2g′′(jh)+O1/n3, for jh<r<(j+1)h, j=0,1,2,…,n-1;(iii)
εn(p)=F^ν(p)-F^ν,n(p)≤M/2n2+O1/n3, where g′′(jh)≤M;Proof.
(i) From (16) and (4), the value of gn(r) at jth nodal point r=jh, j=0,1,2,…,n is given by gn(jh)=∑i=0ng(ih)ψi(jh)=g(jh).
So,g(jh)-gn(jh)=0, for j=0,1,2,…,n.
(ii) Ifr lies between two consecutive integer multiples of h, that is, jh<r<(j+1)h, j=0,1,2,…,n-1, then from (16), we have
(18)gnr=∑i=0ngihψir=gjhψjr+gj+1hψj+1r=gjhj+1h-rh+gjh+hr-jhhUsing3=gjh-jhgjh+h-gjhh+rgjh+h-gjhh=gjh+r-jhgjh+h-gjhh.
As h→0, from (18), we obtain
(19)gnr≅gjh+r-jhg′jh.
By expanding gn(r) in the form of Taylor’s series, in the powers of (r-jh), we have
(20)gr=∑k=0∞r-jhkk!g(k)jh,
where g(k) denotes the kth order derivative of g(r). Using (19) and (20), the error between exact and approximate values of g(r) is given by
(21)gr-gnr=∑k=2∞r-jhkk!g(k)(jh)=r-jh22g′′jh+Or-jh3.
Since (r-jh)<h and nh=1, from (21), we get
(22)gr-gnr≤12n2g′′jh+O1n3.
(iii) The absolute errorεn(p) between exact FHT F^ν(p) and nth approximate FHT F^ν,n(p) is given by
(23)εnp=F^νp-F^ν,np=∫01grJνprdr-∫01gnrJνprdr=∑j=0n-1∫jh(j+1)hgr-gnrJνprdr.
As Jν(pr)≤1, we have
(24)εnp≤∑j=0n-1∫jh(j+1)hgr-gnrdr≤∑j=0n-1∫jhj+1h12n2g′′jh+O1n3drhhhhhhhhhhhhihhhFollowsfrom22.
If g′′(jh)≤M, then it follows that
(25)εnp≤M2n2∑j=0n-1∫jh(j+1)hdr+O1n3∑j=0n-1∫jh(j+1)hdr=M2n2+O1n3,asnh=1.


## 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 byf(r) and the noisy data function fα(r) is obtained by adding a random noise α to f(r) such that fα(ri)=f(ri)+αθi, where ri=ih, i=0,1,2,…,n, nh=30, and θi is a uniform random variable with values in [-1,1] such that max⁡fα(ri)-f(ri)≤α, 0≤i≤n. If rfα(r) is denoted by gα(r) and the approximate HT of the perturbed function is denoted by F^ν,nα(p), then from (17)(26)F^ν,nαp=∑i=0ngαih∫01ψirJνprdr=∑i=0nihfih+αθi∫01ψirJνprdr.
From (26) and (17), we have
(27)F^ν,nαp-F^ν,np=∑i=0ngαih∫01ψirJνprdrmmn-∑i=0ngih∫01ψirJνprdr=∑i=0nihαθi∫01ψirJνprdr≤∑i=0nihαθi∫01ψi(r)Jνprdr≤∑i=0nihαθi∫01ψirdr,AsJνpr≤1.
Let θM=max⁡{θi}; then substituting ψi(r) from (3), we get
(28)F^ν,nαp-F^ν,np≤θMhα×∑i=1n-1i∫(i-1)hihr-i-1hhdr+∫ih(i+1)hi+1h-rhdr+n∫1-h1r-1-hhdr∑i=1n-1i∫(i-1)hihr-i-1hhdr=θMαh2∑i=1n-1i+nh22=θMαn2h22≤α2,AsθM≤1,nh=1.
Thus we have proved that.Theorem 3.
When the input signalf(r) is corrupted with noise α, the proposed algorithm reduces the noise at least by a factor of 1/2 in the output data F^ν,n(p).


## 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 α0=0, α1=0.002, and α2=0.005 (as discussed in the beginning of Section 5). In the following examples, the errors Ej(p), j=0,1,2 are computed and their graphs are sketched for n=100, where Ej(p) = (Exact FHT F^ν(p) − approximate FHT F^ν,nαj(p) obtained from (26) with random noise αj).Further the parameterp ranges between 0 and 30 in steps of 0.2. Figures 3 and 6 depict the graph of F^ν,nαj(p)-F^ν,n(p), j=1,2, 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 Ej(p)2 involved in computations of approximate FHT with noises αj, j=0,1,2, are given in Table 1. These are calculated using the formula
(29)Ejp2=∑i=0nEjpi2n+1,
where pi is taken in steps of 0.2 in the range [0,30].Table 1
Least square errors‖Ej(p)‖2, in different examples for n=100.
‖
E
j
(
p
)
‖
2
Example1
Example2
Example3
Example4
Example5
‖
E
0
(
p
)
‖
2
7.2696
e
-
005
1.1078
e
-
016
4.9160
e
-
006
6.9741
e
-
005
6.6819
e
-
006
‖
E
1
(
p
)
‖
2
6.9980
e
-
005
1.3724
e
-
004
1.5965
e
-
004
7.3362
e
-
005
6.0667
e
-
005
‖
E
2
(
p
)
‖
2
2.2536
e
-
004
3.3312
e
-
004
3.4473
e
-
004
1.7687
e
-
004
1.3147
e
-
004Example 1.
Consider the functionf(r)=(1-r2), 0≤r≤1, given in [14, 20], for which
(30)F1(p)=πJ12(p/2)2p,0<p<∞,0,p=0.
Numerical evaluation of F1(p) has been achieved by Barakat and Sandler [20], by using Filon quadrature philosophy but again the associated error is appreciable for p<1, whereas our method gives almost zero error in that range. Equation (12) is used to obtain the nth approximate HT F^1,n(p). The comparison between exact HT F1(p) and approximate HT F^1,n(p) is shown in Figure 1. In Figure 2, the errors E0(p), E1(p), and E2(p) for n=100 are plotted. In Figure 3, F^1,nα(p)-F^1,n(p) is shown for noises α=0.002, 0.005. This figure is in conformity with Theorem 3.Figure 1
Exact HTF1(p) (solid line) and approximate HT F^1,n(p) (dotted line) for n=100, Example 1.Figure 2
E
0
(
p
), E1(p), E2(p) for n=100, Example 1.Figure 3
Plot ofF^1,nα(p)-F^1,n(p) for α=0.002 and α=0.005 with n=100, Example 1.Example 2 (Sombrero function).
A very important and often used function is theCirc function that can be defined as
(31)Circra=1,r≤a,0,r>a.
The zeroth order HT of Circ(r/a) is the Sombrero function [13, 21], that will be written as S0(p) with the following analytical expression:
(32)S0p=a2J1apap.
We use (12) to obtain the approximation for the FHT F^0,n(p). In Figure 4, we sketch the error E0(p) between exact HT S0(p) and approximate FHT F^0,n(p) of Circ function (without noise) for a=1 and n=10. It is evident that, even for such a small value of n, the error is appreciably small. Figure 5 compares the errors E1(p) and E2(p) for n=100. Figure 6 shows the plot of F^0,nα(p)-F^0,n(p) for α=0.002 and α=0.005 with n=100.Figure 4
E
0
(
p
) for n=10, Example 2.Figure 5
E
1
(
p
) and E2(p) for n=100, Example 2.Figure 6
Plot ofF^0,nα(p)-F^0,n(p) for α=0.002 and α=0.005 with n=100, Example 2.Example 3.
Consider the functionf(r)=2/πarccos(r)-r(1-r2), 0≤r≤1, given in [14, 22] for which zeroth order exact HT is
(33)F0(p)=2J12(p/2)p2,0≤p<∞.
It is a well-known pair which arises in optical diffraction theory [23]. The function f(r) is known as optical transfer function of an aberration-free optical system with a circular aperture, and F0(p) is the corresponding spread function. Barakat and Parshall [22] evaluated F0(p) using Filon quadrature philosophy but the associated error is again appreciable for p<1, whereas our method gives significantly small error in that range. In Figure 7, the error E0(p) is shown. Figure 8 is the comparison between E1(p) and E2(p) for n=100.Figure 7
E
0
(
p
) for n=100 in Example 3.Figure 8
E
1
(
p
) and E2(p) for n=100 in Example 3.Example 4.
Consider the function [24], f(r)=rν(a2-r2)μH(a-r), Re(ν)>-1.
Re
(
μ
)
>
-
1, 0<r<a, where H(r) is a unit step function defined by
(34)Hr=1,r≥00,r<0.
The νth order HT of f(r) is given by
(35)Fνp=2νaμ+ν+1p-μ-1μ+1Jν+μ+1ap.
Using the method of solution, developed in Section 3, the approximate HT F^ν,n(p) has been calculated for ν=3/2, μ=1/2, and a=1. The comparison between F^3/2,n(p) and exact FHT F^3/2(p) is shown in Figure 9. Further in Figure 10, a comparison between the different errors E0(p), E1(p), and E2(p) is shown, for n=100.Figure 9
Exact HTF3/2(p) (solid line) and approximate HT F^3/2,n(p) (dotted line) for μ=1/2 and n=100, Example 4.Figure 10
E
0
(
p
), E1(p), E2(p) for n=100 in Example 4.Example 5.
In this example we choose as a test function the generalized version of the top-hat function, given asf(r)=rν[H(r)-H(r-a)], a>0, and H(r) is the unit step function defined by (34). Then,
(36)Fνp=Jν+1pp.
In [9], authors took a=1 and ν=4 for numerical calculations. We take a=1, ν=5 and observe that the associated errors with and without random noises are quite small. The error E0(p) (without noise) is shown in Figure 11 for n=100. In Figure 12, errors E1(p) and E2(p) for n=100 are shown.Figure 11
E
0
(
p
) for n=100 in Example 5.Figure 12
Comparison betweenE1p and E2p for n=100 in Example 5.


## 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 byf(r) and radiation takes place at the surface into the surrounding medium maintained at zero temperature. We seek a function u(r,t), where r is radius and t is time, (u does not depend upon θ and t) in the following application. The mathematical model of this problem is the diffusion equation in polar coordinates given by
(37)∂2ur,t∂2r+1r∂ur,t∂r-ν2r2ur,t=∂ur,t∂t,0<r<1,0<t<∞
with the following initial and boundary conditions:(i)
ast→0+, u(r,t)→f(r)=rνH(r)-H(r-1), where 0≤r≤1 and H(r) is unit step function given by (34),(ii)
asr→1-, ∂u/∂r+Gu→0 for each fixed t>0 and G>0.When u(r,t) denotes the temperature within the cylinder, G>0 means that heat is being radiated away from the surface of the cylinder.LetΩν,r denote the Bessel differential operator ∂2/∂2r+1/r∂/∂r-ν2/r2. Then differential equation (37) can be written as
(38)Ων,ru=∂u∂t.
Applying the νth order finite Hankel transform operator to both sides of the differential equation (38) and using (8) we get
(39)∫01rΩν,ruJν(βν,sr)dr=∫01r∂u∂tJν(βν,sr)dr,
where βν,s is the sth positive root of equation [24]
(40)GJνz+zJν′z=0,ν≥-12.
The L.H.S. of (39) may be written as
(41)∫01r∂2u∂2r+1r∂u∂r-ν2r2uJν(βν,sr)dr=∫01∂2u∂2rrJνβν,srdr+∫01∂u∂rJνβν,srdr-ν2∫01urJνβν,srdr=lim⁡r→1-∂ur,t∂rJνβν,s-βν,s∫01r∂u∂rJν′βν,srdr-ν2∫01urJνβν,srdr=-GJνβν,slim⁡r→1-ur,t-βν,s∫01∂u∂rrJν′βν,srdr-ν2∫01urJνβν,srdrAsr⟶1-,∂u∂r+Gu⟶0.
Since βν,s is a root of GJν(z)+zJν′(z)=0, we have GJν(βν,s)=-βν,sJν′(βν,s), and thus
(42)∫01r∂2u∂2r+1r∂u∂r-ν2r2uJν(βν,sr)dr=βν,sJν′βν,slim⁡r→1-ur,t-βν,sJν′βν,slim⁡r→1-ur,t+βν,s∫01uJν′βν,sr+βν,srJν′′βν,srdr-ν2∫01urJνβν,srdr.
Using the fact that Jν(x) is a solution of x2d2y/dx2+xdy/dx+(x2-ν2)y=0, (42) may be simplified as
(43)∫01r∂2u∂2r+1r∂u∂r-ν2r2uJν(βν,sr)dr=-βν,s2∫01ur,trJνβν,srdr.
From (43) and (39), we have
(44)-βν,s2∫01ur,trJνβν,srdr=∂∂t∫01ru(r,t)Jν(βν,sr)dr.
If FHT ∫01u(r,t)rJν(βν,sr)dr of u(r,t) is denoted by U-(βν,s,t), then (44) may be written as
(45)∂U-βν,s,t∂t=-βν,s2U-βν,s,t,
whose solution is given by
(46)U-βν,s,t=Aβν,se-βν,s2t.
The constant of integration A(βν,s) in (46) is determined by initial condition. Thus(47)Aβν,s=F^νβν,s=∫01frrJν(βν,sr)dr.
Hence from (47), (46) becomes
(48)U-βν,s,t=F^νβν,se-βν,s2t.
Therefore, by inversion theorem [25], we have
(49)ur,t=lim⁡N→∞∑s=1N2βν,s2F^νβν,se-βν,s2tJνβν,sr(βν,s2-ν2)Jν2(βν,s)+βν,s2Jν′2βν,s=lim⁡N→∞∑s=1N2βν,s2F^νβν,se-βν,s2tJνβν,srβν,s2-ν2Jν2βν,s+G2Jν2βν,s,mmmmmmAsGJνβν,s=-βν,sJν′βν,s.
So
(50)ur,t=lim⁡N→∞∑s=1N2βν,s2F^νβν,se-βν,s2tβν,s2-ν2+G2Jνβν,srJν2βν,s.
We want to prove that u(r,t), 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]:(51)F^νβν,s=O1βν,s3/2βν,s≈πs+14,Jν2βν,s≈2πβν,s,ass⟶∞.
Using the above estimates, we see that series (50) and the series obtained by applying Ων,r and ∂/∂t separately under the summation sign of (50) converge uniformly on 0<r<1 and t>0. Hence by applying Ων,r-∂/∂t on (50) we get
(52)Ων,r-∂∂tu(r,t)=lim⁡N→∞∑s=1N2βν,s4F^ν(βν,s)e-βν,s2t(βν,s2-ν2+G2)mmni·∑s=1N2βν,s4F^ν(βν,s)e-βν,s2t(βν,s2-ν2+G2)βν,s2r2Jν′′βν,sr+βν,srJν′βν,srmmmmnm+βν,s2r2-ν2Jνβν,srmmminm×βν,s2r2Jν2βν,s-1∑s=1N2βν,s4F^ν(βν,s)e-βν,s2t(βν,s2-ν2+G2)=0.
Hence we see that (50) satisfies the differential equation (37).Let us verify the boundary condition (ii); we have(53)lim⁡r→1-∂ur,t∂r+Gur,t=lim⁡r→1-∑s=1N2βν,s2F^νβν,se-βν,s2tβν,s2-ν2+G2Jν2βν,smmmmm×βν,sJν′βν,sr+GJνβν,sr∑s=1N2βν,s2F^νβν,se-βν,s2tβν,s2-ν2+G2Jν2βν,s,
and since the convergence is uniform, we can take the lim⁡r→1- inside the summation sign and arrive at the conclusion, as βν,s’s are the roots of the equation GJν(z)+zJν′(z)=0.The initial condition (i) is already taken care of as we evaluated the constantA(βν,s) 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 ν=5 and G=1. We have used MATLAB routine “fzero” to find zeros of (40) and we have drawn Figures 13–18, by truncating series (50) at N=30. While evaluating the solution u(r,t) from (50), we have evaluated F^5(β5,s) first from its analytical expression given by (36) and denoted the solution thus obtained by u(r,t) in Figures 13–18 and then evaluating nth approximates FHT F^5,n(β5,s) for n=100, using our proposed algorithm for evaluation of the FHT as given by (17). This solution is denoted by ua(r,t) in the above mentioned figures.Figure 13
The initial condition functionf(r) (solid red line) and uar,0  (=lim⁡t→0+ua(r,t)) (*blue).Figure 14
Error betweenf(r) and ua(r,0).Figure 15
The various profiles of the solutionsua(r,t) at fixed times.Figure 16
The absolute error betweenua(r,t) and u(r,t), for different t.Figure 17
The solutionu(r,t) for 0<t<1 and 0<r<1.Figure 18
The approximate solutionua(r,t) using nth approximate FHT F^5,n(β5,s), for 0<t<1 and 0<r<1.Figure13 compares the given initial condition f(r) with ua(r,t) as t→0+ and Figure 14 shows the corresponding error E(r)=ua(r,0)-f(r). Figure 15 depicts the various profiles of ua(r,t) at fixed times tj=0, 1/100, 1/50, 1/25, where the various profiles are denoted by u0,  u1,  u2, and u3. The absolute errors between ua(r,t) and u(r,t), at fixed times tj=0, 1/100, 1/50, 1/25, are shown in Figure 16. As the maximum possible error occurs in the neighborhood of 0 and 0.001, we have restricted t in (0,1] in Figures 17 and 18, representing u(r,t) and ua(r,t), 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 signalrf(r) 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.



---

*Source: 670562-2014-11-13.xml*