Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis

Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis✩ Hong Wang, Hongchao Kang

∗

Department of Mathematics, School of Science, Hangzhou Dianzi University, Hangzhou, Zhejiang 310018, PR China

article

info

Article history: Received 19 July 2019 Received in revised form 5 November 2019 MSC: 65D30 65D32 65R10 41A60 Keywords: Bessel function Singularly oscillatory integrals Complex integration method Modified complex integration method Error analysis

a b s t r a c t This work introduces and analyzes both the complex integration method and the modified complex integration method for two classes of highly oscillatory Bessel transforms with algebraic or algebraic-logarithmic singularities. We first transform the original integrals into the infinite integrals on [1, +∞) by the change of variable t = x−β , where β > 0. With the help of the integral representation of the Bessel function, the integral representation and the asymptotic series of the Whittaker function, we then use the Cauchy’s residue theorem and convert the above infinite integrals into the infinite integrals on [0, +∞). The integrands of the resulting infinite integrals are not oscillatory and decay exponentially fast. Fortunately, they can be efficiently evaluated by using the Gauss-Laguerre quadrature rule. Furthermore, we conduct their error analysis in inverse powers of the frequency ω. Finally, numerical experiments are provided to confirm our theoretical analysis. We can see from theoretical analysis and numerical experiments that the accuracy of the quadrature increases both for the case of growing the number of nodes n and fixed ω, and for the case of an increasing ω and fixed the number of nodes n. The comparison of the proposed two methods and the Filon-type method presented in Kang and Ma (2017), is briefly discussed as well. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Highly oscillatory Bessel transforms arise frequently in many areas such as astronomy, electromagnetic acoustics, scattering problems, physical optics, electrodynamics, seismology image processing and applied mathematics [1–6]. In this paper, we design and analyze efficient and affordable quadrature rules of the highly oscillatory singular Bessel transforms 1

∫ I1 [f ] =

xα f (x)Jm (ω/xβ )dx,

(1.1)

xα ln(x)f (x)Jm (ω/xβ )dx,

(1.2)

0

and 1

∫ I2 [f ] = 0

✩ This research was supported by Zhejiang Provincial Natural Science Foundation of China under Grant Nos. LY18A010009, LZ14A010003, LY17A010029, National Natural Science Foundation of China (Grant Nos. 11301125, 11971138, 11571087, 11171083, 11447005, 11401150), Research Foundation of Hangzhou Dianzi University, PR China (Grant No. KYS075613017). ∗ Corresponding author. E-mail addresses: [email protected] (H. Wang), [email protected], [email protected] (H. Kang). https://doi.org/10.1016/j.cam.2019.112604 0377-0427/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

2

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

where β > 0, α + β > −1 are real numbers, f (x) is analytic in a simply connected and sufficiently (infinitely) large complex region that contains the interval [0, 1], Jm (z) is the Bessel function of the first kind and of order m, the frequency parameter ω is large. Particularly, it should be noticed that transforms (1.1) and (1.2) are integrals with singularities of algebraic or algebraic-logarithmic type, and oscillatory Bessel kernel functions, respectively. In most of cases, such integrals cannot be done analytically and one has to resort to numerical methods. The integrands become highly oscillatory so that it is difficult to calculate by the standard methods when the frequency ω is large. The singularities of algebraic or algebraic-logarithmic type and possible high oscillations of the integrands in (1.1) and (1.2) make the above integrals very difficult to approximate accurately using standard methods, e.g., Gaussian quadrature rules. It is well known that a prohibitively large number of quadrature points is needed if one uses a classic rule such as Gaussian quadratures, or any quadrature method based on (piecewise) polynomial interpolation of the integrands. ∫b In the last few decades, many methods have been developed for calculating the integrals of the form a f (x)Jm (ωg(x))dx without singularity (see [7–17]). Here, it should be noted that Wang in [17] provided asymptotic analysis and gave the Filon-type method for∫computing it. Chen in [7] and [8], presented efficient complex integration methods to calculate the b integrals of the form a f (x)Jm (ωx)dx. In recent years, much progress (the Filon-type method [18], the Clenshaw–Curtis– Filon method [18,19], the orthogonal polynomial ∫ 1 expansion method [20]) has been done in developing numerical schemes for approximating the integrals of the type 0 xα (1 − x)β f (x)Jm (ωx)dx, α, β > −1, with singularities at the two endpoints. Wang, Zhang and Huybrechs in [21] investigated the asymptotics and fast computation of one-sided oscillatory infinite Hilbert transforms. Xu, Xiang and He in [22] studied the fast computation of a class of oscillatory infinite Bessel Hilbert transform. The articles [23,24] devised the Clenshaw–Curtis–Filon method for calculating the oscillatory Airy integrals ∫1 α β x (1 − x) f (x)Ai(−ωx)dx, α, β > −1, with singularities at the two endpoints, and the highly oscillatory finite Hankel 0

∫1

(1)

transform 0 f (x)Hm (ωx)dx, respectively. Kang and Ling in [25] also designed the Clenshaw–Curtis–Filon method for approximating many integrals including different singular oscillatory kernel functions. In fact, Milovanović [26] in 1998, Huybrechs and Vandewalle [27] in 2006, proposed the above-mentioned complex integration method for computing b

∫

f (x)eiωg(x) dx,

(1.3)

a

where f is analytic in a sufficiently large complex region that contains [a, b]. By analytic continuation, the integral (1.3) can be transformed into the two infinite integrals on [0, +∞), where the integrands are not oscillatory and decay exponentially fast. The resulting two infinite integrals can be efficiently computed by using the Gauss–Laguerre quadrature rule. In the following, we would like to mention some related articles [28–34]. The authors of the valuable references [28–31] developed several methods for computing the singular oscillatory integrals of the type 1

∫

β

xα f (x)eiω/x dx,

(1.4)

0

where ω, β > 0, α + β > −1. Gautschi in [29] proposed the modified Chebyshev algorithm for computing the integral (1.4) for the case α = 0, β = 1 and ω = 1. In [30,31], based on the modified Chebyshev algorithm, Hascelik developed an appropriate Gauss quadrature rule to compute the integrals (1.4) for the types of that α = 0, 0 < ω ⩽ 100 and 0.00001 ⩽ β ⩽ 10000. Unfortunately, the methods in [29–31] required the use of high precision arithmetic and the complexity of the modified Chebyshev algorithm in terms of arithmetic operations is O(n2 ) for a n-point Gauss rule. Furthermore, the proposed Gauss quadrature rules in [30,31] were unstable for large ω. Later, Kang and Xiang in [32] devised a numerical steepest descent method for computing the integrals of the form (1.4) for α = 0. Actually, if f is analytical in the complex region G1 containing [0, 1], the numerical steepest descent method can also be applied to (1.4) when α + β > −1. Both the asymptotic method and the Filon-type method were proposed for approximating the integral (1.4) in [28]. Some general methods in [33] were developed by Kang and Xiang, for just sufficiently smooth f on [0, 1], β which relax the strict requirement that f is analytic in G1 in [32]. Further, when the factor eiω/x is changed into the case of β the Bessel function Jm (ω/x ), the integral (1.4) can be seen as the integral (1.1) considered in this paper. For the integrals (1.1) and (1.2), Kang and Ma [34] in 2017 proposed efficient quadrature rules, such as the Filon-type method and the Clenshaw–Curtis–Filon-type method. Their error behaves asymptotically as O(ω−s−2 ), where s is the order of derivatives of the non-oscillatory and smooth part f (x) at the endpoint 1. They can achieve the higher accuracy by either adding derivatives at the endpoint 1 or increasing more interpolation nodes. However, its derivatives are sometimes difficult to Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

3

compute for the complicated function f (x). In addition, only increasing the number of node points cannot improve the error order O(ω−s−2 ). This motivates us to explore more efficient methods so that the obtained error order can be improved greatly by growing the number of node points n. Consequently, our purpose of this paper is to introduce and analyze the efficient and affordable quadrature rules for computing the integrals (1.1) and (1.2). Fortunately, we can devise the two 3 efficient methods in Sections 2–3, whose absolute errors satisfy O(ω−2n− 2 ) or O(ω−λ ), where λ = min{2n + 23 , N + 52 } and N is a positive integer. This is a big advantage that the convergence rate can be speeded up greatly when the number of node points n increases. Now, we summarize the main steps of solution as follows. By change of variable t = x−β , we first transform (1.1) and (1.2) into integrals of the forms I1 [f ] =

1

+∞

∫

β

G(t) t

1

Jm (ωt)dt ,

1 1+ α+ β

(1.5)

and I2 [f ] = −

β

+∞

∫

1 2

G(t) t

1

1 1+ α+ β

ln(t)Jm (ωt)dt ,

(1.6)

−1

where G(t) = f (t β ) is analytic in a simply connected and sufficiently (infinitely) large complex region G = {z ∈ C|1 ⩽ ℜ(z) < +∞, −∞ < ℑ(z) < +∞}. For the convenience of theoretical analysis, we restrict α, β to β > 0, α + β > −1. However, the proposed numerical methods may be available for some special cases of β > 0, α+β ⩽ −1, which are shown in Remark 1 from Section 4. The outline of next sections is as follows. In Section 2 we design the complex integration method and conduct its error analysis. Section 3 presents the modified complex integration method and its error analysis is also discussed here. The above theoretical results and error analyses are verified by numerical experiments in Section 4. Both the methods share the advantageous property that the convergence rate improves greatly as the number of node points n increases for fixed ω. Additionally, the accuracy improves greatly as ω increases for fixed n. Finally, we conclude this paper with some final remarks in Section 5. 2. Complex integration method Before we begin, we mention the preliminaries which will be used below. From [35], the Bessel function Jm (ωt) can be written in the form Jm (ωt) = √

( ω2t )m

1

∫

π Γ (m +

1 ) 2

1

(1 − s2 )m− 2 eiωts ds,

(2.1)

−1

where Γ (z) is the Gamma function. The Whittaker function WhittakerW (k, v, t) in [36, p. 340] is defined by the integral representation t

WhittakerW(k, v, t) =

e− 2 t k

Γ

( 21

+∞

∫

− k + v)

1

z −k− 2 +v (1 + 0

z t

1

)k− 2 +v e−z dz .

(2.2)

From [37, pp. 505 and 508], we know that for −π < arg(t) ⩽ π ,

{ WhittakerW(k, v, t) =

t t k e− 2

1+

N ∑

(−1)

n

n!t n

n=1 1

+O(|t |k−N −1 e− 2 ℜ(t) ),

}

( 12 − k + v )n ( 21 − k − v )n

|t | → ∞.

(2.3)

We can calculate the value of WhittakerW(a, b, z) by either invoking the built-in function ‘whittakerW(a, b, z)’ in Matlab 2017, or truncating the formula (2.3) with a few terms. Next, we start with the derivation of the complex integration method. −1

Theorem 2.1. Suppose that G(z) = f (z β ) is analytic in a simply connected and sufficiently (infinitely) large complex region G = {z ∈ C|1 ⩽ ℜ(z) < +∞, −∞ < ℑ(z) < +∞}. And β > 0, α + β > −1 are real numbers. m is an arbitrary natural number. Then I1 [f ] =

2m

√

m ∑

1

ω β π Γ (m + m

+(−1)j eiω

+∞

∫ 0

1 ) 2 j=0

(2iω)

j

( ){ m j

e−iω

G(1 + ip)T2 (m, j, ω(1 + ip)) (1 + ip)

1 m−j+1+ α+ β

+∞

∫

G(1 − ip)T1 (m, j, ω(1 − ip)) 1 m−j+1+ α+ β

e−ωp dp

(1 − ip)

0

} e−ωp dp ,

(2.4)

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

4

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

I2 [ f ] = −

2m ωm β

{ e−iω

m ∑

π Γ (m + 12 )

j=0

(2iω)j

(1 − ip)

0

m j

1 m−j+1+ α+ β

e−ωp dp

G(1 + ip) ln(1 + ip)T2 (m, j, ω(1 + ip))

+∞

∫

j iω

( )

G(1 − ip) ln(1 − ip)T1 (m, j, ω(1 − ip))

+∞

∫

1

√ 2

+(−1) e

(1 + ip)

0

1 m−j+1+ α+ β

} e

−ωp

dp ,

(2.5)

where T1 (m, j, ωt) and T2 (m, j, ωt) can be written as T1 (m, j, ωt) =

(2ωti)

2m−j−1 2

eωti Γ (

1

2m − j 2m − j

, , 2ωti), 2 2 2m−j−1 1 2m − j 2m − j T2 (m, j, ωt) = (−2ωti) 2 e−ωti Γ ( + 2m − j)WhittakerW(− , , −2ωti). 2 2 2 + 2m − j)WhittakerW(−

2

(2.6) (2.7)

Proof. Putting (2.1) into (1.5) and (1.6) instead of Jm (ωt) we obtain I1 [ f ] =

1

+∞

∫

β

G(t) t

1

√

1 1+ α+ β

( ω2t )m

1

[∫

π Γ (m + 21 )

]

1

(1 − s2 )m− 2 eiωts ds dt ,

(2.8)

−1

and I2 [ f ] = −

1

+∞

∫

β2

G(t) t

1

ln(t) √

1 1+ α+ β

( ω2t )m

1

[∫

π Γ (m + 12 )

]

1

(1 − s2 )m− 2 eiωts ds dt .

(2.9)

−1

By constructing contour integrals, it follows directly from the Cauchy’s residue theorem [38,39] that

∫

1

ie−iωt

1

(1 − s2 )m− 2 eiωts ds =

(ωt)2m

−1

−

]m− 21

u2 + 2iωtu

[

e−u du

0

ieiωt (ω

+∞

∫

+∞

∫

t)2m

u2 − 2iωtu

[

]m− 21

e−u du.

(2.10)

0

Then, using the Binomial theorem, the formula (2.10) can be rewritten as

∫

1

2 m− 21 iωts

(1 − s )

e

( ){ m ∑ (2iω)j m

i

ds =

(ω)2m

−1

j

t 2m−j

j=0 j iωt

+∞

∫

−(−1) e

√

e−iωt

u2m−j

√

u2 + 2iωtu

0

u2m−j

−u

u2 − 2iωtu

0

+∞

∫

e

e−u du

}

du ,

(2.11)

where m is an arbitrary natural number. Now, assume that T1 (m, j, ωt) =

+∞

∫

√ 0

u2m−j u2

+ 2iωtu

e−u du,

(2.12)

e−u du.

(2.13)

and T2 (m, j, ωt) =

+∞

∫

√ 0

u2m−j u2 − 2iωtu

Based on the integral form (2.2) of the WhittakerW (k, v, t), the formulae (2.12) and (2.13) can be transformed into (2.6) and (2.7), respectively. Substituting (2.11) into (2.8) yields that I1 [ f ] =

2m ωm

{∫

√

+∞

πβ Γ (m + 12 ) −iωt ∫ G(t)e

t

1 j

m ∑

i

1 m−j+1+ α+ β

−(−1)

1

t

( ) m j

j=0

+∞

√ 0

G(t)eiωt

+∞

∫

(2iω)j

1 m−j+1+ α+ β

u2m−j u2 + 2iωtu +∞

∫

√ 0

e−u dudt

u2m−j u2 − 2iωtu

} −u

e

dudt

.

(2.14)

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

5

Fig. 1. The integration paths for the oscillator g(t) = t in eiωt (left figure) and the oscillator g(t) = −t in e−iωt (right figure).

Combining (2.12)–(2.14), it is convenient to get I1 [f ] =

√

m ∑

i

(2iω)

j

( ) {∫ m j

ω β π Γ (m + } ∫ +∞ iω t G(t)e T (m , j , ω t) 2 dt . −(−1)j α+1 2m

m

1 ) 2 j=0

t

1

+∞

G(t)e−iωt T1 (m, j, ωt) t

1

1 m−j+1+ α+ β

dt

(2.15)

m−j+1+ β

−1

Since G(z) = f (z β ) is analytic in a simply connected and sufficiently (infinitely) large complex region G = {z ∈ C|1 ⩽ iωz z) ℜ(z) < +∞, −∞ < ℑ(z) < +∞}, the integrand G(z)e T2 (mα+,j,ω is analytic in the region bounded by the closed integration 1 paths h1

⋃

h2

1+A

∫

⋃

z

m−j+1+

β

h3 (see Fig. 1 (left figure)). By the Cauchy residue theorem [38,39], we have

G(t)eiωt T2 (m, j, ωt)

G(z)eiωz T2 (m, j, ωz)

∫

dz , (2.16) 1 m−j+1+ α+ β z where the integration path D includes two parts h2 and h3 defined by the following parametric forms h2 (p) = 1 + ip and h3 (θ ) = 1 + Aeiθ , where p ∈ [0, A], θ : π2 → 0 (see Fig. 1 (left figure)). Then we acquire this series of key formulae 1 m−j+1+ α+ β t

1

1+A

∫

dt =

D

G(t)eiωt T2 (m, j, ωt)

dt 1 m−j+1+ α+ β t ∫ G(z)eiωz T2 (m, j, ωz) dz = 1 m−j+1+ α+ β h2 +h3 z ∫ A G(1 + ip)eiω(1+ip) T2 (m, j, ω(1 + ip)) =i dp 1 m−j+1+ α+ β 0 (1 + ip) 1

π

∫

2

−

iθ

iAeiθ G(1 + Aeiθ )eiω(1+Ae ) T2 (m, j, ω(1 + Aeiθ ))

0

1 m−j+1+ α+ β

(1 + Aeiθ )

dθ.

(2.17)

From (2.3), we have 1

WhittakerW(k, v, t) = O(t k e− 2 t ), |t | → ∞.

(2.18)

Therefore, we can obtain 1

|WhittakerW(k, v, t)| ⩽ K |t |k e− 2 ℜ(t) , |t | → ∞,

(2.19)

where K is some positive constant. From (2.7) and (2.19), it follows that

⏐ ⏐ 2m−j−1 1 ⏐ iω(1+Aeiθ ) ⏐ T2 (m, j, ω(1 + Aeiθ ))⏐ = Γ ( + 2m − j)(2ω|1 + Aeiθ |) 2 ⏐e ⏐2 ⏐ ⏐ ⏐ 2m − j 2m − j iθ ⏐ × ⏐WhittakerW(− , , −2ω(1 + Ae )i)⏐⏐ 2

⩽ KΓ (

1 2

2

1

+ 2m − j)(2ω|1 + Aeiθ |)− 2 e−ωA sin θ .

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

6

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Since G(z) = f (z

− β1

), β > 0, by letting α be argument of the complex number 1 + Aeiθ , we obtain

lim G(1 + Aeiθ ) = lim f ((A2 + 2A cos θ + 1)

A→+∞

− 21β − iβα

e

A→+∞

) = f (0).

This implies that G(1 + Aeiθ ) is bounded, i.e., |G(1 + Aeiθ )| ⩽ C , where C is some positive constant. Based on the inequality [38] sin θ ⩾

2

π

θ,

0⩽θ ⩽

π 2

,

we have

⏐∫ π ⏐ ⏐ 2 iAeiθ G(1 + Aeiθ )eiω(1+Aeiθ ) T (m, j, ω(1 + Aeiθ )) ⏐ ⏐ ⏐ 2 dθ ⏐ ⏐ 1 m−j+1+ α+ ⏐ 0 ⏐ i θ β (1 + Ae ) ⩽

CAK Γ ( 12 + 2m − j) 1 1 m−j+1+ α+ β +2

1 2

(2ω) |A − 1|

CK π Γ ( 12 + 2m − j)

=

1 1 m−j+1+ α+ β +2

3

(2ω) 2 |A − 1|

π

∫

2

2

e − π ω Aθ d θ

0

(1 − e−ωA )

→ 0,

(2.20)

when A → +∞. Therefore, by letting A → +∞ in (2.17), we obtain +∞

∫

G(t)eiωt T2 (m, j, ωt) 1 m−j+1+ α+ β t

1

G(1 + ip)eiω(1+ip) T2 (m, j, ω(1 + ip))

+∞

∫ dt = i

(1 + ip)

0

1 m−j+1+ α+ β

dp.

(2.21)

−1

Since G(z) = f (z β ) is analytic in a simply connected and sufficiently (infinitely) large complex region G = {z ∈ G(z)e−iωz T1 (m,j,ωz) C|1 ⩽ ℜ(z) < +∞, −∞ < ℑ(z) < +∞}, the integrand is analytic in the region bounded by the closed α+1 z

⋃ ⋃

m−j+1+

β

integration paths l1 l2 l3 (see Fig. 1 (right figure)). Here, the integration paths l2 and l3 are defined by the following parametric forms l2 (p) = 1 − ip and l3 (θ ) = 1 + Aeiθ , where p ∈ [0, A], θ : − π2 → 0 in Fig. 1 (right figure). By the similar proof of Eq. (2.21), we get +∞

∫

G(t)e−iωt T1 (m, j, ωt)

1

t

1 m−j+1+ α+ β

+∞

∫ dt = −i

G(1 − ip)e−iω(1−ip) T1 (m, j, ω(1 − ip)) (1 − ip)

0

1 m−j+1+ α+ β

dp.

(2.22)

Thus, combining (2.15), (2.21) with (2.22), it is clear to derive the desired result (2.4). By the similar procedure as used above, the proof of the formula (2.5) can be also attained. ■ Further, we can obtain by letting r = ωp in (2.4),

⎧ ( ) ⎨ −iω ∫ +∞ G(1 − irω )T1 (m, j, ω(1 − irω )) −r 1 e m I1 [ f ] = e dr (2iω)j √ 1 1 j ⎩ ω m−j+1+ α+ 2m ωm β π Γ (m + 2 ) j=0 β 0 (1 − irω ) ⎫ ∫ +∞ ir ir ⎬ iω G(1 + )T (m , j , ω (1 + )) e 2 ω ω +(−1)j e−r dr . α+ 1 m−j+1+ β ⎭ ω 0 (1 + ir ) m ∑

(2.23)

ω

By the Gauss–Laguerre quadrature rule with n nodes ts and weights as , we have the following n-point quadrature formula for (1.1), QnC [G] =

1

2m

ω

m+1

√

β π Γ (m +

m ∑ 1 ) 2 j=0

(2iω)j

⎧ ( )⎨

n

∑ G(1 − its )T1 (m, j, ω(1 − m ω e−iω as 1 j ⎩ m−j+1+ α+ β (1 − its ) s=1

⎫ n ∑ G(1 + itωs )T2 (m, j, ω(1 + itωs )) ⎬ j iω +(−1) e as . 1 m−j+1+ α+ ⎭ β (1 + its ) s=1

its

ω

))

ω

(2.24)

ω

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

7

Similarly, we can get the corresponding n-point quadrature formula for (1.2) as follows,

˜ QnC [G] = − ⎧ ⎨

m ∑

1

2m

ω

e−iω

m+1

n ∑

⎩

(2iω)

√ β π Γ (m + 12 ) 2

its

G(1 −

as

ω

+(−1)j eiω

n ∑

) ln(1 −

as

G(1 +

its

ω

m j

)T1 (m, j, ω(1 −

α+1 its m−j+1+ β )

its

ω

))

ω

its

) ln(1 + ω (1 +

s=1

( )

j=0

(1 −

s=1

j

)T (m, j, ω(1 + ω 2

its

α+1 its m−j+1+ β

ω

)

its

⎫

)) ⎬ ω

.

(2.25)

⎭

Here, we regard the two n-point quadrature formulae (2.24) and (2.25) as the so-called complex integration method. In the sequel, we give the error analysis of the n-point quadrature formulae (2.24) and (2.25) in inverse powers of ω. Theorem 2.2. Under the same conditions as those of Theorem 2.1, it is true for ω ≫ 1 that

|I1 [f ] − QnC [G]| = O( |I2 [f ] − ˜ QnC [G]| = O(

1

ω

2n+ 23

ω

2n+ 23

1

),

(2.26)

).

(2.27)

Proof. From [40], the error formula for the n-point Gauss–Laguerre quadrature rule of the integral as follows, (n!)2 (2n)!

f (2n) (ξ ),

0 < ξ < +∞.

∫ +∞ 0

f (x)e−x dx is given

(2.28)

Combining (2.23), (2.24) and (2.28) yields that

⏐ ⏐ ⏐I1 [f ] − Q C [G]⏐ n

=

2m β

√

j ∑ (2ω) m

1

π Γ (m + 21 )

+(−1)j

+∞

∫

+(−1)j

⩽

2m β

−

√

n ∑

n ∑

−

s=1

)T (m, j, ω(1 + ω 2 α+1 ir m−j+1+ β

ω

)T (m, j, ω(1 + ω 2 α+1 its m−j+1+ β

(1 +

1

j ∑ (2ω)

π Γ (m + 21 )

as

G(1 −

as

ω

m j

α+1 its m−j+1+ β

its

α+1 its m−j+1+ β

)

its

ω

)T1 (m, j, ω(1 −

(1 −

α+1 its m−j+1+ β

ω

its

ω

))

)

⎤⏐ ⏐ ⎦⏐ ⏐ ⏐

)) ⏐ ω

⏐ ⏐ ⏐ ⏐∫ +∞ ⏐ G(1 + irω )T2 (m, j, ω(1 + irω )) −r ⏐+⏐ e dr 1 ⏐ ⏐ m−j+1+ α+ β ⏐ ⏐ 0 (1 + irω )

)) ⏐ ω

)

)T (m, j, ω(1 + ω 2 ω

⎡ n ∑ G(1 − e−r dr ⎦ − ⎣ as

⎧⏐ ⎨⏐⏐∫ +∞ G(1 − ir )T (m, j, ω(1 − ir )) ω 1 ω ⏐ e−r dr α+1 ⎩⏐⏐ 0 ir m−j+1+ β (1 − ω )

ωm+1

its

(1 +

⎤

( )

)T (m, j, ω(1 − ω 1 ω

its

)

its

(1 − G(1 +

j=0

)) ω

s=1

its

m

ir

)

s=1

s=1

n ∑

as

G(1 +

⏐⎡ ⏐ ∫ +∞ ⏐ G(1 − irω )T1 (m, j, ω(1 − irω )) −r ⏐⎣ e dr 1 ⏐ m−j+1+ α+ β ⏐ 0 (1 − irω )

ir

(1 +

0

m j

ωm+1

j=0

G(1 +

( )

its

⏐⎫ ⏐ ⏐ ⏐⎭ ⏐

)) ⏐⎬ ω

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

8

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

=

1

2m β

( ) ⎧⏐ ⎤(2n) ⏐⏐ m ⏐⎡ ⎪ iξ1 ⎨ ⏐ ⏐ j G(1 − ω ) iξ1 ⏐⎣ ⎦ ⏐⏐ T (m , j , ω (1 − )) ⏐ 1 α+ 1 ⏐ ωm+1 ⎪ ω ⎩⏐⏐ (1 − iξω1 )m−j+1+ β ⏐

j m (2ω ) (n!)2 ∑

√ π Γ (m + 12 ) (2n)!

j=0

⏐⎡ ⎤(2n) ⏐⏐⎫ ⏐ iξ2 ⏐⎪ ⎬ ⏐ G(1 + ω ) iξ2 ⏐ ⏐ ⎦ T (m , j , ω (1 + )) + ⏐⎣ ⏐ , α+ 1 2 ⏐⎪ ⏐ (1 + iξ2 )m−j+1+ β ω ⏐⎭ ⏐ ω where ξ1 , ξ2 ∈ (0, +∞). ξ Letting η1 = ω1 , η2 =

ξ2 ω

(2.29)

in (2.29), it follows that

⏐ ⏐ ⏐I1 [f ] − Q C [G]⏐ n

j m (2ω ) (n!)2 ∑

( ) m j

⎧⏐[ ](2n) ⏐⏐ ⎨⏐⏐ ⏐ G(1 − iη1 ) ⏐ ⏐ T (m, j, ω(1 − iη1 )) ⩽ √ α+1 1 ⏐ 1 2n+m+1 ⎩⏐ m m − j + 1 + (2n) ! ω 2 β π Γ (m + 2 ) β ⏐ ⏐ (1 − i η ) j=0 1 ⎫ ⏐ ⏐[ ] (2n) ⏐ ⏐ ⏐⎬ ⏐ G(1 + iη2 ) ⏐ , + ⏐⏐ T (m , j , ω (1 + i η )) 2 2 α+1 ⏐⎭ ⏐ ⏐ (1 + iη2 )m−j+1+ β 1

(2.30)

where η1 , η2 ∈ (0, +∞). From (2.18), we have WhittakerW(− WhittakerW(−

2m − j 2m − j

, 2ωti) = O(|2ωt |−

,

, −2ωti) = O(|2ωt |−

2 2 2m − j 2m − j 2

2m−j 2

,

2

eωℑ(t) ),

2m−j 2

e−ωℑ(t) ).

(2.31) (2.32)

Combining (2.6), (2.7), (2.31) and (2.32), we obtain 1

T1 (m, j, ωt) = O(|2ωt |− 2 ),

(2.33)

and 1

T2 (m, j, ωt) = O(|2ωt |− 2 ).

(2.34)

Following (2.30), (2.33) and (2.34), we can prove that

⏐ ⏐ ⏐I1 [f ] − Q C [G]⏐ n

⎧⏐ ](2n) ⏐⏐ ( ) ⎨⏐[ ⏐ ⏐ G(1 − i η ) 1 m 1 ⏐ ⏐ ⩽ O(|2ω(1 − iη1 )|− 2 ) √ α+1 ⏐ ⏐ 1 2n + m + 1 m j m − j + 1 + ⎩⏐ (1 − iη ) 2 β π Γ (m + 2 ) (2n)! j=0 ω β ⏐ 1 ⏐[ ⏐ ⎫ ] (2n) ⏐ ⏐ ⏐ ⏐⎬ G(1 + iη2 ) − 21 ⏐ + ⏐⏐ ) O( | 2 ω (1 + i η ) | 2 α+ 1 ⏐⎭ ⏐ (1 + iη2 )m−j+1+ β ⏐ ⎧⏐[ ](2n) ⏐⏐ ⏐ ( ) m ⎨⏐ ⏐ 1 (n!)2 ∑ m G(1 − iη1 ) 1 ⏐ ⏐ = O(|2(1 − iη1 )|− 2 ) √ 1 ⏐ j ⎩⏐ 1 (2n)! 2n+ 32 m−j+1+ α+ β ⏐ ⏐ ω β π Γ (m + 2 ) (1 − i η ) j=0 1 ⏐[ ⏐ ⎫ ] (2n) ⏐ ⏐ ⏐ ⏐⎬ G(1 + iη2 ) − 12 ⏐ + ⏐⏐ O( | 2(1 + i η ) | ) 2 α+1 ⏐⎭ ⏐ (1 + iη2 )m−j+1+ β ⏐ 1

3

= O(ω−2n− 2 ),

m (n!)2 ∑ (2ω)m

as ω → +∞.

Therefore, we complete the proof of (2.26). The proof of (2.27) can be obtained by the analogous procedure used above. ■ Unfortunately, the produced quadrature formulae (2.24) and (2.25) only hold for all natural number m since the transform from (2.10) to (2.11) is based on the Binomial theorem. Here, m is the order of the Bessel function. In the next section, we shall present a modified method to overcome this weakness. Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

9

3. Modified complex integration method In this part, a new method which possesses a wider and more flexible application for all real number m will be proposed. Before we begin, we illustrate the preliminaries which will be used in the following. The Bessel function Jm (t) [37, p. 358] can be expressed by Jm (t) =

] 1 [ (1) (2) (t) , Hm (t) + Hm 2

(1)

(3.1)

(2)

(1)

where the Hm (t) and Hm (t) denote Hankel functions of the first and second kinds, respectively. Moreover, the Hm (t) (2) and Hm (t) can be written by the Whittaker function WhittakerW(k, v, t) [37, pp. 505 and 510] in the forms,

√ (1) (t) = Hm

2

π

√πt 2

(2) Hm (t) =

πt

1

e−i 2 (m+ 2 ) WhittakerW(0, m, −2it), π

(3.2)

1

ei 2 (m+ 2 ) WhittakerW(0, m, 2it).

(3.3)

From (2.3), the asymptotic series of WhittakerW(k, v, t) can be rewritten as,

{ k − 2t

WhittakerW(k, v, t) = t e

+∞ ∑

1+

(−1)

n

( 12 − k + v )n ( 21 − k − v )n n!t n

n=1 t

= t k e− 2 Φk,v (t),

}

|t | → ∞,

(3.4)

where

Φk,v (t) = 1 +

+∞ ∑

(−1)n

( 21 − k + v )n ( 12 − k − v )n n!t n

n=1

.

Letting N ∑

Φk,v,N (t) = 1 +

(−1)n

( 12 − k + v )n ( 12 − k − v )n n!t n

n=1

,

(3.5)

this implies that

Φk,v (t) − Φk,v,N (t) = O(|t |−N −1 ),

|t | → ∞.

(3.6)

Then we give a significant theorem for approximation of integrals (1.1) and (1.2) as follows. −1

Theorem 3.1. Suppose that G(z) = f (z β ) is analytic in a simply connected and sufficiently (infinitely) large complex region G = {z ∈ C|1 ⩽ ℜ(z) < +∞, −∞ < ℑ(z) < +∞}. And β > 0, α + β > −1, m, are real numbers. Then it follows that I1 [f ] = I2 [f ] =

1[ 2 1[ 2

]

MC Q MC (1) [G] + Q (2) [G] , Hm

(3.7)

Hm

]

˜ ˜MC Q MC (1) [G] + Q (2) [G] , Hm

(3.8)

Hm

where Q MC (1) [G] = Hm

1

√

β

Q MC (2) [G] = − Hm

1

β

2

π

πω √

1

e−i 2 (m+ 2 ) ieiω

2

πω √

+∞

∫

G(1 + ip)Φ0,m (−2iω(1 + ip)) (1 + ip)

0 π

1

ei 2 (m+ 2 ) ie−iω

+∞

∫

3 1 + α+ β 2

G(1 − ip)Φ0,m (2iω(1 − ip)) (1 − ip)

0

3 1 + α+ β 2

e−ωp dp,

(3.9)

e−ωp dp,

(3.10)

∫ 1 2 −i π (m+ 1 ) iω +∞ G(1 + ip) ln(1 + ip)Φ0,m (−2iω(1 + ip)) −ωp ˜ 2 ie Q MC e 2 e dp, (1) [G] = − 2 3 Hm + α+1 β πω 0 (1 + ip) 2 β √ ∫ 1 2 i π (m+ 1 ) −iω +∞ G(1 − ip) ln(1 − ip)Φ0,m (2iω(1 − ip)) −ωp MC ˜ 2 2 Q (2) [G] = 2 e ie e dp. 3 Hm + α+1 β πω 0 (1 − ip) 2 β

(3.11)

(3.12)

Proof. From (3.2)–(3.4) we can get +∞

∫ 1

G(t)

(1) Hm ( 1 1+ α+ β t

√ ωt)dt =

2

πω

π

1

e−i 2 (m+ 2 )

+∞

∫ 1

G(t) t

3 1 + α+ β 2

Φ0,m (−2iωt)eiωt dt ,

(3.13)

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

10

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

and +∞

∫

√

G(t)

(2) Hm ( 1 1+ α+ β t

1

ωt)dt =

2

π

πω

+∞

∫

1

ei 2 (m+ 2 )

G(t) t

1

Φ0,m (2iωt)e−iωt dt .

1 3 + α+ β 2

(3.14)

Then using the Cauchy’s residue theorem [38,39] (see Fig. 1), and according to the proof of Theorem 2.1, we have

∫ G(z)z

1 − 32 − α+ β

Φ0,m (−2iωz)eiωz dz → 0,

h3

and

∫ G(z)z

1 − 23 − α+ β

Φ0,m (2iωz)e−iωz dz → 0,

l3

as A → +∞. Therefore, we obtain +∞

∫

G(t) t

1

∫

Φ0,m (−2iωt)eiωt dt =

3 1 + α+ β 2

G(z)z

1 − 23 − α+ β

Φ0,m (−2iωz)eiωz dz

h2

= ieiω

+∞

∫

G(1 + ip)Φ0,m (−2iω(1 + ip))e−ωp 3

(1 + ip) 2

0

1 + α+ β

dp,

(3.15)

and +∞

∫

G(t) t

1

∫

Φ0,m (2iωt)e−iωt dt =

3 1 + α+ β 2

G(z)z

1 − 23 − α+ β

Φ0,m (2iωz)e−iωz dz

l2

= −ie−iω

+∞

∫

G(1 − ip)Φ0,m (2iω(1 − ip))e−ωp 3

(1 − ip) 2

0

1 + α+ β

dp.

(3.16)

Then plugging (3.15) and (3.16) into (3.13) and (3.14) derives (3.9) and (3.10), respectively. The proof of (3.11) and (3.12) is similar to the above procedure. Therefore, we complete the proof. ■ Further, setting r = ωp in (3.9) and (3.10) yields that Q

MC

(1) Hm

[G] =

√

1

β

2

πω

−i π2 (m+ 12 ) ie

e

iω

+∞

∫

ω

G(1 +

ir

ω

)Φ0,m (−2iω(1 +

(1 +

0

ir

) ω

3 1 + α+ β 2

ir

ω

))

e−r dr ,

(3.17)

e−r dr .

(3.18)

and Q MC (2) [G] = − Hm

1

√

β

2

π

πω

1

ei 2 (m+ 2 )

ie−iω

+∞

∫

ω

G(1 −

ir

ω

)Φ0,m (2iω(1 −

(1 −

0

ir

) ω

3 1 + α+ β 2

ir

ω

))

Hence, by using the Gauss–Laguerre quadrature rule with n nodes tk and weights ak in (3.17) and (3.18), we obtain Q

MC

(1) n,Hm

[G] =

1

√

β

2

π

πω

1

e−i 2 (m+ 2 )

n ieiω ∑

ω

ak

G(1 +

itk

ω

)Φ0,m (−2iω(1 +

(1 +

k=1

itk

ω

)

3 1 + α+ β 2

itk

ω

))

,

(3.19)

.

(3.20)

and Q MC (2) [G] = − n,Hm

1

√

β

2

π

πω

n ie−iω ∑

1

ei 2 (m+ 2 )

ω

ak

G(1 −

itk

ω

)Φ0,m (2iω(1 −

(1 −

k=1

itk

3 1 + α+ β 2

)

ω

itk

ω

))

By truncating after the first N + 1 terms of Φ0,m (t) in (3.19) and (3.20), we obtain Q

MC

(1) n,N ,Hm

[G] =

1

√

β

2

π

πω

1

e−i 2 (m+ 2 )

n ieiω ∑

ω

ak

G(1 +

itk

ω

)Φ0,m,N (−2iω(1 + (1 +

k=1

itk

ω

)

3 1 + α+ β 2

itk

ω

))

,

(3.21)

.

(3.22)

and Q

MC

(2) n,N ,Hm

[G] = −

1

β

√

2

πω

π

1

ei 2 (m+ 2 )

n ie−iω ∑

ω

k=1

ak

G(1 −

itk

ω

)Φ0,m,N (2iω(1 −

(1 −

itk

ω

)

3 1 + α+ β 2

itk

ω

))

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

11

Combining (3.7), (3.21) and (3.22), we have the n-point quadrature formula for (1.1) as follows, QnMC ,N [G] =

1 2β

−e

√

⎡

iω

2

⎣e−i π2 (m+ 12 ) ie πω ω

i π2 (m+ 21 )

n ie−iω ∑

ω

ak

n ∑

ak

G(1 +

itk

(1 +

k=1

G(1 −

itk

ω

itk

ω

)Φ0,m,N (2iω(1 −

(1 −

k=1

)Φ0,m,N (−2iω(1 +

ω

itk

ω

3

)2

1 + α+ β

)

3 1 + α+ β 2

itk

ω

))

itk

ω

))

⎤ ⎦.

(3.23)

Similarly, we can get the following n-point quadrature formula for (1.2), 1 ˜ QnMC ,N [G] = − 2β 2

−e

√

⎡

iω

2

⎣e−i π2 (m+ 12 ) ie πω ω

i π2 (m+ 21 )

n ie−iω ∑

ω

ak

k=1

G(1 −

n ∑

ak

itk

G(1 +

ω

) ln(1 +

ω

ω

(1 +

k=1 itk

itk

) ln(1 − (1 −

itk

ω

)Φ0,m,N (−2iω(1 + itk

ω

)

3 1 + α+ β 2

)Φ0,m,N (2iω(1 −

itk

ω

3

)2

1 + α+ β

itk

ω

))

itk

ω

))

⎤ ⎦.

(3.24)

Here we call Eqs. (3.23) and (3.24) as the modified complex integration method, which holds for all real number m. The error analysis of the above method is then described by Theorem 3.2. Theorem 3.2. Under the same conditions as those of Theorem 3.1, it follows for ω ≫ 1 that −λ |I1 [f ] − QnMC ,N [G]| = O(ω ), −λ |I2 [f ] − ˜ QnMC ,N [G]| = O(ω ),

(3.25) (3.26)

where λ = min{2n + 23 , N + 25 }. Proof. Using the formulae (2.28) and (3.6), we can derive an expression for the absolute error that

⏐ ⏐ ⏐ ⏐ MC [ G ] ⏐ ⏐QH (1) [G] − QnMC (1) ,N ,H m m ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ MC ⏐ MC MC ⩽ ⏐Q MC (1) [G]⏐ (1) [G] − Q (1) [G]⏐ + ⏐Q (1) [G] − Q n,N ,Hm n,Hm n,Hm Hm ⎧⏐ √ ir ir ⎨⏐∫ 1 2 1 ⏐ +∞ G(1 + ω )Φ0,m (−2iω(1 + ω )) −r ⏐ = e dr 3 α+ 1 ⏐ + β π ω ω ⎩⏐ 0 (1 + irω ) 2 β ⏐ ⏐ ⏐ ⏐ n n it it it ∑ G(1 + ωk )Φ0,m (−2iω(1 + ωk )) ⏐ ⏐∑ G(1 + ωk )Φ0,m (−2iω(1 + ⏐+⏐ a − ak k 3 3 ⏐ ⏐ + α+1 + α+1 it it ⏐ ⏐ k=1 (1 + ωk ) 2 β (1 + ωk ) 2 β k=1 ⏐⎫ ⏐ n it it ∑ G(1 + ωk )Φ0,m,N (−2iω(1 + ωk )) ⏐⎬ ⏐ − ak 3 ⏐⎭ + α+1 it ⏐ (1 + ωk ) 2 β k=1 ⏐ ⎧ ⏐⎡ ⎤(2n) ⏐ ⏐ √ ⎪ n iξ iξ ⏐ ⎨ 2 ⏐ ∑ G(1 + ω )Φ0,m (−2iω(1 + ω )) 2 1 (n!) ⏐ ⏐ ⎦ = ak ⏐⎣ ⏐ 3 α+ 1 3 iξ 2 + β ⏐ ⏐ β πω ⎪ (2n) ! ⎩ (1 + ω ) ⏐ k=1 ξ ∈(0,+∞) ⏐ ⏐ n ⏐} ⏐∑ ⏐ itk itk −N − 5 − α+1 ⏐ −N −1 ⏐ β 2 +⏐ O((2ω) )⏐ . ak G(1 + )(1 + ) ⏐ ⏐ ω ω

itk

ω

))

(3.27)

k=1

Letting η =

ξ ω

in (3.27), it follows that

⏐ ⏐ ⏐ MC ⏐ ⏐QH (1) [G] − QnMC (1) [G]⏐ ,N ,Hm m ⎧ ⏐[ ⏐ ](2n) √ ⏐ n ⎨ 1 (n!)2 ⏐⏐ ∑ ⏐ 1 2 G(1 + iη)Φ0,m (−2iω(1 + iη)) ⏐ ⏐ ⩽ a k 3 α+1 ⏐ ⏐ 3 2n + β π ω ⎩ ω (2n)! ⏐ β 2 ⏐ (1 + iη) k=1 η∈(0,+∞) ⏐ n ⏐} ⏐∑ ⏐ itk −N − 5 − α+1 itk ⏐ ⏐ β O((2ω )−N −1 )⏐ . 2 +⏐ ak G(1 + )(1 + ) ⏐ ⏐ ω ω

(3.28)

k=1

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

12

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

[ The derivative

Φ0,m (−2iωt) 3 + α+1

t2

[

Φ0,m (−2iωt) t

3 1 + α+ β 2

β

](2q) required in the second to the last line of (3.28) can be written as t =1+iη

](2q) =

∞ ∑ (−1)n+2q ( 21 + m)n ( 12 − m)n

t =1+iη

∏2q−1 j=0

3 2

(n +

+

α+1 β

+ j)

(−2iω)n n!

n=0

(1 + iη)

1 −n− 23 − α+ β −2q

,

⏐[ ](2q) ⏐⏐ ⏐ ⏐ Φ0,m (−2iωt) ⏐ where (a)0 = 1. This implies that ⏐ ⏐ = O(1) as ω → ∞. Therefore, we can directly from (3.28) obtain ⏐ t 32 + α+β 1 t =1+iη ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ MC ⏐ ⏐ MC ⏐ ⏐ MC ⏐ MC MC ⏐QH (1) [G] − QnMC (1) [G]⏐ ⩽ ⏐Q (1) [G] − Q (1) [G]⏐ + ⏐Q (1) [G] − Q (1) [G]⏐ ,N ,Hm Hm n,Hm n,Hm n,N ,Hm m ] 3 [ = ω− 2 O(ω−2n ) + O(ω−N −1 ) , ω → +∞. Similarly, we have

⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ MC ⏐ ⏐ MC ⏐ ⏐ MC ⏐ MC MC ⏐QH (2) [G] − QnMC (2) [G]⏐ ⩽ ⏐Q (2) [G] − Q (2) [G]⏐ + ⏐Q (2) [G] − Q (2) [G]⏐ ,N ,Hm Hm n,Hm n,Hm n,N ,Hm m ] 3 [ = ω− 2 O(ω−2n ) + O(ω−N −1 ) , ω → +∞. Therefore, it is true from (3.7) that

⏐ ⏐ [ ] ⏐I1 [f ] − Q MC [G]⏐ ⩽ ω− 32 O(ω−2n ) + O(ω−N −1 ) n,N 3

5

= O(ω−2n− 2 ) + O(ω−N − 2 ),

ω → +∞.

(3.29)

This implies that the proof of (3.25) is complete. By the similar procedure used above, the proof of (3.26) is also attained. ■ 4. Numerical examples The methods proposed in Sections 2 and 3 are extremely efficient for the calculation of highly oscillatory integrals of the forms (1.1) and (1.2). Their error analyses (2.26), (2.27), (3.25) and (3.26) in Theorems 2.2 and 3.2 demonstrate that the accuracy can be rapidly improved by either adding the node points for fixed ω or increasing ω for fixed n. More particularly, we can see from Theorems 2.2 and 3.2 that the convergence rate would be improved by increasing the number of nodes n for fixed ω. All our numerical experiments are carried out in Matlab 2017. 4.1. Complex integration method Based on the above discussion for the complex integration method, we briefly outline the steps of solution as follows. Method 1 Step 1 Transform the integrals (1.1) and (1.2) by change of variable t = x−β into (1.5) and (1.6), respectively. Step 2 Plugging the integral form (2.1) of the Bessel function into (1.5) and (1.6) yields (2.8) and (2.9). Step 3 Convert the integral form (2.1) into two infinite integrals (2.11) on [0, +∞) by the Cauchy’s residue theorem and the Binomial theorem. Step 4 The two infinite integrals in (2.15) are expressed as (2.21) and (2.22) by means of the Cauchy’s residue theorem. Step 5 Applying the n-point Gauss–Laguerre quadrature rule to (2.23) and the related formula, derives the desired n-point quadrature formulae (2.24) and (2.25) for (1.1) and (1.2). In the following, on the basis of the steps of solution, we show numerical experiments to verify the complex integration method. Example 1. Let us consider displayed in Fig. 2.

∫1 0

1

1

x− 8 cos xJ3 (ω/x 6 )dx by using the complex integration method and numerical results are

Example 2. We focus on the calculation of error is illustrated in Fig. 3.

∫1 0

1

1

x− 3 ln(x) arctan xJ2 (ω/x 3 )dx by the complex integration method and its

∫1

1

1

Example 3. We consider the computation of 0 x− 5 arctan xJ1 (200/x 5 )dx and complex integration method. The numerical results are shown in Fig. 4.

∫1 0

1

1

x− 6 ln(x) sin xJ4 (200/x 6 )dx by the

We can verify Theorems 2.1 and 2.2 by Examples 1–3. Figs. 2–3 show that the complex integration method is very efficient because of only using a few nodes, such as n = 2, 3. We can also observe from Figs. 2–3 that the error order 3 O(ω−2n− 2 ) is consistent with the error analysis of Theorem 2.2. In addition, for fixed number of nodes n, the required Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 2. The absolute error (left figure) and absolute error scaled by ω7.5 (right figure) of Q3C [G], for the integral

∫1 0

13

1

1

x− 8 cos xJ3 (ω/x 6 )dx with n = 3.

accuracy increases rapidly as the frequency ω increases. Moreover, we can see from the left (right) figure of Fig. 4 that the desired accuracy would be improved enormously when adding the nodes but n ⩽ 4 (n ⩽ 6) for fixed ω. When we continue to increase number of nodes n, the complex integration method can achieve machine precision for computing the integrals in Example 3 with ω = 200. 4.2. Modified complex integration method Based on the above discussion for the modified complex integration method, we summarize the main steps of solution as follows. Method 2 Step 1 Transform the integral (1.1) and (1.2) by change of variable t = x−β into (1.5) and (1.6). Step 2 With the help of the key expression (3.1) of the Jm (t) in terms of the Hankel functions, we rewrite (1.5) as two related infinite integrals about the Hankel functions. Step 3 Express the two kinds of the Hankel functions by the Whittaker functions in (3.2) and (3.3), respectively. Substitute the formulae (3.2) and (3.3) into the left sides of Eqs. (3.13) and (3.14), respectively. Step 4 Convert the two infinite integrals on the right sides of (3.13) and (3.14) into (3.15) and (3.16) by using the Cauchy’s residue theorem, respectively. Step 5 Applying the Gauss–Laguerre quadrature rule to the integrals in (3.17) and (3.18) and truncating the first N + 1 terms of Φ0,m (t), derives the n-point quadrature formula (3.23) for (1.1). By the similar process, the n-point quadrature formula (3.24) for (1.2) is also obtained. In the following, on the basis of the steps of solution, we now test numerical examples to illustrate the quality of the approximations provided by the modified complex integration method. Example 4. We apply the modified complex integration method for computing

∫1 0

1

1

x− 5 (1 + x)4 J 1 (ω/x 2 )dx, and the error 3

is presented in Fig. 5. Here, the order of the Bessel function can be chosen as a non-integer m = 13 . In fact, the modified complex integration method is available for all positive real number m. This is a big improvement for the complex integration method, which holds only for all natural number m. Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

14

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 3. The absolute error (left figure) and absolute error scaled by ω5.5 (right figure) of ˜ Q2C [G], for the integral n = 2.

Fig.

∫1 0

4. The − 16

x

absolute

error

of

the

complex

integration

method

for

the

integral

∫1 0

1

∫1 0

1

1

x− 3 ln(x) arctan xJ2 (ω/x 3 )dx with

1

x− 5 arctan xJ1 (200/x 5 )dx

(left

figure)

and

1 6

ln(x) sin xJ4 (200/x )dx (right figure) with n from 2 to 20.

Example 5. In Fig. 6, we illustrate the error of the integral

∫1 0

1

1

x− 2 sin x ln(x)J6 (ω/x 5 )dx by using the modified complex

integration method. Example 6. Let us consider

∫1 0

1

1

x− 5 sin xJ3 (ω/x 2 )dx by using the modified complex integration method and show the

absolute error and scaled absolute error in Fig. 7. Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 5. The absolute error (left figure) and absolute error scaled by ω5.5 (right figure) of Q2MC [G], for the integral

∫1 0

15

1

1

x− 5 (1 + x)4 J 1 (ω/x 2 )dx with 3

n = 2 and N = 2n − 1.

∫1

1

1

2

∫1

5

1

Example 7. We consider the computation of 0 x− 3 sinexx J 1 (500/x 4 )dx and 0 x− 6 ln(x) arctan x3 J2 (500/x 2 )dx by the 5 modified complex integration method. The numerical results are illustrated in Fig. 8. Figs. 5–8 show the high accuracy and efficiency of the modified complex integration method. In Examples 4–5, the 3 5 decay rate of their absolute error is O(ω−2n− 2 ) for N = 2n − 1 as shown in Figs. 5–6, while the decay rate is O(ω−N − 2 ) 3 for N < 2n − 1 in Fig. 7, which can confirm our error analysis in Theorem 3.2. Actually, the decay rate can be O(ω−2n− 2 ) as long as N ⩾ 2n − 1. Moreover, Fig. 8 in Example 7 shows that the modified complex integration method can achieve machine precision for approximating the related integrals with ω = 500 by only using four nodes. Also, we can see that for fixed n and N, the higher ω, the higher precision in Figs. 5–7. Remark 1. We devise these two methods under the condition that α + β > −1, which can make the convergence of the integrals (1.5) and (1.6) guaranteed well. But for some other special cases that α, β satisfy α + β ⩽ −1, our presented methods may be also efficient. For instance, we can see from Fig. 9 that the modified complex integration method for ∫1 6 ∫1 1 1 x 37 approximating the integrals 0 x− 5 arctan J (ω/x 7 )dx(α + β = − 35 < −1) and 0 x−2 arctan x3 J3 (ω/x 2 )dx(α + β = − 23 < 1+x 1 −1) can also reach high precision by only using two or three nodes. 4.3. Numerical comparison Here, we make a comparison for the complex integration method, the modified complex integration method and the Filon-type method proposed in [34] by numerical examples.

∫1

1

1

Example 8. We consider the approximation of the integral 0 x− 6 sin xJ2 (ω/x 6 )dx. The comparison between the two methods proposed in this paper is shown in Fig. 10. Tables 1–3 illustrate absolute errors and computation times of the three different methods. In Example 8, the same integral is calculated by the three methods. As shown in Fig. 10 and Tables 1–3, we can observe that the two proposed methods can achieve the same accuracy. Nevertheless, compared with the modified Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

16

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 6. The absolute error (left figure) and absolute error scaled by ω9.5 (right figure) of ˜ Q4MC [G], for the integral n = 4 and N = 2n − 1.

∫1 0

1

1

x− 2 sin x ln(x)J6 (ω/x 5 )dx with

Table 1 ∫1 1 1 Absolute errors and times for computing the integral 0 x− 6 sin xJ2 (ω/x 6 )dx by QnC [G] (method 1), QnMC ,2n−1 [G] (method 2) for n = 2, and Filon-type method with the nodes {−1, 1} and multiplicities {0, 0} for ω = 20, 40, 80, 160, 320, 640. Method 1

ω

Absolute error

20 40 80 160 320 640

4.62 2.33 3.60 1.34 2.92 6.19

× × × × × ×

10−4 10−5 10−7 10−8 10−10 10−12

Method 2 Times

Absolute error

2.09 2.02 1.90 1.85 1.95 1.84

4.62 2.33 3.60 1.34 2.92 6.19

s s s s s s

× × × × × ×

10−4 10−5 10−7 10−8 10−10 10−12

Filon-type method Times

Absolute error

1.57 1.54 1.53 1.62 1.61 1.62

6.18 3.81 1.74 1.77 1.98 1.23

s s s s s s

× × × × × ×

10−3 10−4 10−4 10−5 10−6 10−7

Times 1.96 1.97 1.93 2.01 2.15 2.03

s s s s s s

Table 2 ∫1 1 1 Absolute errors and times for computing the integral 0 x− 6 sin xJ2 (ω/x 6 )dx by QnC [G] (method 1), QnMC ,2n−1 [G] (method 2) for n = 3, and Filon-type method with the nodes {−1, 1} and multiplicities {0, 1} for ω = 20, 40, 80, 160, 320, 640. Method 1

ω

Absolute error

20 40 80 160 320 640

4.27 4.46 1.68 1.94 1.05 5.83

× × × × × ×

10−4 10−6 10−8 10−10 10−12 10−15

Method 2 Times

Absolute error

1.87 1.84 1.88 1.89 1.84 1.85

4.27 4.46 1.68 1.94 1.05 5.83

s s s s s s

× × × × × ×

10−4 10−6 10−8 10−10 10−12 10−15

Filon-type method Times

Absolute error

1.55 1.55 1.55 1.58 1.65 1.60

1.12 5.84 1.05 1.64 4.70 1.34

s s s s s s

× × × × × ×

10−3 10−5 10−6 10−6 10−7 10−8

Times 2.20 2.11 2.12 2.23 2.40 2.27

s s s s s s

complex integration method, the complex integration method takes much more computation time due to its complicated procedure. Moreover, the accuracy level of the Filon-type method is far lower than that of the proposed two methods Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 7. The absolute error (left figure) and absolute error scaled by ω4.5 (right figure) of Q4MC [G], for the integral and N = 2.

Fig. 8. The absolute error of the modified complex integration method for the integral

∫1 0

5 6

∫1 0

1

x− 3

∫1

sin x2 ex

0

17

1

1

x− 5 sin xJ3 (ω/x 2 )dx with n = 4

1

J 1 (500/x 4 )dx (left figure) and 5

1 2

x− ln(x) arctan x3 J2 (500/x )dx (right figure) with n from 2 to 20 and N = 2n − 1.

under the same conditions. Meanwhile, the Filon-type method takes nearly most computation time but obtains the lowest accuracy level. However, as mentioned above, the Filon-type method is valid for all smooth f on [0, 1], which relaxes the Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

18

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 9. The absolute error of the modified complex integration method for the integral and

∫1 0

∫1 0

6

x− 5

arctan x J ( 1+x 1

1

ω/x 7 )dx with n = 2, N = 2n − 1 (left figure)

1 2

x−2 arctan x3 J3 (ω/x )dx with n = 3, N = 2n − 1 (right figure).

Fig. 10. The comparison of the absolute errors for the complex integration method (left) and modified complex integration method (right). Here, n = 4 and N = 2n − 1 in the modified complex integration method. Table 3 ∫1 1 1 Absolute errors and times for computing the integral 0 x− 6 sin xJ2 (ω/x 6 )dx by QnC [G] (method 1),

QnMC ,2n−1 [G] (method 2) for n = 4, and Filon-type method with the nodes {−1, 0, 1} and multiplicities {0, 0, 1} for ω = 20, 40, 80, 160, 320, 640. Here, ‘MP’ indicates that the produced numerical results have reached machine precision in Matlab. Method 1

ω 20 40 80 160 320 640

Method 2

Absolute error 1.53 8.56 5.67 1.36 1.19 MP

× × × × ×

−4

10 10−7 10−10 10−12 10−15

Times 1.88 1.89 1.88 1.95 1.90 2.07

s s s s s s

Filon-type method

Absolute error 1.53 8.56 5.67 1.36 1.19 MP

× × × × ×

−4

10 10−7 10−10 10−12 10−15

Times 1.60 1.63 1.61 1.58 1.59 1.58

s s s s s s

Absolute error 1.32 8.40 2.13 1.93 3.62 1.56

× × × × × ×

−3

10 10−5 10−6 10−6 10−7 10−8

Times 2.41 2.34 2.41 2.48 2.85 2.62

s s s s s s

strict requirement that f is analytic in G1 . Therefore, these methods have their respective advantages and disadvantages, and can be applied to different situations. It should be also mentioned that the complex integration method holds only for the Bessel functions of the natural number order m (see Figs. 2–4). Fortunately, the modified complex integration method is very efficient for an arbitrary Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

19

non-negative real number m (see Figs. 5–8). Therefore the modified complex integration method shares big advantages due to its simplicity, efficiency and a wide range of applications. In fact, they have the same error order when N ⩾ 2n − 1. 5. Conclusion In this paper, we have proposed two approaches to calculate the highly oscillatory integrals of the forms (1.1) and (1.2). We first transform the integrals (1.1) and (1.2) into (1.5) and (1.6) by change of variable t = x−β , β > 0. Then based on the well-known Cauchy’s residue theorem, the integrals (1.5) and (1.6) can be transformed into the infinite integrals on [0, +∞), where the integrands are not oscillatory and decay exponentially fast. The resulting infinite integrals can be efficiently computed by using the Gauss–Laguerre quadrature rule. Moreover, we provide error analysis in inverse powers of ω in Theorems 2.2 and 3.2. Furthermore, the accuracy can be improved greatly by adding more nodes for fixed ω due to 3 its error order satisfying O(ω−2n− 2 ). Meanwhile, the two methods can produce higher accuracy as ω increases for fixed n. Additionally, numerical examples are shown to confirm our theoretical analysis. We also observe from numerical examples that the two methods proposed in this paper are more accurate than the Filon-type method at the same computational cost. Acknowledgments The authors would like to express their most sincere thanks to the referees and editors for their very helpful comments and suggestions, which greatly improved the quality of this paper. References [1] G. Arfken, Mathematical Methods for Physicists, third ed., Academic Press, Orlando, FL, 1985. [2] G. Bao, W. Sun, A fast algorithm for the electromagnetic scattering from a large cavity, SIAM J. Sci. Comput. 27 (2005) 553–574. [3] H. Brunner, Open problems in the computational solution of Volterra functional equations with highly oscillatory kernels, in: HOP 2007: Effective Computational Methods for Highly Oscillatory Solutions, Isaac Newton Institute. [4] H. Brunner, On the numerical solution of first-kind Volterra integral equations with highly oscillatory kernels, in: Highly Oscillatory Problems: From Theory to Applications, HOP 13–17, September, 2010, Isaac Newton Institute. [5] P.J. Davies, D.B. Duncan, Stability and convergence of collocation schemes for retarded potential integral equations, SIAM J. Numer. Anal. 42 (2004) 1167–1188. [6] D. Huybrechs, S. Vandewalle, A sparse discretization for integral equation formulations of high frequency scattering problems, SIAM J. Sci. Comput. 29 (2007) 2305–2328. [7] R. Chen, Numerical approximations to integrals with a highly oscillatory Bessel kernel, Appl. Numer. Math. 62 (2012) 636–648. [8] R. Chen, On the evaluation of Bessel transformations with the oscillators via asymptotic series of Whittaker functions, J. Comput. Appl. Math. 250 (2013) 107–121. [9] R. Chen, Numerical approximations for highly oscillatory Bessel transforms and applications, J. Math. Anal. Appl. 421 (2015) 1635–1650. [10] G.A. Evans, K.C. Chung, Some theoretical aspects of generalised quadrature methods, J. Complexity 19 (2003) 272–285. [11] D. Levin, Procedures for computing one-and-two dimensional integrals of functions with rapid irregular oscillations, Math. Comp. 38 (1982) 531–538. [12] S. Olver, Numerical approximation of vector-valued highly oscillatory integrals, BIT 47 (2007) 637–655. [13] R. Piessens, M. Branders, Modified Clenshaw–Curtis method for the computation of Bessel function integrals, BIT 23 (1983) 370–381. [14] S. Xiang, Numerical analysis of a fast integration method for highly oscillatory functions, BIT 47 (2007) 469–482. [15] S. Xiang, H. Wang, Fast integration of highly oscillatory integrals with exotic oscillators, Math. Comp. 79 (2010) 829–844. [16] S. Xiang, Y. Cho, H. Wang, H. Brunner, Clenshaw–Curtis–Filon-type methods for highly oscillatory Bessel transforms and applications, IMA J. Numer. Anal. 31 (4) (2011) 1281–1314. [17] H. Wang, A unified framework for asymptotic analysis and computation of finite Hankel transform, J. Math. Anal. Appl. 483 (2019) 123640. [18] R. Chen, C. An, On the evaluation of infinite integrals involving Bessel functions, Appl. Math. Comput. 235 (2014) 212–220. [19] Z. Xu, S. Xiang, Computation of Highly Oscillatory Bessel Transforms with Algebraic Singularities, Technical Report, 27-03-2014, Central South University, 2014. [20] S. Lewanowicz, Evaluation of Bessel function integrals with algebraic singularities, J. Comput. Appl. Math. 37 (1991) 101–112. [21] H. Wang, L. Zhang, D. Huybrechs, Asymptotic expansions and fast computation of oscillatory Hilbert transforms, Numer. Math. 123 (2013) 709–743. [22] Z. Xu, S. Xiang, G. He, Efficient evaluation of oscillatory Bessel Hilbert transforms, J. Comput. Appl. Math. 258 (2014) 57–66. [23] Z. Xu, S. Xiang, Numerical evaluation of a class of highly oscillatory integrals involving Airy functions, Appl. Math. Comput. 246 (2014) 54–63. [24] Z. Xu, S. Xiang, On the evaluation of highly oscillatory finite Hankel transformusing special functions, Numer. Algorithms 72 (2016) 37–56. [25] H. Kang, C. Ling, Computation of integrals with oscillatory singular factors of algebraic and logarithmic type, J. Comput. Appl. Math. 285 (2015) 72–85. [26] G.V. Milovanović, Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures, Comput. Math. Appl. 36 (1998) 19–39. [27] D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (2006) 1026–1048. [28] H. Mo, S. Xiang, On the asymptotic order of Filon-type methods for highly oscillatory integrals with an algebraic singularity, Appl. Math. Comput. 217 (22) (2011) 9105–9110. [29] W. Gautschi, Computing polynomials orthogonal with respect to densely oscillating and exponentially decaying weight functions and related integrals, J. Comput. Appl. Math. 184 (2005) 493–504. [30] A.I. Hascelik, On numerical computation of integrals with integrands of the form f (x) sin(ω/xr) on [0, 1], J. Comput. Appl. Math. 223 (2009) 399–408. [31] A.I. Hascelik, Suitable Gauss and Filon-type methods for oscillatory integrals with an algebraic singularity, Appl. Numer. Math. 59 (2009) 101–118.

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.

20

H. Wang and H. Kang / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

[32] H. Kang, S. Xiang, On the calculation of highly ocillatory integrals with an algebraic singularity, Appl. Math. Comput. 217 (8) (2010) 3890–3897. [33] H. Kang, S. Xiang, Efficient quadrature of highly oscillatory integrals with algebraic singularities, J. Comput. Appl. Math. 237 (2013) 576–588. [34] H. Kang, J. Ma, Quadrature rules and asymptotic expansions for two classes of oscillatory Bessel integrals with singularities of algebraic or logarithmic type, Appl. Numer. Math. 118 (2017) 277–291. [35] G.N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge University Press, 1952. [36] E.T. Whittaker, G.N. Watson, A Course in Modern Analysis, fourth ed., Cambridge University Press, 1990. [37] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. [38] M.J. Ablowitz, A.S. Fokas, Complex Variables: Introduction and Applications, Cambridge University Press, 1997. [39] P. Henrici, Applied and Computational Complex Analysis, Vol. I, Wiley and Sons, New York, 1974. [40] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second ed., Academic Press, 1984.

Please cite this article as: H. Wang and H. Kang, Numerical methods for two classes of singularly oscillatory Bessel transforms and their error analysis, Journal of Computational and Applied Mathematics (2019) 112604, https://doi.org/10.1016/j.cam.2019.112604.