A Numerical Algorithm for Arbitrary Real-Order Hankel Transform

The Hankel transform is widely used to solve various engineering and physics problems, such as the representation of electromagnetic field components in the medium, the representation of dynamic stress intensity factors, vibration of axisymmetric infinite membrane and displacement intensity factors which all involve this type of integration. However, traditional numerical integration algorithms cannot be used due to the high oscillation characteristics of the Bessel function, so it is particularly important to propose a high precision and efficient numerical algorithm for calculating the integral of high oscillation. In this paper, the improved Gaver-Stehfest (G-S) inverse Laplace transform method for arbitrary real-order Bessel function integration is presented by using the asymptotic characteristics of the Bessel function and the accumulation of integration, and the optimized G-S coefficients are given. The effectiveness of the algorithm is verified by numerical examples. Compared with the linear transformation accelerated convergence algorithm, it shows that the G-S inverse Laplace transform method is suitable for arbitrary real order Hankel transform, and the time consumption is relatively stable and short, which provides a reliable calculation method for the study of electromagnetic mechanics, wave propagation, and fracture dynamics.


Introduction
The integration of high oscillation function generally exists in the application fields of aeronautics, seismic imaging, electromagnetic mechanics, and so on. The integral with Bessel function is widely used in geological exploration, nuclear magnetic physics, hydrodynamics, electromagnetic response, signal processing, acoustic scattering, droplet wetting, and fracture mechanics [1][2][3][4][5][6][7] . This kind of integral is becoming one of the core problems of high oscillation integrals, thus its efficient numerical solution has been a core research topic. In addition, higher-order Bessel beams, higher-order Bessel vortex electromagnetic wave, application of multiple Bessel Gaussian beams in turbulent transmission and communication, higher-order Hankel transform and its application in beams transmission involve higher-order Bessel integration [8][9][10][11] . Therefore, it is meaningful and valuable in engineering practice to study efficient algorithms of arbitrary real-order Hankel transform.
The method of dividing an oscillatory integral at its zeros, forming a sequence of partial sums, and using extrapolation to accelerate convergence have been found to be the most efficient technique available [12][13][14][15] . Chave [13] used the accumulation of integration interval to write the integral to be calculated as the sum of the integrals between the intercell, and the integral at each intercell was written in the form of continued fractions to speed up the convergence. The continued fraction method is only slightly effective for rapidly convergent integrals, but for slowly convergent integrals the summation behavior is dramatic. Lucas et al [15] compared various extrapolation techniques as well as the choice of endpoints in dividing the integral, and established the efficient method for evaluating infinite integrals including Bessel functions of arbitrary order. However, this method requires that the non-oscillating term must be monotonic, and its accuracy is related to the segmentation interval, which is not efficient. Two zero-order Hankel transform filters (61-point and 120-point operators) and two first-order Hankel transform filters (47-point and 140-point operators) have been designed [16] . For these transforms, the error of new filters is much lower compared with all other known filters. However, the numerical approximation algorithm for arbitrary order Hankel transform is not given. Filon [17] proposed a Filon-type method to solve the high oscillation function integral in a finite interval. However, it is very difficult to calculate the integral of the highly oscillating function in the infinite interval with the Filon-type method. Yu [18] modified the Filon-type method based on special functions and gave a high-precision solution of strictions for ( ) f  and its derivative are required in the algorithm. Chen [19] divided the infinite interval into two intervals for discussion, as the Filon-type method was used on the finite interval, and the asymptotic method was constructed on the infinite interval, which overcame the defects of the Filon-type method, where oscillation factor 1 r  , and the larger the r, the higher the accuracy. Chen [20] analyzed the integral and its analytical extension, where v is an arbitrary positive real number and r is large enough. Then the integral was transformed into a form that the integrand does not oscillate and decay exponentially quickly on [0, )  , which could be efficiently computed by using the Gauss-Laguerre quadrature rule. Zhang et al [21] added a new basis function in the interpolation process and proposed an improved algorithm of the Pravin method based on Hankel transform. The improved method is of high precision, high efficiency, and good adaptability to kernel function. However, the order range of this method is small, i.e., it is only applicable to the zero-order and first-order Hankel transformations. Niu [22] gave an efficient numerical algorithm for the integration with Bessel kernel based on the idea of the Filon-type method and the steepest descent method [23] . Cohen et al [24] proposed a linear transformation accelerated convergence algorithm (LTACA) to solve the alternating series. Wang et al [25] applied the LTACA to the numerical solution of the Hankel transform. The algorithm is simple to implement and can achieve high-accuracy, but it can only be used to solve non-negative integer order Hankel transform, and the efficiency of solving high-order Hankel transform is not high. Based on the variable upper bound integrals, the asymptotic quadrature rule of high frequency Hankel transform was derived [26] . However, the order of Bessel function is required to be non-negative, and the accuracy is low. Kisselev [27] proved the Hankel transform can be expressed by the absolutely and uniformly convergent series in reciprocal powers of parameter r, and can be used to obtain analytic expressions for the Hankel transform of an arbitrary integer order with the positive parameter r. But the conditions for the non-oscillating item are very stringent. The G-S inverse Laplace transform method (G-SILTM) uses a pure real number operation, and only a small number of Laplace transform variable values need to be calculated, so the efficiency of the algorithm is relatively high. However, the precision of the G-S algorithm is closely related to the finite word length of the computer. Therefore, this algorithm has higher requirements for computer hardware and software [28,29] .
The above numerical algorithms of the Hankel transform have disadvantages such as complex processing, low computational efficiency, higher requirements for the integrand function (excluding the Bessel factor), and a small range of application of the Hankel transform order. However, the G-SILTM has the advantages of high efficiency, simplicity, wide application range, and low requirements for the integrand function (excluding Bessel factor), so we develop an algorithm to calculate the arbitrary real-order Hankel transform.

Numerical Analysis of Hankel Transform
For  -order Hankel transform complex function, the real and imaginary parts can be reduced to the form of equation (1), respectively [30] .
, where  is a positive integer, and  makes  satisfy the large argument approximate condition of the Hankel function, i.e.
where the selection of m should be as small as possible on the premise of satisfying the large argument approximate condition.

Linear Transformation Accelerated Convergence Algorithm
For the alternating series , the linear transformation accelerates the convergence of the alternating series as follows [24,25] : Input: number of sum terms N; Initialize: (3 8) ; The algorithm is simple to implement and has high accuracy, and it is still effective to apply this method to the numerical calculation of the integer-order Hankel transform [25] . For equation (1) where m z is the value of the m th non-negative zero point of the  -order Bessel function ( ) J   of the first-kind normalized by the spatial distance  , and the first non-negative zero point is 0.
When ( ) F  is continuous and monotonic, due to the oscillation property of the Bessel function, equation (11) is the alternating series of the interval integral term, and equation (1) can be solved by applying the LTACA to equation (11). However, the application range of Hankel transformation order is relatively small, seeing the numerical examples for details.

G-S Inverse Laplace Transform Method
The expression of the G-SILTM is [33] L 1 where L F is the Laplace transform of f , ( , ) d j N are the filtering coefficients of the G-SILTM, [( 1) 2] j  is the integer part of ( 1) 2 j  , and min( , 2) j N is the minimum value of j and 2 N . The G-SILTM is mainly used to solve the transient sounding forward problem, and here we introduce it to solve arbitrary real-order Hankel transform.
The G-SILTM of sine and cosine transforms are Substituting equations (15) and (16) into equation The number of coefficient points of the G-SILTM is closely related to the word length of the computer. As far as the current general computer, we find that the 18-point coefficients are the best.
The influence of the value of the partition point  on the numerical results is shown in Ref. [31]. The following numerical examples are calculated with the interval partition point 10   , and 18 N  for the G-SILTM. At present, the LTACA has high accuracy in solving the Hankel transform, especially when the order of the Hankel transform is small. In order to illustrate the superiority of our algorithm, a detailed comparison is made in example 3.
Example 1 The calculation results are shown as Fig.1 and Table 1. Figure 1 shows the comparisons between the numerical solutions and the analytical solutions, and gives the relative errors at 0.01,0.1,1 b  . It can be seen that the relative error increases with the increase of b.

Fig. 6 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=3,4,5
It can be seen from Figs. 3-6 that the LTACA has a small relative error for the non-negative integer order Hankel transform. For 0,1, 2,3, 4,5 n  , the relative error is 14 (10 )   , and the relative error of the G-SILTM is larger and basically stable at 5 (10 )   . It can be seen from Table 2 that the time consumption of the LTACA increases significantly with the increase of the order of the Hankel function, while the time consumption of the G-SILTM basically remains unchanged with the increase of n. When the Hankel transform order reaches 4, the latter algorithm takes less time.  . It can be obtained that when the order of the Bessel function increases to a certain extent, the error decreases with the increase of variable a. That is, when n and a are larger, the relative error is smaller. As can be seen from Figs. 3-8, the LTACA has higher accuracy than the G-SILTM. It can be seen from Tables 2 and 3 that the time consumption of the LTACA increases significantly with the increase of order, but the time-consumption of the G-SILTM is more stable and will not increase significantly with the increase of n.  In addition, the LTACA needs to calculate the zero point of the Bessel function, and it is hard to obtain the zero of a fractional Bessel function at present, so it is not clear whether the algorithm is suitable for solving the fractional Hankel transform. Therefore, the LTACA can only be used to solve non-negative integer Hankel transform at present. However, the G-SILTM can be used to fractional Hankel transform, as described below.  As we can see from Figs. 9-13, the G-SILTM can solve not only the Hankel transform with order of rational number, but also the Hankel transform with order of irrational number. Therefore, the algorithm is suitable for the case of arbitrary real-order Hankel transform. The relative error decreases with the increase of the order.
It can be seen from Tables 1-4 that the G-SILTM can be used to solve the arbitrary real-order Hankel transform and the time consumption of the algorithm is relatively stable and short. Figure 14 (a) and (b) show the variations of relative error with the order of Hankel transform and the oscillation factor of the G-SILTM. It can be seen that the relative error is the largest when the order is the largest and the oscillation factor is the smallest. The relative errors at other positions are smaller. When the relative error is constant and the oscillation factor increases, the order of the Hankel transform will increase.
Based on the comparisons between the numerical solutions and the analytical solutions of the above examples, it can be seen that the G-SILTM can be employed to calculate the integral with arbitrary order Bessel function, and the algorithm is efficient. In addition, the algorithm can be used to calculate non-integer order Hankel transform.

Application in Vibration of Axisymmetric Infinite Membrane
Circular membranes are important parts of drums, pumps, microphones, and other devices. This accounts for their great importance in engineering. When the circular membrane is plane and its material is elastic and offers no resistance to bending (this excludes thin metal-lic membranes), the vibrations of the circular membrane is given in the form of two-dimensional wave equation in cylindrical coordinate system, ( , , ) r z  , i.e. [ The boundary conditions are Using zero order Hankel transform for equations (18) and (19), we can obtain a linear ordinary differential equation with boundary conditions General solution of equation (20) is Substituting equation (21) into (22) leads to The solution of equation (20) is Using inverse Hankel transform for equation (23), the solution of the problem is There is the integral result in Ref. [36] as follows: Equation (25) and the numerical solution of equation (24) are compared as Fig. 15.
It can be seen from Fig. 15 that the G-SILTM used to solve the vibration of axisymmetric infinite membrane is effective, and the relative error is 5 (10 )   .

Conclusion
In this paper, the numerical algorithm for the n th -order Hankel transform is studied and the G-SILTM algorithm is given. The numerical solutions of the examples are in good agreement with the analytical solutions, which fully demonstrates the correctness of the algorithm, but the algorithm requires convergence of the in- The comparisons of the numerical examples show that the LTACA has high accuracy, but it only works for non-negative integer Hankel transformations and the time consumption increases with n. The G-SILTM can overcome these drawbacks and compute arbitrary real-order Hankel transformations. The time consumption is relatively stable and short, but the accuracy is slightly lower. In addition, we use the G-SILTM to solve the vibration of axisymmetric infinite membrane.