A Numerical Algorithm for Arbitrary Real-Order Hankel Transform

Yonglin YANG; Xing LI; Shenghu DING; Wenshuai WANG

doi:10.1051/wujns/2022271026

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 Journal of Natural Sciences 2022, Vol.27 No.1, 26-34

Mathematics

CLC number: O242.2

A Numerical Algorithm for Arbitrary Real-Order Hankel Transform

Yonglin YANG, Xing LI, Shenghu DING and Wenshuai WANG^†

School of Mathematics and Statistics, Ningxia University, Yinchuan 750021, Ningxia, China

^† To whom correspondence should be addressed. E-mail: wws@nxu.edu.cn

Received: 5 November 2021

Abstract

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.

Key words: Hankel transform / large argument approximate expression of the Bessel function / linear transformation accelerated convergence algorithm (LTACA) / G-S inverse Laplace transform method (G-SILTM)

Biography: YANG Yonglin, male, Ph. D. candidate, research direction: singular integral equation and its application in new material mechanics. E-mail: yang_yongl@yeah.net

Foundation item: Supported by the National Natural Science Foundation of China (42064004, 12062022, 11762017, 11762016)

© Wuhan University 2022

This 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 $\int_{a}^{\infty} f (λ) J_{v} (λ r) d λ,$ but more stringent restrictions 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 $r ≫ 1$ , and the larger the r, the higher the accuracy. Chen^[20] analyzed the integral $\int_{a}^{\infty} f (λ) J_{v} (λ r) d λ,$ $(a \geq 0)$ 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, + \infty)$ , 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 $H_{ν} (α) = \int_{0}^{\infty} λ F (λ) J_{ν} (λ α) d λ$ (1)where $α \lor 0$ , $J_{ν} (λ α)$ is a Bessel function of the first kind of $ν (\geq - 1 / 2)$ -order. This article only discusses the case that $F (λ)$ is a real function. If $F (λ)$ is a complex function, the real and imaginary parts can be reduced to the form of equation (1), respectively^[30].

The integral interval $[0, \infty)$ is divided as $[0, ξ] \cup [ξ, \infty)$ , where $ξ$ is a positive integer, and $ξ$ makes $ξ α$ satisfy the large argument approximate condition of the Hankel function, i.e. $ξ α \geq m ≫ 1$ (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 $H_{ν} (α) = H_{ν}^{1} + H_{ν}^{2}$ (3)with $H_{ν}^{1} = \int_{0}^{ξ} λ F (λ) J_{ν} (λ α) d λ$ (4)and $H_{ν}^{2} = \int_{ξ}^{\infty} λ F (λ) J_{ν} (λ α) d λ$ (5)where $H_{ν}^{1}$ 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 $H_{ν}^{2}$ , the large argument approximate expression of Bessel function of real variable is ^[31,32] $J_{ν} (x) \approx \sqrt{\frac{2}{π x}} \cos [x - \frac{π}{2} (ν + \frac{1}{2})], x ≫ 1$ (6)Substituting equation (6) into (5) leads to $H_{ν}^{2} \approx \sqrt{\frac{2}{π α}} \int_{ξ}^{\infty} \sqrt{λ} F (λ) \cos [λ α - \frac{π}{2} (ν + \frac{1}{2})] d λ$ (7)where $\begin{array}{l} \cos [λ α - \frac{π}{2} (ν + \frac{1}{2})] = \cos (λ α) \cos [\frac{π}{2} (ν + \frac{1}{2})] \\ + \sin (λ α) \sin [\frac{π}{2} (ν + \frac{1}{2})] \end{array}$ (8)Substituting equation (8) into (7) leads to $\begin{array}{l} H_{ν}^{2} \approx \frac{1}{\sqrt{π α}} \cos [\frac{π}{2} (ν + \frac{1}{2})] \int_{ξ}^{\infty} \sqrt{λ} F (λ) \cos (λ α) d λ \\ + \sin [\frac{π}{2} (ν + \frac{1}{2})] \int_{ξ}^{\infty} \sqrt{λ} F (λ) \sin (λ α) d λ \end{array}$ (9)At this time, equation (9) is the following Fourier sine and cosine transform forms, $\begin{array}{l} f_{c} (k) = \int_{0}^{\infty} F (x) \cos (k x) d x \\ f_{s} (k) = \int_{0}^{\infty} F (x) \sin (k x) d x \end{array}$ (10)

2 Linear Transformation Accelerated Convergence Algorithm

For the alternating series $S = \sum_{k = 0}^{\infty} {(- 1)}^{k} a_{k}$ , the linear transformation accelerates the convergence of the alternating series as follows ^[24,25]:

Input: number of sum terms N;

Initialize: $\begin{array}{l} d = (3 + \sqrt{8})^{N}; \\ d = (d + 1 / d) / 2; \\ b = 1; c = - d; s = 0; \end{array}$

Loop: $\begin{array}{l} for k = 0 : N - 1 \\ c = b - c; \\ s = s + c \cdot a_{k}; \\ b = (k + N) (k - N) / ((k + 1 / 2) (k + 1)); \\ s / d; \\ end \end{array}$

Output: $S = s / d$ .

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 $H_{n} (α) \approx S = \sum_{m = 0}^{N} p_{m}$ (11)with $p_{m} = \int_{z_{m}}^{z_{m + 1}} λ F (λ) J_{ν} (λ α) d λ$ (12)where $z_{m}$ 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.

3 G-S Inverse Laplace Transform Method

The expression of the G-SILTM is ^[33] $f (t) \approx \frac{\ln 2}{t} \sum_{j = 1}^{N} d (j, N) F_{L} (j \frac{\ln 2}{t})$ (13)with $\begin{array}{l} d (j, N) = (- 1)^{j + N / 2} \\ \cdot \sum_{k = [(j + 1) / 2]}^{\min (j, N / 2)} \frac{k^{N / 2} (2 k)!}{(N / 2 - k)! k! (k - 1)! (j - k)! (2 k - j)!} \end{array}$ (14)where $F_{L}$ is the Laplace transform of $f$ , $d (j, N)$ are the filtering coefficients of the G-SILTM, $[(j + 1) / 2]$ is the integer part of $(j + 1) / 2$ , and $\min (j, N / 2)$ is the minimum value of j and $N / 2$ . 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 $\begin{array}{l} f_{c} (t) = \frac{\ln 2}{t} \sum_{j = 1}^{N} d (j, N) F_{Lc} (j \frac{\ln 2}{t}) \\ = {(\frac{\ln 2}{t})}^{2} \sum_{j = 1}^{N} j d (j, N) \int_{0}^{\infty} \frac{F (s) d s}{s^{2} + {(j \frac{\ln 2}{t})}^{2}} \end{array}$ (15) $\begin{array}{l} f_{s} (t) = \frac{\ln 2}{t} \sum_{j = 1}^{N} d (j, N) F_{Ls} (j \frac{\ln 2}{t}) \\ = \frac{\ln 2}{t} \sum_{j = 1}^{N} d (j, N) \int_{0}^{\infty} \frac{s F (s) d s}{s^{2} + {(j \frac{\ln 2}{t})}^{2}} \end{array}$ (16)

Let $F (s)$ make the integral $\int_{0}^{\infty} \frac{s F (s) d s}{s^{2} + C},$ $\int_{0}^{\infty} \frac{F (s) d s}{s^{2} + C}$ convergent, where C is an arbitrary constant.

Substituting equations (15) and (16) into equation (9), the G-SILTM of $H_{ν}^{2}$ is $\begin{array}{l} H_{ν}^{2} \approx \frac{1}{\sqrt{π α}} \cos [\frac{π}{2} (ν + \frac{1}{2})] {{(\frac{\ln 2}{α})}^{2} \sum_{j = 1}^{N} j d (j, N) \int_{0}^{\infty} \frac{\sqrt{s} F (s) d s}{s^{2} + {(j \frac{\ln 2}{α})}^{2}} - \int_{0}^{ξ} \sqrt{λ} F (λ) \cos (λ α) d λ} \\ + \frac{1}{\sqrt{π α}} \sin [\frac{π}{2} (ν + \frac{1}{2})] {\frac{\ln 2}{α} \sum_{j = 1}^{N} d (j, N) \int_{0}^{\infty} \frac{s \sqrt{s} F (s) d s}{s^{2} + {(j \frac{\ln 2}{α})}^{2}} - \int_{0}^{ξ} \sqrt{λ} F (λ) \sin (λ α) d λ} \end{array}$ (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 $ξ = 10$ , and $N = 18$ 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 $f (a) = \int_{0}^{\infty} λ \exp (- a λ^{2}) J_{0} (λ b) d λ = \frac{1}{ 2 a} \exp (- \frac{b^{2}}{4 a})$ The calculation results are shown as Fig.1 and Table 1.

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 $b = 0.01, 0.1, 1$ . It can be seen that the relative error increases with the increase of b.

Example 2 $f (b) = \int_{0}^{\infty} λ \exp (- a λ) J_{1} (λ b) d λ = \frac{b}{ {(a^{2} + b^{2})}^{3 / 2}}$

Let $a = 2$ , the calculation results are as Fig.2.

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 $10^{- 4}$ .

Example 3 $\begin{array}{l} f (a) = \int_{0}^{\infty} \exp (- μ λ) J_{n} (λ a) d λ \\ = \frac{1}{\sqrt{μ^{2} + a^{2}}} {(\frac{a}{\sqrt{μ^{2} + a^{2}} + μ})}^{n} \end{array}$ where n is an arbitrary real.

Let $μ = 2$ , $n =0,1,2,3,4,5$ , the comparison results of the G-SILTM and the LTACA are as shown in Figs. 3-6.

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

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

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

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 $n = 0, 1, 2, 3, 4, 5$ , the relative error is $Ο (10^{- 14})$ , and the relative error of the G-SILTM is larger and basically stable at $Ο (10^{- 5})$ .

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 $n = 20, 30$ . 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.

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

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 $μ = 2$ , $n = - \frac{1}{2}, \frac{1}{2}, - \frac{1}{4}, \frac{4}{5}, \frac{16}{3}, \frac{15}{2}, - \frac{\sqrt{2}}{4}, - \frac{\sqrt{7}}{7}, \frac{\sqrt{2}}{3},$ $\sqrt{46}$ , the calculation results of the G-SILTM are shown as Figs.9-13.

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

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

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

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

n = \frac{- \sqrt{2}}{4}, \frac{- \sqrt{7}}{7}

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

n = \frac{\sqrt{2}}{3}, \sqrt{46}

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.

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 $ξ = 10$ and N=18

Table 2

The time consumption of the two algorithms with $ξ = 10$ , N=18 and $μ = 2$ 　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, $(r, θ, z)$ , i.e.^[35] $\begin{array}{l} \frac{1}{r^{2}} \frac{\partial^{2} u}{\partial θ^{2}} + \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial u}{\partial r}) = \frac{1}{a^{2}} \frac{\partial^{2} u}{\partial t^{2}}, \\ r \in (0, \infty), t \lor 0, θ \in [0, 2 π] \end{array}$

Let the membrane initially be the bell shape, the axisymmetric infinite vibration problem can be simplified as $\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial u}{\partial r}) = \frac{1}{a^{2}} \frac{\partial^{2} u}{\partial t^{2}}, r \in (0, \infty), t \lor 0$ (18)The boundary conditions are $u (r, 0) = u_{0} (r) = \frac{1}{\sqrt{r^{2} + 1}}, \frac{\partial u (r, 0)}{\partial t} = u_{1} (r) = 0$ (19)Using zero order Hankel transform for equations (18) and (19), we can obtain a linear ordinary differential equation $\frac{\partial^{2} \bar{u}}{\partial t^{2}} + a^{2} λ^{2} \bar{u} = 0$ (20)with boundary conditions ${\bar{u}}_{0} (r) = \frac{1}{λ} \exp (- λ), {\bar{u}}_{1} (r) = 0$ (21)where $\bar{u} (λ, t) = \int_{0}^{\infty} r J_{0} (λ r) u (r, t) d r$ .

General solution of equation (20) is ${\begin{cases} \bar{u} (λ, t) = c_{1} \cos (a λ t) + c_{2} \sin (a λ t) \\ \frac{\partial}{\partial t} \bar{u} (λ, t) = - a λ c_{1} \sin (a λ t) + a λ c_{2} \cos (a λ t) \end{cases}$ (22)Substituting equation (21) into (22) leads to $c_{1} = {\bar{u}}_{0} (r), c_{2} = \frac{{\bar{u}}_{1} (r)}{a λ}$ The solution of equation (20) is $\bar{u} (λ, t) = \frac{1}{λ} \exp (- λ) \cos (a λ t)$ (23)Using inverse Hankel transform for equation (23), the solution of the problem is $u (r, t) = \int_{0}^{\infty} J_{0} (λ r) \exp (- λ) \cos (a λ t) d λ$ (24)There is the integral result in Ref. [36] as follows: $\begin{array}{l} \int_{0}^{\infty} J_{0} (c λ) \exp (- b λ) \cos (a λ) d λ \\ = \frac{{[\sqrt{{(b^{2} + c^{2} - a^{2})}^{2} + 4 a^{2} b^{2}} + b^{2} + c^{2} - a^{2}]}^{1 / 2}}{\sqrt{2} \sqrt{{(b^{2} + c^{2} - a^{2})}^{2} + 4 a^{2} b^{2}}} \end{array}$ Let $a = a t, b = 1, r = c$ , the analytical solution of the problem is $u (r, t) = \frac{{[\sqrt{{(1 + r^{2} - a^{2} t^{2})}^{2} + 4 a^{2} t^{2}} + 1 + r^{2} - a^{2} t^{2}]}^{1 / 2}}{\sqrt{2} \sqrt{{(1 + r^{2} - a^{2} t^{2})}^{2} + 4 a^{2} t^{2}}}$ (25)

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

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

a = 1, t = 1, r \in [0.01, 150]

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 $Ο (10^{- 5})$ .

6 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 integrals $\int_{0}^{\infty} \frac{s F (s) d s}{s^{2} + C}$ and $\int_{0}^{\infty} \frac{F (s) d s}{s^{2} + C}$ .

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

Kaufman A A, Keller G V. Frequency and Transient Soundings [M]. Beijing: Geological Publishing House, 1987. [Google Scholar]
Fang W Z, Li Y G, Li X. Theory of TEM Sounding [M]. Xian: Press of North-west Industry University, 1993(Ch). [Google Scholar]
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]
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]
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]
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]
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]
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]
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]
Wang W J. Propagation of the Bessel Gaussian Beam in Turbulence and Application [D]. Xi’an: Xidian University, 2019(Ch). [Google Scholar]
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]
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]
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]
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]
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]
Guptasarma D, Singh B. New digital linear filters for Hankel J₀ and J₁ transforms [J]. Geophysical Prospecting, 1997, 45(5): 745-762. [CrossRef] [Google Scholar]
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]
Yu R. Two Types of Bessel Transform Its Numerical Integration Methods [D]. Changsha: Central South University, 2007(Ch). [Google Scholar]
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]
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]
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]
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]
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]
Cohen H, Villegas F R, Zagier D. Convergence acceleration of alternating series [J]. Experimental Mathematics, 2000, 9(1): 3-12. [CrossRef] [MathSciNet] [Google Scholar]
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]
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]
Kisselev A V. Exact expansions of Hankel transforms and related integrals [J]. The Ramanujan Journal, 2020, 55: 349-367. [Google Scholar]
Raiche A P. Transient electromagnetic field computations for polygonal loops on layered earth [J]. Geophysics, 1987, 526: 785-793. [NASA ADS] [CrossRef] [Google Scholar]
Villinger H. Solving cylindrical geothermal problems using Gaver-Stehfest inverse Laplace transform [J]. Geophysics, 1985, 50(10): 1581-1587. [NASA ADS] [CrossRef] [Google Scholar]
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]
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]
Zeidler E. Teubner-Taschenbuch der Mathematik [M]. Berlin: Vieweg+Teubner Verlag, 2012. [Google Scholar]
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]
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]
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]
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 $ξ = 10$ and N=18

In the text

Table 2

The time consumption of the two algorithms with $ξ = 10$ , N=18 and $μ = 2$ 　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

	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

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

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

	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

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

	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

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

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

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

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

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

	Fig. 12 Comparisons and relative errors between the analytical and the numerical solutions at $n = \frac{- \sqrt{2}}{4}, \frac{- \sqrt{7}}{7}$
In the text

	Fig. 13 Comparison and relative error between the analytical and the numerical solutions at $n = \frac{\sqrt{2}}{3}, \sqrt{46}$
In the text

	Fig. 14 Variations of relative error with n and a
In the text

	Fig. 15 Comparisons and relative errors between the analytical and the numerical solutions at $a = 1, t = 1, r \in [0.01, 150]$
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.

[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 J₀ and J₁ 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]