EDP Sciences logo
Submit your paper
Search
Menu
All issues Volume 27 / No 1 (March 2022) Wuhan Univ. J. Nat. Sci., 27 1 (2022) 26-34 Full HTML
Open Access
Issue
Wuhan Univ. J. Nat. Sci.
Volume 27, Number 1, March 2022
Page(s) 26 - 34
DOI https://doi.org/10.1051/wujns/2022271026
Published online 16 March 2022

© Wuhan University 2022

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

0 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-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-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-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 the integral but more stringent restrictions for 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 , 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 , 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.

1 Numerical Analysis of Hankel Transform

For -order Hankel transform(1)where , is a Bessel function of the first kind of -order. This article only discusses the case that is a real function. If is a complex function, the real and imaginary parts can be reduced to the form of equation (1), respectively[30].

The integral interval is divided as , where is a positive integer, and makes satisfy the large argument approximate condition of the Hankel function, i.e.(2)where the selection of m should be as small as possible on the premise of satisfying the large argument approximate condition.

Using the accumulation of integral, equation (1) is rewritten as(3)with(4)and(5)where is a normal integral, and it can be solved by any general numerical algorithm. The Simpson algorithm is used in the numerical examples of Section 4. For , the large argument approximate expression of Bessel function of real variable is [31,32](6)Substituting equation (6) into (5) leads to(7)where(8)Substituting equation (8) into (7) leads to(9)At this time, equation (9) is the following Fourier sine and cosine transform forms,(10)

2 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:

Loop:

Output: .

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), let(11)with(12)where is the value of the mth non-negative zero point of the -order Bessel function of the first-kind normalized by the spatial distance , and the first non-negative zero point is 0.

When 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.

3 G-S Inverse Laplace Transform Method

The expression of the G-SILTM is [33](13)with(14)where is the Laplace transform of , are the filtering coefficients of the G-SILTM, is the integer part of , and is the minimum value of j and . 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(15)(16)

Let make the integral convergent, where C is an arbitrary constant.

Substituting equations (15) and (16) into equation (9), the G-SILTM of is(17)

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.

4 Numerical Examples

The numerical examples are from Refs. [16,21,34]. 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 , and 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.

thumbnail Fig. 1 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at b=0.01, 0.1, 1

Figure 1 shows the comparisons between the numerical solutions and the analytical solutions, and gives the relative errors at . It can be seen that the relative error increases with the increase of b.

Example 2  

Let , the calculation results are as Fig.2.

thumbnail Fig. 2 Comparison and relative error between the analytical and the numerical solutions of the G-SILTM

Figure 2 illustrates the comparison between the numerical solution and the analytical solution of the G-SILTM, as well as its relative error. The numerical solution of the example is consistent with the analytical solution, and the relative error is in the range of .

Example 3  where n is an arbitrary real.

Let , , the comparison results of the G-SILTM and the LTACA are as shown in Figs. 3-6.

thumbnail Fig. 3 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=0,1,2

thumbnail Fig. 4 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=0,1,2

thumbnail Fig. 5 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=3,4,5

thumbnail 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 , the relative error is , and the relative error of the G-SILTM is larger and basically stable at .

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.

Figures 7 and 8 show the comparisons and relative errors between the numerical solutions and the analytical solutions of the G-SILTM and the LTACA of the example at . 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.

thumbnail Fig. 7 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=20, 30

thumbnail Fig. 8 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=20, 30

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.

Let , , the calculation results of the G-SILTM are shown as Figs.9-13.

thumbnail Fig. 9 Comparisons and relative errors between the analytical and the numerical solutions at n=-1/2, 1/2

thumbnail Fig. 10 Comparisons and relative errors between the analytical and the numerical solutions at n=-1/4, 4/5

thumbnail Fig. 11 Comparisons and relative errors between the analytical and the numerical solutions at n=16/3, 15/2

thumbnail Fig. 12 Comparisons and relative errors between the analytical and the numerical solutions at

thumbnail Fig. 13 Comparison and relative error between the analytical and the numerical solutions at

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.

thumbnail Fig. 14 Variations of relative error with n and a

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.

Table 1

The time consumption of the main body of the algorithm with and N=18

Table 2

The time consumption of the two algorithms with , N=18 and 　s

Table 3

The time consumption of the two algorithms　s

Table 4

The time consumption of the main body of the G-SILTM

5 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 metallic membranes), the vibrations of the circular membrane is given in the form of two-dimensional wave equation in cylindrical coordinate system, , i.e.[35]

Let the membrane initially be the bell shape, the axisymmetric infinite vibration problem can be simplified as(18)The boundary conditions are(19)Using zero order Hankel transform for equations (18) and (19), we can obtain a linear ordinary differential equation(20)with boundary conditions(21)where .

General solution of equation (20) is(22)Substituting equation (21) into (22) leads toThe solution of equation (20) is(23)Using inverse Hankel transform for equation (23), the solution of the problem is(24)There is the integral result in Ref. [36] as follows:Let , the analytical solution of the problem is(25)

Equation (25) and the numerical solution of equation (24) are compared as Fig. 15.

thumbnail Fig. 15 Comparisons and relative errors between the analytical and the numerical solutions at

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 .

6 Conclusion

In this paper, the numerical algorithm for the nth-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 integrals and .

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.

References

  1. Kaufman A A, Keller G V. Frequency and Transient Soundings [M]. Beijing: Geological Publishing House, 1987. [Google Scholar]
  2. Fang W Z, Li Y G, Li X. Theory of TEM Sounding [M]. Xian: Press of North-west Industry University, 1993(Ch). [Google Scholar]
  3. Wen A H, Wang X Q. Utilizing direct integration to enhance calculation accuracy of 1D electromagnetic response for current dipole source [J]. Northwestern Seismological Journal, 2003, 25(3): 193-197(Ch). [Google Scholar]
  4. Ji G, Zhang W K, Lu X P. Numerical integration on Green source function with free water surface [J]. Journal of Naval University of Engineering, 2004, 16(6): 89-92(Ch). [Google Scholar]
  5. He J S, Bao L Z. EM field of vertical wire electric source and its practical meaning [J]. Journal of Central South University (Science and Technology), 2011, 42(1): 130-135(Ch). [Google Scholar]
  6. Singh B M, Rokne J, Dhaliwal R S. Diffraction of antiplane shear waves by a finite crack in a piezoelectric material [J]. ZAMM-Journal of Applied Mathematics and Mechanics /Zeitschrift Für Angewandte Mathematik Und Mechanik, 2011, 91(11): 866-874. [CrossRef] [MathSciNet] [Google Scholar]
  7. Yang Y L, Li X, Wang W S. Wettability of semispherical droplets on layered elastic gradient soft substrates [J]. Scientific Reports, 2021, 11(1): 2236. [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  8. Zhang Q A, Wu F T, Zheng W T, et al. Self-reconstructing properties of high-order Besssel-Gauss beam [J]. Scientia Sinica Physica, Mechanica & Astronomica, 2011, 41(10): 1131-1137. [NASA ADS] [CrossRef] [Google Scholar]
  9. You K M, Lin Y L, Wang Y W, et al. High order Hankel transform based on Dini expansion and its applications in beam propagation [J]. Acta Physica Sinica, 2013, 62(14): 30-35. [Google Scholar]
  10. Wang W J. Propagation of the Bessel Gaussian Beam in Turbulence and Application [D]. Xi’an: Xidian University, 2019(Ch). [Google Scholar]
  11. Meng X S. Electromagnetic Vortex Wave Generation and Target Near-Field Scattering Based on Artificial Electromagnetic Metasurface [D]. Xi’an: Xidian University, 2019(Ch). [Google Scholar]
  12. Blakemore M, Evans G, Hyslop J. Comparison of some methods for evaluating infinite range oscillatory integrals [J]. Journal of Computational Physics, 1976, 22(3): 352-376. [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  13. Chave A D. Numerical integration of related Hankel transforms by quadrature and continued fraction expansion [J]. Geophysics, 1983, 48(12): 1671-1686. [NASA ADS] [CrossRef] [Google Scholar]
  14. Cree M J, Bones P J. Algorithms to numerically evaluate the Hankel transform [J]. Computers & Mathematics with Applications, 1993, 26(1): 1-12. [CrossRef] [MathSciNet] [Google Scholar]
  15. Lucas S K, Stone H A. Evaluating infinite integrals involving Bessel functions of arbitrary order [J]. Journal of Computational and Applied Mathematics, 1995, 64(3): 217-231. [CrossRef] [MathSciNet] [Google Scholar]
  16. Guptasarma D, Singh B. New digital linear filters for Hankel J0 and J1 transforms [J]. Geophysical Prospecting, 1997, 45(5): 745-762. [CrossRef] [Google Scholar]
  17. Filon L N G. On a quadrature formula trigonometric integrals [J]. Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 1928, 49: 38-47. [Google Scholar]
  18. Yu R. Two Types of Bessel Transform Its Numerical Integration Methods [D]. Changsha: Central South University, 2007(Ch). [Google Scholar]
  19. Chen R Y. Numerical analysis on a class of integrals involving Bessel function [J]. Journal of Chongqing Institute of Technology (Natural Science), 2008, 22(11): 83-88(Ch). [Google Scholar]
  20. Chen R Y. Numerical approximations to integrals with a highly oscillatory Bessel kernel [J]. Applied Numerical Mathematics, 2012, 62(5): 636-648(Ch). [CrossRef] [MathSciNet] [Google Scholar]
  21. Zhang H Q, Chen Y, Nie X. Fast Hankel transforms algorithm based on kernel function interpolation with exponential functions [J]. Journal of Applied Mathematics, 2014, 2014: 1-7. [Google Scholar]
  22. Niu F F. A Numerical Method for Computing Highly Oscillatory Integrals with a Bessel Kernel [D]. Wuhan: Huazhong University of Science &Technology, 2015(Ch). [Google Scholar]
  23. Majidian H. Numerical approximation of highly oscillatory integrals on semi-finite intervals by steepest descent method [J]. Numerical Algorithms, 2012, 63(3): 537-548. [Google Scholar]
  24. Cohen H, Villegas F R, Zagier D. Convergence acceleration of alternating series [J]. Experimental Mathematics, 2000, 9(1): 3-12. [CrossRef] [MathSciNet] [Google Scholar]
  25. Wang H J, Li G. Hankel transform to accelerate the convergence of the numerical integration algorithm [J]. Technological Development of Enterprise, 2012, 31(16): 8-9+44(Ch). [Google Scholar]
  26. Gao J Z, Chen R Y. On computation of high frequency Hankel transforms [J]. Alexandria Engineering Journal, 2019, 58(3): 1033-1037. [CrossRef] [Google Scholar]
  27. Kisselev A V. Exact expansions of Hankel transforms and related integrals [J]. The Ramanujan Journal, 2020, 55: 349-367. [Google Scholar]
  28. Raiche A P. Transient electromagnetic field computations for polygonal loops on layered earth [J]. Geophysics, 1987, 526: 785-793. [NASA ADS] [CrossRef] [Google Scholar]
  29. Villinger H. Solving cylindrical geothermal problems using Gaver-Stehfest inverse Laplace transform [J]. Geophysics, 1985, 50(10): 1581-1587. [NASA ADS] [CrossRef] [Google Scholar]
  30. Li D Z, Yang Z Q, Liang Z. Numerical computation of Bessel function for complex arguments [J]. Journal of University of Electronic Science and Technology of China, 1996, 25(7): 125-128(Ch). [Google Scholar]
  31. Hua J, Jiang Y S, Wang W B. The numerical integration of dual Hankel transformation [J]. Coalgeology and Exploration, 2001, 29(3): 58-62. [Google Scholar]
  32. Zeidler E. Teubner-Taschenbuch der Mathematik [M]. Berlin: Vieweg+Teubner Verlag, 2012. [Google Scholar]
  33. Luo H G. The Study about 1D Forward Modeling of Large-Fixed Loop TEM [D]. Beijing:China University of Geosciences (Beijing), 2012(Ch) . [Google Scholar]
  34. Anderson W L. Fast Hankel transforms using related and lagged convolutions [J]. ACM Transactions on Mathematical Software, 1982, 8(4): 344-368. [CrossRef] [Google Scholar]
  35. Kasemsuwan J, Sabau S V, Somboon U. Differential transformation method for circular membrane vibrations [J]. Bulletin of the Transilvania University of Brasov, Series III: Mathematics, Informatics, Physics, 2020, 61(12): 333-350 [CrossRef] [Google Scholar]
  36. Gradshteyn I S, Ryzhik I M. Table of Integrals, Series, and Products [M]. Seventh Edition. Amsterdam: Elsevier, 2007: 743. [MathSciNet] [Google Scholar]

All Tables

Table 1

The time consumption of the main body of the algorithm with and N=18

In the text
Table 2

The time consumption of the two algorithms with , N=18 and 　s

In the text
Table 3

The time consumption of the two algorithms　s

In the text
Table 4

The time consumption of the main body of the G-SILTM

In the text

All Figures

thumbnail Fig. 1 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at b=0.01, 0.1, 1
In the text
thumbnail Fig. 2 Comparison and relative error between the analytical and the numerical solutions of the G-SILTM
In the text
thumbnail Fig. 3 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=0,1,2
In the text
thumbnail Fig. 4 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=0,1,2
In the text
thumbnail Fig. 5 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=3,4,5
In the text
thumbnail Fig. 6 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=3,4,5
In the text
thumbnail Fig. 7 Comparisons and relative errors between the analytical and the numerical solutions of the LTACA at n=20, 30
In the text
thumbnail Fig. 8 Comparisons and relative errors between the analytical and the numerical solutions of the G-SILTM at n=20, 30
In the text
thumbnail Fig. 9 Comparisons and relative errors between the analytical and the numerical solutions at n=-1/2, 1/2
In the text
thumbnail Fig. 10 Comparisons and relative errors between the analytical and the numerical solutions at n=-1/4, 4/5
In the text
thumbnail Fig. 11 Comparisons and relative errors between the analytical and the numerical solutions at n=16/3, 15/2
In the text
thumbnail Fig. 12 Comparisons and relative errors between the analytical and the numerical solutions at
In the text
thumbnail Fig. 13 Comparison and relative error between the analytical and the numerical solutions at
In the text
thumbnail Fig. 14 Variations of relative error with n and a
In the text
thumbnail Fig. 15 Comparisons and relative errors between the analytical and the numerical solutions at
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.