# Fast Hankel Transforms Algorithm Based on Kernel Function Interpolation with Exponential Functions

Authors: Huaiqing Zhang; Yu Chen; Xin Nie
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The Pravin method for Hankel transforms is based on the decomposition of kernel function with exponential function. The defect of such method is the difficulty in its parameters determination and lack of adaptability to kernel function especially using monotonically decreasing functions to approximate the convex ones. This thesis proposed an improved scheme by adding new base function in interpolation procedure. The improved method maintains the merit of Pravin method which can convert the Hankel integral to algebraic calculation. The simulation results reveal that the improved method has high precision, high efficiency, and good adaptability to kernel function. It can be applied to zero-order and first-order Hankel transforms.


---

## Body

## 1. Introduction

The Hankel transforms (HT) arise naturally in the discussion of problems posed in cylindrical coordinates and hence [1], as a result of separation of variables, involving Bessel functions. The HT are frequently used as a tool for solving numerous scientific problems. For example, the stratified model and cylindrical coordinates are widely used in geophysical research, and the HT arise in forward and inverse calculation with zero or first order. The general Hankel transforms pair with the kernel being J
ν is defined as
(1)
F
ν
(
r
)
=
∫
0
∞
λ
f
(
λ
)
J
ν
(
λ
r
)
d
λ
and HT, being self-reciprocal, inverse is given by
(2)
f
(
λ
)
=
∫
0
∞
r
F
ν
(
r
)
J
ν
(
λ
r
)
d
r
,
where J
ν is the νth-order Bessel function of first kind. Analytical evaluations are rare and their numerical computations are difficult because of the oscillatory behavior of the Bessel function and the infinite length of the interval.To overcome these difficulties, various different techniques are available in the literature. The first is the fast Hankel transforms as proposed by Siegman in [2]. Here, via an exponential change of variables, the problem is transformed in the space of the logarithmic coordinates and the fast Fourier transform in that space. The disadvantage is requiring sampling over an exponential grid, thereby leading to important errors for functions with an oscillating tail. Moreover, it is sensitive to the smoothness of the functions in that space [3]. The second is the back-projection and projection-slice methods [4], which carry out the HT as a double integral by means of one of the standard integral representations of the Bessel functions. But these methods generally require the efficient implementation of Chebyshev and Abel transforms. The computational complexity is unfortunately O
(
N
2
) [5]. In [6], the authors used Filon quadrature philosophy to evaluate zero-order Hankel transforms. They separated the integrand into the product of slowly varying component and a rapidly oscillating one. This method works quite well for computing F
0
(
r
) for r
⩾
1, but the error is appreciable for r
<
1. And the calculation of inverse HT is more difficult as F
0
(
r
) is no longer a smooth function but a rapidly oscillating one. In 1998, Yu et al. [7] gave another method to compute zero-order quasi-discrete HT by approximating the input function by a Fourier-Bessel series over a finite integration interval. It leads to a symmetric transformation matrix. And later in 2004, Guizar-Sicairos and Gutiérrez-Vega [8] obtained a powerful scheme to calculate the HT of order n by extending the zero-order HT algorithm of Yu to higher orders. Their algorithm is based on the orthogonality properties of Bessel functions.More recently, Postnikov [9] and Zykov and Postnikov [10] proposed, for the first time, a novel and powerful method for computing zero- and first-order HT by using Haar basis and piecewise-linear basis. After that, Singh et al. adopted the linear Legendre multiwavelets [3], rationalized Haar (RH) wavelets [11], and the hybrid method of block-pulse and Legendre polynomials [1] for small random perturbation in data function situation [12]. Above methods can be cast into a general class as expansion of integrand by wavelets.Similarly, Gupta et al. [13] have proposed to transform the integral kernel function using orthogonal exponential expansion. This expansion reduces the integral to a simple algebraic sum because the analytical solution of HT of an exponential function is readily available. The essence of Pravin method is adopting exponential functions for fitting or interpolating basis [14]. However, only the monotonic exponential base function is selected, and, as a result, it is difficult to approximate the convex functions. So in this paper, we extend this idea further and obtain an improvement method by adding a new base function for fast calculation of HT. Numerical examples and comparison were given to illustrate the proposed algorithms.


## 2. Calculation Principles Based on Kernel Function Decomposition

### 2.1. Kernel Function Interpolation Based on the Exponential Function

In Pravin’s innovative idea of exponential expansion, the kernel functionG
(
λ
) (namely, λ
f
(
λ
)) which is defined as the part to be multiplied with Bessel function is decomposed by a series of exponential functions. In other words, the kernel function represents a linear combination of exponential functions and the coefficient of every exponential function can be determined by interpolation. Because the analytical solution of HT of an exponential function is readily available in literature, so we can finally obtain the integral through a simple algebraic sum. Since the order of Bessel function in the majority of physical problems is an integer as zero or one, this thesis focuses on zero- and first-order Hankel transforms. The calculation principles based on the exponential function expansion can be summarized as the following.(1)
Decomposition of kernel function using a linear combination form is as the following:(3)
G
(
λ
)
=
λ
f
(
λ
)
=
∑
i
=
1
N
a
i
e
-
ε
i
λ
,
whereε
i is the approximation parameter of the exponential function and it is determined by orthogonalization in Pravin method and a
i is the coefficient and it can be determined by collocation method. The interpolation form of above formula can be expressed briefly in matrix as G
=
M
a, so coefficients a
=
M
-
1
G, where matrix element is M
j
i
=
e
-
ε
i
λ
j.
(2)
Calculation of every component’s integral is as follows: the integral for the product of every base function with Bessel function is given in explicit form as(4)
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
=
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
v
.
(3)
HT calculation procedure is as follows: based on collocation and explicit formula, we can obtain(5)
F
ν
(
r
)
=
∫
0
∞
λ
f
(
λ
)
J
ν
(
λ
r
)
d
λ
=
∑
i
=
1
N
a
i
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
v
.And it can be expressed as F
=
N
a, where the matrix element is N
j
i
=
(
1
/
ε
i
2
+
r
j
2
)
[
r
j
/
ε
i
2
+
r
j
2
+
ε
i
]
v.It is obvious that the Pravin method is simple and concise. And the calculation error is only caused by the interpolation. Nevertheless, the base functions are monotonically decreasing exponential functions. So it is difficult to approximate convex function, such as approximating the kernel function whose value is zero at the origin. Moreover, to determine the parameterε is also hard in order to ensure the method’s accuracy and stability. Thus, this paper proposed to add a new type base function, which has a convex form and its value is zero at the origin.


### 2.2. Improved Kernel Function Interpolation Theory

Based on the idea of exponential expansion, the new type base function should satisfy the following two conditions. Specifically, the function has convexity and explicit form of HT. After systematically analyzing, we choose the following function as(6)
ϕ
(
λ
)
=
λ
ν
+
1
e
-
ε
λ
2
,
where ν is zero or one, correspondingly to the order of Bessel function. And the improved interpolation formula is as the following:
(7)
G
(
λ
)
=
λ
f
(
λ
)
=
∑
i
=
1
N
a
i
e
-
ε
i
λ
+
∑
i
=
1
N
b
i
λ
ν
+
1
e
-
ε
i
λ
2
.
It can be written as G
=
Ф
w
=
[
Ф
1
,
Ф
2
]
[
a
,
b
]
T, and coefficients w
=
[
a
,
b
]
T
=
[
Ф
1
,
Ф
2
]
-
1
G; the elements are Φ
j
i
1
=
e
-
ε
i
λ
j and Φ
j
i
2
=
λ
j
v
+
1
e
-
ε
i
λ
j
2. So the coefficients a
,
b can be determined by collocation. The HT of (7) can be obtained as
(8)
F
ν
(
r
)
=
∑
i
=
1
N
a
i
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
+
∑
i
=
1
N
b
i
∫
0
∞
λ
ν
+
1
e
-
ε
i
λ
2
J
ν
(
λ
r
)
d
λ
.
The integrals of the second component can be given by the following formula:
(9)
∫
0
∞
λ
ν
+
1
e
-
ε
λ
2
J
ν
(
λ
r
)
d
λ
=
r
ν
(
2
ε
)
ν
+
1
e
-
(
r
2
/
4
ε
)
.
For v
=
0, referring to formulas (4) and (9), we can obtain zero-order HT as
(10)
F
0
(
r
)
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
+
∑
i
=
1
N
b
i
1
2
ε
i
e
-
(
r
2
/
4
ε
i
)
.
And for v
=
1, we can obtain the first-order HT as
(11)
F
1
(
r
)
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
+
∑
i
=
1
N
b
i
r
(
2
ε
i
)
2
e
-
(
r
2
/
4
ε
i
)
.
So the improvement method maintains the simple and concise characteristics. The HT calculation can be obtained by an algebraic summation.


## 2.1. Kernel Function Interpolation Based on the Exponential Function

In Pravin’s innovative idea of exponential expansion, the kernel functionG
(
λ
) (namely, λ
f
(
λ
)) which is defined as the part to be multiplied with Bessel function is decomposed by a series of exponential functions. In other words, the kernel function represents a linear combination of exponential functions and the coefficient of every exponential function can be determined by interpolation. Because the analytical solution of HT of an exponential function is readily available in literature, so we can finally obtain the integral through a simple algebraic sum. Since the order of Bessel function in the majority of physical problems is an integer as zero or one, this thesis focuses on zero- and first-order Hankel transforms. The calculation principles based on the exponential function expansion can be summarized as the following.(1)
Decomposition of kernel function using a linear combination form is as the following:(3)
G
(
λ
)
=
λ
f
(
λ
)
=
∑
i
=
1
N
a
i
e
-
ε
i
λ
,
whereε
i is the approximation parameter of the exponential function and it is determined by orthogonalization in Pravin method and a
i is the coefficient and it can be determined by collocation method. The interpolation form of above formula can be expressed briefly in matrix as G
=
M
a, so coefficients a
=
M
-
1
G, where matrix element is M
j
i
=
e
-
ε
i
λ
j.
(2)
Calculation of every component’s integral is as follows: the integral for the product of every base function with Bessel function is given in explicit form as(4)
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
=
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
v
.
(3)
HT calculation procedure is as follows: based on collocation and explicit formula, we can obtain(5)
F
ν
(
r
)
=
∫
0
∞
λ
f
(
λ
)
J
ν
(
λ
r
)
d
λ
=
∑
i
=
1
N
a
i
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
v
.And it can be expressed as F
=
N
a, where the matrix element is N
j
i
=
(
1
/
ε
i
2
+
r
j
2
)
[
r
j
/
ε
i
2
+
r
j
2
+
ε
i
]
v.It is obvious that the Pravin method is simple and concise. And the calculation error is only caused by the interpolation. Nevertheless, the base functions are monotonically decreasing exponential functions. So it is difficult to approximate convex function, such as approximating the kernel function whose value is zero at the origin. Moreover, to determine the parameterε is also hard in order to ensure the method’s accuracy and stability. Thus, this paper proposed to add a new type base function, which has a convex form and its value is zero at the origin.


## 2.2. Improved Kernel Function Interpolation Theory

Based on the idea of exponential expansion, the new type base function should satisfy the following two conditions. Specifically, the function has convexity and explicit form of HT. After systematically analyzing, we choose the following function as(6)
ϕ
(
λ
)
=
λ
ν
+
1
e
-
ε
λ
2
,
where ν is zero or one, correspondingly to the order of Bessel function. And the improved interpolation formula is as the following:
(7)
G
(
λ
)
=
λ
f
(
λ
)
=
∑
i
=
1
N
a
i
e
-
ε
i
λ
+
∑
i
=
1
N
b
i
λ
ν
+
1
e
-
ε
i
λ
2
.
It can be written as G
=
Ф
w
=
[
Ф
1
,
Ф
2
]
[
a
,
b
]
T, and coefficients w
=
[
a
,
b
]
T
=
[
Ф
1
,
Ф
2
]
-
1
G; the elements are Φ
j
i
1
=
e
-
ε
i
λ
j and Φ
j
i
2
=
λ
j
v
+
1
e
-
ε
i
λ
j
2. So the coefficients a
,
b can be determined by collocation. The HT of (7) can be obtained as
(8)
F
ν
(
r
)
=
∑
i
=
1
N
a
i
∫
0
∞
e
-
ε
i
λ
J
ν
(
λ
r
)
d
λ
+
∑
i
=
1
N
b
i
∫
0
∞
λ
ν
+
1
e
-
ε
i
λ
2
J
ν
(
λ
r
)
d
λ
.
The integrals of the second component can be given by the following formula:
(9)
∫
0
∞
λ
ν
+
1
e
-
ε
λ
2
J
ν
(
λ
r
)
d
λ
=
r
ν
(
2
ε
)
ν
+
1
e
-
(
r
2
/
4
ε
)
.
For v
=
0, referring to formulas (4) and (9), we can obtain zero-order HT as
(10)
F
0
(
r
)
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
+
∑
i
=
1
N
b
i
1
2
ε
i
e
-
(
r
2
/
4
ε
i
)
.
And for v
=
1, we can obtain the first-order HT as
(11)
F
1
(
r
)
=
∑
i
=
1
N
a
i
1
ε
i
2
+
r
2
[
r
ε
i
2
+
r
2
+
ε
i
]
+
∑
i
=
1
N
b
i
r
(
2
ε
i
)
2
e
-
(
r
2
/
4
ε
i
)
.
So the improvement method maintains the simple and concise characteristics. The HT calculation can be obtained by an algebraic summation.


## 3. Numerical Examples

As illustrative numerical examples which possess exact solutions, consider the following functions. Firstly, we carry out the kernel function interpolation simulation which can reveal the approximation ability of improved method. Then, numerical examples of zero-order and first-order HT are implemented, and the method’s accuracy and computational efficiency are analyzed.


### 3.1. Kernel Function Interpolation Simulation

In this section, we focus on the interpolation ability of added base functions, namely,ϕ
1
(
λ
)
=
λ
e
-
ε
λ
2 and ϕ
2
(
λ
)
=
λ
2
e
-
ε
λ
2 which are adopted, respectively, in zero- and first-order HT. Select two typical kernel functions as
(12)
G
1
(
λ
)
=
e
-
2
λ
,
λ
∈
[
0,5
]
,
G
2
(
λ
)
=
λ
2
e
-
λ
2
,
λ
∈
[
0,5
]
.
The first is a monotonic function and the second is a convex one. In numerical calculations, the interpolation nodes are uniformly distributed with interval h
=
0.1 and the test nodes interval h
t
=
0.01. In Pravin’s method, the approximation parameter ε
i
=
0.5
*
(
1.2
)
i
-
1
,
i
=
1,2
,
…
,
51, the first value is 0.5, the common ratio is 1.2, and the last 51th point is 4550.2. So the parameter varies in a wide range with exponential form which maybe results in a poor condition number for interpolation matrix. In this paper, we choose the parameter ε with linear variation as ε
i
=
1.0
+
(
i
-
1
)
h. The added constant is to avoid overflow in (10) or (11). Comparisons between Pravin and the improved method were listed in Tables 1 and 2. The Max_Abs represents the maximum absolute error and RMS means the root mean square error which is defined as
(13)
RMS
=
∑
i
=
1
M
(
abs
_
erro
r
i
)
2
M
.Table 1
Interpolation error comparison between Pravin and improved method(
ε
i
=
0.5
*
(
1.2
)
i
-
1, i
=
1,2
,
…
,
51
).
Method
G
1(λ)
G
2(λ)
Max_Abs
RMS
Max_Abs
RMS
Pravin
4.9469
E
-
3
2.5993
E
-
4
1.3781
E
+
5
7.2755
E
+
3
Φ
1
(
λ
)
2.1563
E
-
3
1.1404
E
-
4
2.9663
E
-
1
2.5646
E
-
2
Φ
2
(
λ
)
1.5147
E
-
3
9.2647
E
-
5
4.9784
E
-
1
3.0318
E
-
2Table 2
Interpolation error comparison between Pravin and improved method(
ε
i
=
1.0
+
(
i
-
1
)
h
,
i
=
1,2
,
…
,
51
).
Method
G
1
(
λ
)
G
2
(
λ
)
Max_Abs
RMS
Max_Abs
RMS
Pravin
6.2607
E
-
5
9.8167
E
-
6
1.0523
E
-
4
1.0746
E
-
5
Φ
1
(
λ
)
5.6133
E
-
6
1.6543
E
-
6
9.3835
E
-
6
2.5233
E
-
6
Φ
2
(
λ
)
2.9375
E
-
6
6.4956
E
-
7
2.8484
E
-
6
1.2230
E
-
6From Tables1 and 2, we can make the following conclusions. Firstly, the Pravin method can obtain good approximation precision for the kernel function G
1
(
λ
) in two kinds of parameter ε settings. Because the G
1
(
λ
) and the interpolation basis e
-
ε
i
λ are all monotonically decreasing functions. However, for function G
2
(
λ
) which is a convex one, the Pravin method is unstable when parameter ε is exponential variation and is still stable for linear variation. And secondly, the improved method not only maintains the stability for both types of functions, but also acquires higher accuracy whose interpolation error is lesser by one or two magnitudes than the Pravin’s. The numerical results showed that the condition number of interpolation matrix decreased dramatically in improved method. For example, the condition number was reduced from 1E38 to 1E17 when parameter ε isexponential variation and was reduced from 1E20 to 1E18 in linear variation. The only deficiency of improved method lies in it being time-consuming, that is, 0.1032 s (for Pravin) and 0.1557 s (for improved method).In addition, the extrapolation ability has been investigated. We choose the region of[
0,3
] for interpolation. The extrapolation test is defined in region [
0,5
]. Numerical results show that the improved method’s maximum absolute errors were 3.6249
E
-
6
(
G
1
(
λ
)
) and 2.6540
E
-
6
(
G
2
(
λ
)
), while the errors of Pravin were 7.9155
E
-
5 and 0.5713 correspondingly. So we can conclude that the improved method is feasible to approximate the kernel function and therefore to calculate HT.


### 3.2. Hankel Transforms Calculation

Here, we have chosen four analytical integrals given by Anderson [15]. They involve the zero- and first-order HT with monotonic and convex types. The analytical integral pairs were
(14)
∫
0
∞
λ
exp
⁡
(
-
λ
2
)
J
0
(
λ
r
)
d
λ
=
1
2
exp
⁡
(
-
r
2
4
)
Integral
1
,
∫
0
∞
exp
⁡
(
-
2
λ
)
J
0
(
λ
r
)
d
λ
=
1
4
+
r
2
Integral
2
,
∫
0
∞
λ
2
exp
⁡
(
-
λ
2
)
J
1
(
λ
r
)
d
λ
=
r
4
exp
⁡
(
-
r
2
4
)
Integral
3
,
∫
0
∞
exp
⁡
(
-
λ
)
J
1
(
λ
r
)
d
λ
=
1
+
r
2
-
1
r
1
+
r
2
Integral
4
.In order to compare with above interpolation simulation easily, all the parameters are chosen the same as Section3.1, only changing the collocation nodes interval h
=
0.25
(
v
=
1
) or 0.2 (v
=
0). Using (7) in the improved method, we can obtain the coefficients and then substitute into (10) or (11) to calculate the unknown values. Relative error comparisons of four HT pairs among Anderson, Pravin, and the improved method proposed in this paper are listed in Tables 3 and 4.Table 3
Relative error comparison among Anderson, Pravin, and improved method for zero-order HT.
r
Integral 1
Integral 2
Anderson
Pravin
Improved
Anderson
Pravin
Improved
1
.
00
E
-
04
1.07
E
-
06
6.38
E
-
06
2.73
E
-
08
2.81
E
-
05
2.76
E
-
10
3.04
E
-
07
1.00
E
-
03
4.77
E
-
07
6.38
E
-
06
2.73
E
-
08
2.62
E
-
06
2.76
E
-
10
3.04
E
-
07
5.00
E
-
03
1.07
E
-
07
6.37
E
-
06
2.73
E
-
08
5.96
E
-
07
2.75
E
-
10
3.04
E
-
07
1.00
E
-
02
1.10
E
-
07
6.33
E
-
06
2.73
E
-
08
7.15
E
-
07
2.73
E
-
10
3.04
E
-
07
5.00
E
-
02
1.19
E
-
07
5.14
E
-
06
2.73
E
-
08
2.38
E
-
07
2.05
E
-
10
3.04
E
-
07
1.00
E
-
01
1.20
E
-
07
2.42
E
-
06
2.73
E
-
08
1.19
E
-
07
5.06
E
-
10
3.04
E
-
07
5.00
E
-
01
1.27
E
-
07
1.09
E
-
06
2.74
E
-
08
0.00
E
+
00
2.79
E
-
10
3.14
E
-
07
1.00
E
+
00
3.83
E
-
07
1.79
E
-
06
2.84
E
-
08
2.67
E
-
07
6.50
E
-
10
3.50
E
-
07
2.00
E
+
00
1.86
E
-
06
3.80
E
-
06
4.08
E
-
08
8.43
E
-
07
7.16
E
-
10
4.65
E
-
07Table 4
Relative error comparison among Anderson, Pravin, and improved method for first-order HT.
r
Integral 3
Integral 4
Anderson
Pravin
Improved
Anderson
Pravin
Improved
1
.
00
E
-
04
2.18
E
-
07
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
3.71
E
-
08
5.11
E
-
08
1.00
E
-
03
0.00
E
+
00
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
4.08
E
-
08
5.1
1
E
-
08
5.00
E
-
03
2.79
E
-
07
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
4.06
E
-
08
5.11
E
-
08
1.00
E
-
02
1.86
E
-
07
4.37
E
-
04
4.60
E
-
09
3.64
E
-
07
4.04
E
-
08
5.11
E
-
08
5.00
E
-
02
2.24
E
-
07
3.80
E
-
04
4.61
E
-
09
3.64
E
-
07
3.50
E
-
08
5.11
E
-
08
1.00
E
-
01
1.49
E
-
07
2.37
E
-
04
4.61
E
-
09
3.64
E
-
07
2.14
E
-
08
5.13
E
-
08
5.00
E
-
01
1.27
E
-
07
5.11
E
-
06
4.86
E
-
09
3.64
E
-
07
3.30
E
-
10
5.62
E
-
08
1.00
E
+
00
2.45
E
-
06
1.39
E
-
06
5.83
E
-
09
3.64
E
-
07
1.14
E
-
10
6.63
E
-
08
2.00
E
+
00
5.91
E
-
06
7.36
E
-
07
1.21
E
-
08
3.64
E
-
07
2.71
E
-
11
7.65
E
-
08From Tables3 and 4, we can conclude that the proposed method can improve the HT accuracy especially for integrals 1 and 3. In above four integrals, the worst accuracy is about the magnitude of 10−7, so the algorithm can obtain good approximation for the theoretical value of zero- or first-order HT. The following Figures 1, 2, 3, and 4 illustrated the absolute error of improved method.Figure 1
Absolute error of Integral 1.Figure 2
Absolute error of Integral 2.Figure 3
Absolute error of Integral 3.Figure 4
Absolute error of Integral 4.Numerical results showed that the maximum absolute errors of the four integrals were1.3639
E
-
8, 1.7500
E
-
7, 5.5460
E
-
9, and 2.2042
E
-
8 separately. And the root mean square errors were corresponding to 8.4464
E
-
9, 1.6560
E
-
7, 3.2052
E
-
9, and 1.5964
E
-
8. The time consuming of above examples was almost 0.05 s.


### 3.3. Application in Electromagnetic Sounding

Considering a four-hierarchical-layer model in DC sounding, the parameters areρ
1
=
3, h
1
=
10, ρ
2
=
30, h
2
=
10, ρ
3
=
1, h
3
=
300, and ρ
4
=
100. The resistivity conversion functionT
(
λ
)
=
T
1
(
λ
) complies with the following recurrence formula as
(15)
T
4
(
λ
)
=
ρ
4
,
i
=
4
,
T
i
(
λ
)
=
ρ
i
ρ
i
(
1
-
e
-
2
λ
h
i
)
+
T
i
+
1
(
1
+
e
-
2
λ
h
i
)
ρ
i
(
1
+
e
-
2
λ
h
i
)
+
T
i
+
1
(
1
-
e
-
2
λ
h
i
)
,
i
=
1,2
,
3
.
The apparent resistivity function is the Hankel transforms of T
(
λ
):
(16)
ρ
s
(
r
)
r
2
=
∫
0
∞
T
(
λ
)
J
1
(
λ
r
)
λ
d
λ
,
T
(
λ
)
=
∫
0
∞
ρ
s
(
r
)
r
2
J
1
(
λ
r
)
r
d
r
.So the resistivity conversion function is a convex one in above parameters; then the Pravin method and the improved method were adopted to approximate theT
(
λ
) and the kernel function λ
T
(
λ
). The interpolation nodes were set in domain [0.01,0.5] with interval h
=
0.01; the test nodes’ interval was h
t
=
0.001. Choosing approximation parameter ε as ε
i
=
15.0
+
(
i
-
1
)
h, the approximation results were shown in Figures 5 and 6.Figure 5
Comparison ofT
(
λ
) approximation.Figure 6
Comparison ofλ
T
(
λ
) approximation.It is obvious that the resistivity conversion function and kernel function can be approximated much better by the improved method. Numerical results inλ
T
(
λ
) approximation showed that the maximum absolute error and root mean square error of approximation were 7.3090
E
-
1 and 1.0405
E
-
1 (for Pravin), while the data were 4.2557
E
-
2 and 3.2775
E
-
3 (for improved method).Finally, the Hankel transforms were performed and compared. The improved method also showed its significant superiority which can be revealed from Figure7. The apparent resistivity function achieved a good approximation on the whole except for the oscillation at the beginning. However, the Pravin method was unstable. The maximum absolute error and root mean square error were 1.8632
E
-
1 and 6.3126 (for Pravin method), and corresponding to improved method, the error was 4.1050 and 8.3735
E
-
1.Figure 7
Comparison of apparent resistivity function by the Hankel transforms.


## 3.1. Kernel Function Interpolation Simulation

In this section, we focus on the interpolation ability of added base functions, namely,ϕ
1
(
λ
)
=
λ
e
-
ε
λ
2 and ϕ
2
(
λ
)
=
λ
2
e
-
ε
λ
2 which are adopted, respectively, in zero- and first-order HT. Select two typical kernel functions as
(12)
G
1
(
λ
)
=
e
-
2
λ
,
λ
∈
[
0,5
]
,
G
2
(
λ
)
=
λ
2
e
-
λ
2
,
λ
∈
[
0,5
]
.
The first is a monotonic function and the second is a convex one. In numerical calculations, the interpolation nodes are uniformly distributed with interval h
=
0.1 and the test nodes interval h
t
=
0.01. In Pravin’s method, the approximation parameter ε
i
=
0.5
*
(
1.2
)
i
-
1
,
i
=
1,2
,
…
,
51, the first value is 0.5, the common ratio is 1.2, and the last 51th point is 4550.2. So the parameter varies in a wide range with exponential form which maybe results in a poor condition number for interpolation matrix. In this paper, we choose the parameter ε with linear variation as ε
i
=
1.0
+
(
i
-
1
)
h. The added constant is to avoid overflow in (10) or (11). Comparisons between Pravin and the improved method were listed in Tables 1 and 2. The Max_Abs represents the maximum absolute error and RMS means the root mean square error which is defined as
(13)
RMS
=
∑
i
=
1
M
(
abs
_
erro
r
i
)
2
M
.Table 1
Interpolation error comparison between Pravin and improved method(
ε
i
=
0.5
*
(
1.2
)
i
-
1, i
=
1,2
,
…
,
51
).
Method
G
1(λ)
G
2(λ)
Max_Abs
RMS
Max_Abs
RMS
Pravin
4.9469
E
-
3
2.5993
E
-
4
1.3781
E
+
5
7.2755
E
+
3
Φ
1
(
λ
)
2.1563
E
-
3
1.1404
E
-
4
2.9663
E
-
1
2.5646
E
-
2
Φ
2
(
λ
)
1.5147
E
-
3
9.2647
E
-
5
4.9784
E
-
1
3.0318
E
-
2Table 2
Interpolation error comparison between Pravin and improved method(
ε
i
=
1.0
+
(
i
-
1
)
h
,
i
=
1,2
,
…
,
51
).
Method
G
1
(
λ
)
G
2
(
λ
)
Max_Abs
RMS
Max_Abs
RMS
Pravin
6.2607
E
-
5
9.8167
E
-
6
1.0523
E
-
4
1.0746
E
-
5
Φ
1
(
λ
)
5.6133
E
-
6
1.6543
E
-
6
9.3835
E
-
6
2.5233
E
-
6
Φ
2
(
λ
)
2.9375
E
-
6
6.4956
E
-
7
2.8484
E
-
6
1.2230
E
-
6From Tables1 and 2, we can make the following conclusions. Firstly, the Pravin method can obtain good approximation precision for the kernel function G
1
(
λ
) in two kinds of parameter ε settings. Because the G
1
(
λ
) and the interpolation basis e
-
ε
i
λ are all monotonically decreasing functions. However, for function G
2
(
λ
) which is a convex one, the Pravin method is unstable when parameter ε is exponential variation and is still stable for linear variation. And secondly, the improved method not only maintains the stability for both types of functions, but also acquires higher accuracy whose interpolation error is lesser by one or two magnitudes than the Pravin’s. The numerical results showed that the condition number of interpolation matrix decreased dramatically in improved method. For example, the condition number was reduced from 1E38 to 1E17 when parameter ε isexponential variation and was reduced from 1E20 to 1E18 in linear variation. The only deficiency of improved method lies in it being time-consuming, that is, 0.1032 s (for Pravin) and 0.1557 s (for improved method).In addition, the extrapolation ability has been investigated. We choose the region of[
0,3
] for interpolation. The extrapolation test is defined in region [
0,5
]. Numerical results show that the improved method’s maximum absolute errors were 3.6249
E
-
6
(
G
1
(
λ
)
) and 2.6540
E
-
6
(
G
2
(
λ
)
), while the errors of Pravin were 7.9155
E
-
5 and 0.5713 correspondingly. So we can conclude that the improved method is feasible to approximate the kernel function and therefore to calculate HT.


## 3.2. Hankel Transforms Calculation

Here, we have chosen four analytical integrals given by Anderson [15]. They involve the zero- and first-order HT with monotonic and convex types. The analytical integral pairs were
(14)
∫
0
∞
λ
exp
⁡
(
-
λ
2
)
J
0
(
λ
r
)
d
λ
=
1
2
exp
⁡
(
-
r
2
4
)
Integral
1
,
∫
0
∞
exp
⁡
(
-
2
λ
)
J
0
(
λ
r
)
d
λ
=
1
4
+
r
2
Integral
2
,
∫
0
∞
λ
2
exp
⁡
(
-
λ
2
)
J
1
(
λ
r
)
d
λ
=
r
4
exp
⁡
(
-
r
2
4
)
Integral
3
,
∫
0
∞
exp
⁡
(
-
λ
)
J
1
(
λ
r
)
d
λ
=
1
+
r
2
-
1
r
1
+
r
2
Integral
4
.In order to compare with above interpolation simulation easily, all the parameters are chosen the same as Section3.1, only changing the collocation nodes interval h
=
0.25
(
v
=
1
) or 0.2 (v
=
0). Using (7) in the improved method, we can obtain the coefficients and then substitute into (10) or (11) to calculate the unknown values. Relative error comparisons of four HT pairs among Anderson, Pravin, and the improved method proposed in this paper are listed in Tables 3 and 4.Table 3
Relative error comparison among Anderson, Pravin, and improved method for zero-order HT.
r
Integral 1
Integral 2
Anderson
Pravin
Improved
Anderson
Pravin
Improved
1
.
00
E
-
04
1.07
E
-
06
6.38
E
-
06
2.73
E
-
08
2.81
E
-
05
2.76
E
-
10
3.04
E
-
07
1.00
E
-
03
4.77
E
-
07
6.38
E
-
06
2.73
E
-
08
2.62
E
-
06
2.76
E
-
10
3.04
E
-
07
5.00
E
-
03
1.07
E
-
07
6.37
E
-
06
2.73
E
-
08
5.96
E
-
07
2.75
E
-
10
3.04
E
-
07
1.00
E
-
02
1.10
E
-
07
6.33
E
-
06
2.73
E
-
08
7.15
E
-
07
2.73
E
-
10
3.04
E
-
07
5.00
E
-
02
1.19
E
-
07
5.14
E
-
06
2.73
E
-
08
2.38
E
-
07
2.05
E
-
10
3.04
E
-
07
1.00
E
-
01
1.20
E
-
07
2.42
E
-
06
2.73
E
-
08
1.19
E
-
07
5.06
E
-
10
3.04
E
-
07
5.00
E
-
01
1.27
E
-
07
1.09
E
-
06
2.74
E
-
08
0.00
E
+
00
2.79
E
-
10
3.14
E
-
07
1.00
E
+
00
3.83
E
-
07
1.79
E
-
06
2.84
E
-
08
2.67
E
-
07
6.50
E
-
10
3.50
E
-
07
2.00
E
+
00
1.86
E
-
06
3.80
E
-
06
4.08
E
-
08
8.43
E
-
07
7.16
E
-
10
4.65
E
-
07Table 4
Relative error comparison among Anderson, Pravin, and improved method for first-order HT.
r
Integral 3
Integral 4
Anderson
Pravin
Improved
Anderson
Pravin
Improved
1
.
00
E
-
04
2.18
E
-
07
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
3.71
E
-
08
5.11
E
-
08
1.00
E
-
03
0.00
E
+
00
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
4.08
E
-
08
5.1
1
E
-
08
5.00
E
-
03
2.79
E
-
07
4.39
E
-
04
4.60
E
-
09
3.64
E
-
07
4.06
E
-
08
5.11
E
-
08
1.00
E
-
02
1.86
E
-
07
4.37
E
-
04
4.60
E
-
09
3.64
E
-
07
4.04
E
-
08
5.11
E
-
08
5.00
E
-
02
2.24
E
-
07
3.80
E
-
04
4.61
E
-
09
3.64
E
-
07
3.50
E
-
08
5.11
E
-
08
1.00
E
-
01
1.49
E
-
07
2.37
E
-
04
4.61
E
-
09
3.64
E
-
07
2.14
E
-
08
5.13
E
-
08
5.00
E
-
01
1.27
E
-
07
5.11
E
-
06
4.86
E
-
09
3.64
E
-
07
3.30
E
-
10
5.62
E
-
08
1.00
E
+
00
2.45
E
-
06
1.39
E
-
06
5.83
E
-
09
3.64
E
-
07
1.14
E
-
10
6.63
E
-
08
2.00
E
+
00
5.91
E
-
06
7.36
E
-
07
1.21
E
-
08
3.64
E
-
07
2.71
E
-
11
7.65
E
-
08From Tables3 and 4, we can conclude that the proposed method can improve the HT accuracy especially for integrals 1 and 3. In above four integrals, the worst accuracy is about the magnitude of 10−7, so the algorithm can obtain good approximation for the theoretical value of zero- or first-order HT. The following Figures 1, 2, 3, and 4 illustrated the absolute error of improved method.Figure 1
Absolute error of Integral 1.Figure 2
Absolute error of Integral 2.Figure 3
Absolute error of Integral 3.Figure 4
Absolute error of Integral 4.Numerical results showed that the maximum absolute errors of the four integrals were1.3639
E
-
8, 1.7500
E
-
7, 5.5460
E
-
9, and 2.2042
E
-
8 separately. And the root mean square errors were corresponding to 8.4464
E
-
9, 1.6560
E
-
7, 3.2052
E
-
9, and 1.5964
E
-
8. The time consuming of above examples was almost 0.05 s.


## 3.3. Application in Electromagnetic Sounding

Considering a four-hierarchical-layer model in DC sounding, the parameters areρ
1
=
3, h
1
=
10, ρ
2
=
30, h
2
=
10, ρ
3
=
1, h
3
=
300, and ρ
4
=
100. The resistivity conversion functionT
(
λ
)
=
T
1
(
λ
) complies with the following recurrence formula as
(15)
T
4
(
λ
)
=
ρ
4
,
i
=
4
,
T
i
(
λ
)
=
ρ
i
ρ
i
(
1
-
e
-
2
λ
h
i
)
+
T
i
+
1
(
1
+
e
-
2
λ
h
i
)
ρ
i
(
1
+
e
-
2
λ
h
i
)
+
T
i
+
1
(
1
-
e
-
2
λ
h
i
)
,
i
=
1,2
,
3
.
The apparent resistivity function is the Hankel transforms of T
(
λ
):
(16)
ρ
s
(
r
)
r
2
=
∫
0
∞
T
(
λ
)
J
1
(
λ
r
)
λ
d
λ
,
T
(
λ
)
=
∫
0
∞
ρ
s
(
r
)
r
2
J
1
(
λ
r
)
r
d
r
.So the resistivity conversion function is a convex one in above parameters; then the Pravin method and the improved method were adopted to approximate theT
(
λ
) and the kernel function λ
T
(
λ
). The interpolation nodes were set in domain [0.01,0.5] with interval h
=
0.01; the test nodes’ interval was h
t
=
0.001. Choosing approximation parameter ε as ε
i
=
15.0
+
(
i
-
1
)
h, the approximation results were shown in Figures 5 and 6.Figure 5
Comparison ofT
(
λ
) approximation.Figure 6
Comparison ofλ
T
(
λ
) approximation.It is obvious that the resistivity conversion function and kernel function can be approximated much better by the improved method. Numerical results inλ
T
(
λ
) approximation showed that the maximum absolute error and root mean square error of approximation were 7.3090
E
-
1 and 1.0405
E
-
1 (for Pravin), while the data were 4.2557
E
-
2 and 3.2775
E
-
3 (for improved method).Finally, the Hankel transforms were performed and compared. The improved method also showed its significant superiority which can be revealed from Figure7. The apparent resistivity function achieved a good approximation on the whole except for the oscillation at the beginning. However, the Pravin method was unstable. The maximum absolute error and root mean square error were 1.8632
E
-
1 and 6.3126 (for Pravin method), and corresponding to improved method, the error was 4.1050 and 8.3735
E
-
1.Figure 7
Comparison of apparent resistivity function by the Hankel transforms.


## 4. Conclusions

Based on the decomposition of kernel function with exponential base functions in Pravin method, this thesis raised the problem of using monotonically decreasing functions to approximate the convex functions. Hence, we proposed an improved scheme by adding a new base function in interpolation which has convexity and explicit form of HT. The improvement method maintains the simple and concise characteristics that convert the Hankel integral to algebraic calculation without any additional numerical processing. The kernel function approximation simulation showed that the improved method acquires higher accuracy whose interpolation error is lesser by one or two magnitudes than the Pravin’s. And the proposed linear variation scheme of parameterε is much more stable than exponential variation scheme. Nevertheless, how to choose the constant and interval h in linear variation scheme is still a discussion problem. In DC sounding example, the improved method also has good precision and high speed even with fewer nodes. So it can be applied to zero-order and first-order HT computation.



---

*Source: 105469-2014-05-05.xml*