Numerical approximations for highly oscillatory Bessel transforms and applications

J. Math. Anal. Appl. 421 (2015) 1635–1650

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa

Numerical approximations for highly oscillatory Bessel transforms and applications ✩ Ruyun Chen College of Science, Guangdong Ocean University, Zhanjiang, Guangdong 524088, China

a r t i c l e

i n f o

Article history: Received 21 April 2012 Available online 19 August 2014 Submitted by K. Driver Keywords: Bessel transforms Highly oscillatory integral equation Complex integration method Steepest descent method Numerical solution

a b s t r a c t This paper presents an efficient numerical method for approximating highly oscillatory Bessel transforms. Based on analytic continuation, we transform the integrals into the problems of integrating the forms on [0, +∞) with the integrand that does not oscillate and decays exponentially fast, which can be efficiently computed by using Gauss–Laguerre quadrature rule. We then derive the error of the method depending on the frequency and the node number. Moreover, we apply the scheme for studying the approximations of the solutions of two kinds of highly oscillatory integral equations. Preliminary numerical results show the efficiency and accuracy of numerical approximations. © 2014 Elsevier Inc. All rights reserved.

1. Introduction The problem of evaluating the highly oscillatory Bessel transforms a I[f ] =

f (x)Jv (ωx)dx,

ω1

(1.1)

0

occurs in many areas of science and technology, for example, in astronomy, optics, electro-magnetics, seismology image processing and applied mathematics (see [3,4,10,19,31,32]). Here, f (x) is a smooth real-valued function, v is an arbitrary nonnegative integer and Jv (ωx) is the Bessel function of the first kind and of order v. In most of the cases, the integrals cannot be done analytically and one has to rely on numerical methods. However, for large ω, the integrands become highly oscillatory and present serious difficulties in obtaining numerical convergence of the integrations. For example, classical quadrature rules of the Newton– Cotes type or Gaussian type are based on polynomial interpolation. It is well known that polynomials are ✩

The work is supported partly by NSF of Guangdong Ocean University (No. C13443). E-mail address: [email protected].

http://dx.doi.org/10.1016/j.jmaa.2014.08.021 0022-247X/© 2014 Elsevier Inc. All rights reserved.

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1636

not suited for the approximation of oscillatory functions, and the integration error of these quadrature rules increases rapidly with increasing ω. The theoretical and numerical aspects of the integrals have received considerable attention in the last few years and many efficient methods have been developed. Modified Clenshaw–Curtis methods [29] are 1 applicable to the integral 0 f (x)Jv (ωx)dx for non-negative order v, which is presented by replacing the smooth factor of the integrand by a truncated Chebyshev series approximation. Levin-type methods (see [22, b 23,28,38,41]) are developed for evaluating a f (x)Jv (ωx)dx under the assumption that the weight functions satisfies certain differential conditions for 0 ∈ / [a, b]. Based on Lagrange’s identity    zLS − SMz = Z S(ω, x), z(x) x

(1.2)

by choosing the differential operator L to satisfy LS = 0 for the weight function S(ω, x), where M is the adjoint operator of L, generalized quadrature rules are presented under the condition that 0 ∈ / [a, b]. This method is efficient for more general irregular oscillatory weight functions S(ω, x). But a disadvantage is that the approach only works for 0 ∈ / [a, b] [9,13,40]. Based on He [16,17] and Watson’s work [35,36], Chen and Liang proposed a homotopy perturbation method to consider the asymptotic expression of the Bessel, Anger and Weber transformations [7,8]. All these methods share an advantageous property that the accuracy greatly improves when ω increases. Complex integration method [27,37] and steepest descent method [18] for the form 1 f (x)eiωx dx

(1.3)

−1

have been developed. If f is analytic in a sufficiently large complex region G that contains [−1, 1], the integrals can be transformed into the problems of integrating two integrals on [0, +∞), where the integrands are not oscillatory and decay exponentially fast [18,21,27]. In the present work, an entirely different method is presented to compute the integrals (1.1). We consider the numerical computation of the integral following Wong [37], Milovanović [27] and Huybrechs and Vandewalle’s work [18]. The numerical solution of the integral equation of the first kind with a highly oscillatory Bessel kernel x

  J0 ω(x − t) y(t)dt = f (x),

x ∈ [0, T ],

(1.4)

0

has attracted much attention during the last few years [6,10,30,33], where y(x) is the unknown function whose value is to be determined in the interval [0, T ], f (x) is a given sufficiently smooth function and ω is a large parameter. One feature of the Volterra integral equation (1.4) is of particular note: when ω  1, the kernel function J0 (ω(x − t)) would become highly oscillatory. This means that some general numerical methods may not be immediately applicable to this equation. When f (x) ∈ C 1 [0, T ], the unique solution of Eq. (1.4) can be given by [30, Eq. 1.8.1] as

y(x) =

1 ω



d2 + ω2 dx2

2 x

  (x − t)J1 ω(x − t) f (t)dt,

x ∈ [0, T ].

(1.5)

0

At the same time, we will consider the numerical solution of the Fredholm integral equation with a highly oscillatory Bessel kernel

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1637

a y(x) −

Jv (ωt)y(t)dt = f (x),

(1.6)

0

where the exact solutions may be denoted by the form [30, Eq. 4.8.2]  1−

y(x) = f (x) +

−1 a

a Jv (ωt)dt 0

f (t)Jv (ωt)dt.

(1.7)

0

Although the exact solutions of the two kinds of equations are available, they involve a highly oscillatory Bessel integral, which cannot be integrated in closed form. Therefore, one needs to consider the use of reliable numerical schemes to provide the approximations. For the form x K(x, t)y(t)dt = f (x),

x ∈ [a, b]

(1.8)

a

where K : [a, b] × [a, b] → R, some numerical approaches have been constructed by replacing the integral by a numerical quadrature formula (see [5,6,12,15,24,26]). However, these approaches cannot applied to (1.4) and (1.6) since the kernel Jv (ωx) is highly oscillatory for large values of ω. The computation of the integral containing Jv (ωx) by standard quadrature methods is exceedingly difficult and the cost steeply increases with ω [20,22]. Other aim of the article is to apply the results of evaluating the integrals (1.1) to solve these equations. An outline of this paper is organized as follows. In the next section, we present the details of the proposed method for computing the integrals involving Bessel function Jv (ωt) and give the error analysis of the proposed method. In Section 3, we first transform the exact solution (1.5) into an equivalent form containing the kernel J0 (ωt). We then apply the results of Section 2 for the equivalent form and the exact solution (1.7). At the same time, numerical examples are provided to demonstrate the effectiveness of the proposed methods. 2. Numerical evaluations of the integrals involving Bessel function We rewrite (1.1) as a I[f ] =

F (x)Jv (ωx)dx +

j=0

0

2n+1 where F (x) = f (x) − j=0 p. 707] and [2, Eq. 11.1.1]) a xj Jv (ωx)dx = 0

2n+1

f (j) (0) j j! x

and

) 2j Γ ( j+v+1 2 ω j+1 Γ ( j−v+1 ) 2

+

a 0

f (j) (0) j!

a xj Jv (ωx)dx = I1 [f ] + I2 [f ],

(2.1)

0

xj Jv (ωx)dx can be represented as follows ([25, p. 44], [14,

a a(j + v − 1) (2) (2) Jv (ωa)sj−1,v−1 (ωa) − j Jv−1 (ωa)sj,v (ωa), ωj ω

(2.2)

(2)

where Γ is the gamma function and sj,v (z) denotes the Lommel function of the second kind [25,34]. The (2)

computation of sj,v (z) admits the following asymptotic expansion [25, pp. 351–352]: (2) sj,v (z)

=z

j−1



2 2 2 2   p {(j − 1) − v } · · · {(j − 2p + 1) − v } 1 + · · · + (−1) + O z j−2p−2 , 2p z

(2.3)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1638

where    j−2p−2  22p−1  Γ (p + 1 − j−v 2 )Γ (p + 1 − ≤  O z  1−j+v 2p−j z Γ ( 2 )Γ ( 1−j−v ) 2

j+v  2 )

see [39].



(2.4)

(2)

Therefore, sj,v (z) can be efficiently approximated by a few terms of

2 2 2 2 (j − 1)2 − v 2 p {(j − 1) − v } · · · {(j − 2p + 1) − v } , z j−1 1 − + · · · + (−1) z2 z 2p

(2.5)

when z is large values. Remark 1. Since the Whittaker W function Wμ,ν (x) has the following formula [2, p. 507]:  Wμ,ν (x) =



 1 μ 1 − Wμ,ν (x) − Wμ−1,ν (x), 2 x x

(2.6)

(2n) 1 , k = 0, 1, · · · , 2n. Therefore, it is true that the expression xk 2n+1 f (j) (0) j f (x)− j=0 x j! [ W− 2v−k , 2v−k (±2ωxi)](2n) is a sufficiently smooth function. This expression will be used x1/2 2 2

then Wμ,ν (x) is the combination of in later sections.

From [34], Bessel function Jv (ωx) may be written as the form ( ωx )v Jv (ωx) = √ 2 πΓ (v + 12 )

1

 v− 12 iωxτ e dτ. 1 − τ2

(2.7)

−1

By substituting (2.7) into (2.1), we obtain ( ω )v I[f ] = √ 2 πΓ (v + 12 )

 1

a v

x F (x)



1−τ

 1 2 v− 2 iωxτ e

 dτ dx.

(2.8)

−1

0

1

In the integrals (2.8), since (1 −τ 2 )v− 2 eiωxτ is analytic in the half-strip of the complex plane: −1 ≤ (τ ) ≤ 1 and (τ ) ≥ 0, by the Cauchy Residue Theorem [1], then we have 



1 − τ2

v− 12

eiωxτ dτ = 0,

(2.9)

Γ 1+Γ 2+Γ 3+Γ 4

with all contours taken in the counterclockwise direction which is depicted in Fig. 1. Noting that 

 v− 12 iωxτ 1 − τ2 e dτ =

Γ2





1 − τ2

v− 12

eiωxτ dτ,

(2.10)

−1

Γ1



1

1−τ

 1 2 v− 2 iωxτ e

iωx

dτ = ie

R 

p2 − 2ip

v− 12

e−ωxp dp,

(2.11)

0

      1 2 v− 2 iωxτ   → 0 as R → +∞, 1 − τ e dτ   Γ3

(2.12)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1639

Fig. 1. Illustration of the integration path for the kernel eiωxτ .

and 



1−τ

 1 2 v− 2 iωxτ e

−iωx

dτ = −ie

R 

p2 + 2ip

v− 12

e−ωxp dp,

(2.13)

0

Γ4

we get 

 = − lim

R→+∞ Γ2

Γ1

 − lim

R→+∞ Γ4

,

(2.14)

that is, 1



1−τ

 1 2 v− 2 iωxτ e

+∞  v− 12 v− 12  −ωxp  −iωx  2  dτ = − ieiωx p2 − 2ip dp. ie p + 2ip e

−1

(2.15)

0

By letting q = ωxp as x = 0, we have 1



1−τ

 1 2 v− 2 iωxτ e

+∞ 

dτ =

−1

ie−iωx ωx



q ωx

2

2iq + ωx

v− 12

ieiωx − ωx



q ωx

2

2iq − ωx

v− 12  e−q dq.

(2.16)

0

Replacing (2.16) into (2.8) yields

1 ik+1 v! I[f ] = √ 1 v−k k!(v − k)!ω v−k πΓ (v + 2 ) k=0 2 v

+∞ 



− (−1)k eiωx 0

a 0

 +∞  F (x) −iωx q 2v−k  e−q dq e v−k 2 x q + 2iqωx 0



q 2v−k q 2 − 2iqωx

e−q dq dx.

(2.17)

Let 1 M1 (x) = √ π

∞  0

q 2v−k q 2 + 2iωxq

e

−q

dq,

1 M2 (x) = √ π

By means of the definition of Whittaker W function [34]

∞ 0

q 2v−k  e−q dq. q 2 − 2iωxq

(2.18)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1640



 1−2u+2v 1−2u−2v ( ) ( ) n n 2 2 1+ , (−1)n n n!x n=1

Wu,v (x) = xu e−x/2

+∞

x → ∞,

(2.19)

we have 2v−k

2v−k−1 (2l − 1) (2iωx) 2 eωxi W−u,u (2ωxi) 22v−k

l=1

M1 (x) =

(2.20)

and 2v−k l=1

M2 (x) = where u =

2v−k 2 .

(2l − 1)

22v−k

(−2iωx)

2v−k−1 2

e−ωxi W−u,u (−2ωxi),

(2.21)

Hence, there exists

 a 2v−k k+1 v

k−1 iv+ 2 v! l=1 (2l − 1) 1 I[f ] = x 2 F (x)eωxi W−u,u (2ωxi)e−iωx dx 3k−1 k−1 1 Γ (v + 2 ) k=0 k!(v − k)!2− 2 ω − 2 0

v+ k−1 2

a

− (−1)

x

k−1 2

F (x)e

−ωxi

 iωx

W−u,u (−2ωxi)e

dx .

(2.22)

0

By using the techniques of (2.9), the following equations are true: a 0

i 1−k W−u,u (2ωxi)dx = ω 2 x F (x)

∞

F (a − (a −

0

∞

i − ω

0

ip ω)

ip 1−k 2 ω)

 e W−u,u

F (− ip ω) (− ip ω)

p

1−k 2

  ip 2ω a − i e−p dp ω

 p

e W−u,u

  ip 2ω − i e−p dp ω

(2.23)

and a 0

i 1−k W−u,u (−2ωxi)dx = ω x 2 F (x)

−

∞ 0

i ω

   ip p −2ω i e−p dp e W −u,u ip 1−k ω ( ) 2 F ( ip ω) ω

∞ 0

   ip p −2ω a + i e−p dp. e W −u,u ip 1−k ω 2 (a + ω ) F (a +

ip ω)

(2.24)

Applying a Gauss–Laguerre quadrature rule with n points xk and weights wk , we have  n 2v−k     k+3 v

iv+ 2 v! l=1 (2l − 1)

ixs ixs 1 xs − G1 − Qn [F ] = ws e G 1 a − k−3 ω ω Γ (v + 12 ) k=0 k!(v − k)!2− 3k−1 2 ω− 2 s=1      2n+1 n

f (j) (0) ixs ixs v+ k−1 x s 2 − G2 a + + ws e G 2 − (−1) ω ω j! s=1 j=0 

× where G1 (x) = x

k−1 2

) 2j Γ ( j+v+1 2 ω j+1 Γ ( j−v+1 ) 2

+

  a (2) (2) (j + v − 1)J (ωa)s (ωa) − J (ωa)s (ωa) , v v−1 j−1,v−1 j,v ωj

F (x)W− 2v−k , 2v−k (2ωxi) and G2 (x) = x 2

2

k−1 2

F (x)W− 2v−k , 2v−k (−2ωxi). 2

2

(2.25)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1641

Table 1 Truncated errors for W− 2v−k , 2v−k (2ωxi) with N = 25, v = 5 and x = 1. 2

2

ω

k=0

k=1

k=2

k=3

k=4

k=5

50 100 150 200 500

5.64×−29 2.75×−38 9.64×−44 1.29×−47 6.00×−60

1.51×−28 1.04×−37 4.47×−43 6.93×−47 5.07×−59

3.74×−28 3.64×−37 1.90×−42 3.42×−46 3.95×−58

8.40×−28 1.15×−36 7.40×−42 1.53×−46 2.80×−57

1.68×−27 3.26×−36 2.56×−41 6.12×−45 1.77×−56

2.95×−27 8.05×−36 7.76×−41 2.14×−44 9.78×−56

Remark 2. The efficient implementations of Wμ,κ (x) are based on the asymptotic formula (2.19). W− 2v−k , 2v−k (2ωxi) and W− 2v−k , 2v−k (−2ωxi) can be efficiently approximated by truncating (2.19) with 2 2 2 2 a few terms once ω is large. Table 1 shows the efficiency of the approximation. To show the high accuracy of the approximation we use Maple 14 with 500-digit to compute the approximations. Throughout the numerical results in this paper, we evaluate W·,· (·) by truncating after the first 25 terms of (2.19). Theorem 2.1. Suppose that f is an analytic function in the complex plane {0 ≤ (z) ≤ a}. If f (k) (x) are a bounded for k = 0, 1, 2, · · · , 2n, then the error of (2.25) for approximating 0 f (x)Jv (ωx)dx is given by the formula  En [F ] = I[F ] − Qn [F ] = O

 1 (n!)2 , (2n)! ω 2n+3/2

ω  1.

(2.26)

Proof. If f (2n) (x) is a bounded value for x ∈ (0, ∞), then the error formula for the n-point Gauss–Laguerre  +∞ quadrature rule applied to the integrals 0 f (x)e−x dx can be given by the formula in [11] (n!)2 (2n) f (ξ), (2n)!

0 < ξ < +∞.

(2.27)

Due to 2v−k k+3 v 1 (n!)2 iv+ 2 v! l=1 (2l − 1) En [f ] = k−3 Γ (v + 12 ) (2n)! k=0 k!(v − k)!2− 3k−1 2 ω− 2    k−1    (2n)   iξ1 iξ1 2 iξ1 ξ1 e W− 2v−k , 2v−k 2ω a − i × a− F a− 2 2 ω ω ω    k−1    (2n)  iξ2 iξ2 2 iξ2 ξ2 e W− 2v−k , 2v−k 2ω − i − − F − 2 2 ω ω ω   k−1    (2n)   k−1 iξ3 2 iξ3 iξ3 ξ3 e W− 2v−k , 2v−k −2ω i F − (−1)v+ 2 2 2 ω ω ω  k−1      (2n)   iξ4 iξ4 2 iξ4 ξ4 e W− 2v−k , 2v−k −2ω a + i , − a+ F a+ 2 2 ω ω ω where ξ1 , ξ2 , ξ3 , ξ4 ∈ (0, +∞). Let L(x) =

 (2n) L(x) =

e−aωi ω 2n (2ωi)

2v−k 2

F (a−

d(2n)

iξ1 ω

)eξ1 W− 2v−k , 2v−k (2ω(a− 2 iξ (a− ω1

2 1−k ) 2

iξ1 ω

)i)

, then

)n ( 1 )n   F (a−ip){1++∞ (−1)n ( 1+4v−2k 2 2 } n=1

(a−ip) dp2n

n!(2ω(a−ip)i)n 1−k 2

    

. p=ξ1 /ω

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1642

 Fig. 2. Numerical quadrature of 01 xex J0 (ωx)dx by the method (2.25): The absolute error scaled by ω 5.5 of Q2 [f ] (left figure) and 7.5 scaled by ω of Q3 [f ] (right figure).

 5.5 1 Fig. 3. Numerical quadrature of 01 1+x of Q2 [f ] (left figure) and 2 J2 (ωx)dx by the method (2.25): The absolute error scaled by ω 7.5 scaled by ω of Q3 [f ] (right figure).

By Leibniz’s Theorem, [1 + that F (a −

+∞

1+4v−2k )n ( 12 )n (k) n( 2 n=1 (−1) n!(2ω(a−ip)i)n ]

iξ1 ξ1 (2ω(a , 2v−k ω )e W− 2v−k 2 2 iξ1 1−k (a − ω ) 2

Similarly, other three derivatives satisfy O(

En [f ] =

−

1 ω 2n (2ωi)

2v−k 2

are bounded for k = 0, 1, · · · , 2n. Then it is true

2n iξ1 ω )i)

 =O



1 ω 2n (2ωi)

2v−k 2

.

) too. Therefore,

2v−k     v v! l=1 (2l − 1) (n!)2 1 1 1 (n!)2

= O . O k−3 2v−k (2n)! ω 2n+ 32 Γ (v + 12 ) (2n)! k=0 k!(v − k)!2− k−1 2 ω− 2 ω 2n (2ω) 2

This proves the results (2.26). 2 Example 1. Let us consider the quadrature method (2.25) for Example 2. Let us consider the quadrature method (2.25) for

1 0

1

xex J0 (ωx)dx (see Fig. 2).

1 J (ωx)dx 0 1+x2 2

(see Fig. 3).

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1643

3. Applications to the approximations of the solutions of the equations In this section, we consider the applications to the following integral equations: x

  J0 ω(x − t) y(t)dt = f (x),

x ∈ [0, T ]

(3.1)

0

where the exact solutions may be denoted by the form ([30, Eq. 1.8.1] and [33]) 1 y(x) = ω



d2 + ω2 dx2

2 x

  (x − t)J1 ω(x − t) f (t)dt,

x ∈ [0, T ]

(3.2)

0

and a y(x) −

Jv (ωt)y(t)dt = f (x),

(3.3)

0

where the exact solutions may be denoted by the form [30, Eq. 4.8.2]  1−

y(x) = f (x) +

−1 a

a Jv (ωt)dt 0

f (t)Jv (ωt)dt.

(3.4)

0

3.1. The approximations of the solutions of Eq. (3.1) Lemma 3.1. (See [33].) The exact solution of (3.1) can be transformed into the following equivalent form: x



y(x) = f (x) + ω 0

J1 (ω(x − t)) f (t)dt. x−t

(3.5)

Although an explicit integral expression of the exact solution is derived, the solution is still not available because the formula (3.5) involves a highly oscillatory Bessel integral which has a removable singularity x at the point t = x. From [11], the integral 0 J1 (ω(x−t)) f (t)dt has one practical difficulty, that is, it is x−t oscillatory. Wang and Xiang obtained [33]



y(x) = f (x) + ω

x 2

x f1 (t)J0 (ωt)dt + ω

0

f2 (t)J1 (ωt)dt,

(3.6)

0

where f1 (t) = f (x − t),

f2 (t) = −f  (x − t).

(3.7)

Thus, the numerical solution of (3.1) is converted into the computation of two highly oscillatory Bessel integrals. We rewrite (3.6) as

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1644

x 



y(x) = f (x) + ω

f1 (t) −

2

2n+1

k=0

0

(k) f1 (0) k t J0 (ωt)dt k!

x x  2n+1 2n+1

f (k) (0)

f (k) (0)  k 2 2 1 t J1 (ωt)dt + ω +ω f2 (t) − tk J0 (ωt)dt k! k! k=0

0 2n+1

+ω

k=0

(k)

f2 (0) k!

k=0

0

x tk J1 (ωt)dt.

(3.8)

0

Noticing that −

1 dJ0 (ωx) = J1 (ωx), ω dx

(3.9)

then  

y(x) = f (x) −

f2 (x) −

2n+1

k=0

x  f2 (t) −

+

k=0

0

+

2n+1

k=0

2n+1

(k)

f2 (0) (k − 1)!

x  2n+1 (k)

f (k) (0) f2 (0) k 2 k 1 x J0 (ωx) + ω t J0 (ωt)dt f1 (t) − k! k!

(k)

f2 (0) k t k!

k=0

0

 J0 (ωt)dt + ω

2

2n+1

k=0

(k)

f1 (0) k!

x tk J0 (ωt)dt 0

x tk−1 J0 (ωt)dt,

(3.10)

0

x

where the moments I[k, ω, x] = 0 tk J0 (ωt)dt can be evaluated by (2.2). From (2.25) and (3.10), we have the following formula for evaluating the solutions (3.2) yn (x) =

n F1 (− ixωs )exs W0,0 (2xs ) F1 (x − ω

w − √ s 1 1 2 2 π s=1 xs2 F1 ( ixωs )exs W0,0 (−2ω( ixωs )i)

F1 (x +

ixs xs ω )e W0,0 (2ω(x 1 (iω(x − ixωs )) 2

ixs xs ω )e W0,0 (−2ω(x 1 (−iω(x + ixωs )) 2

− 1 xs2  ixs xs n

F (− ω )e W0,0 (2xs ) F2 (x − 1 ws 2 − + 1√ 1 2 2 πω s=1 xs2

+

+

F2 ( ixωs )exs W0,0 (−2ω( ixωs )i) 1

−

xs2 + ω2

2n+1

k=0

where F1 (x) = f1 (x) −

F2 (x +

+

−

ixs ω )i)

ixs ω )i)

ixs xs ω )e W0,0 (2ω(x 1 (iω(x − ixωs )) 2

ixs xs ω )e W0,0 (−2ω(x 1 (−iω(x + ixωs )) 2

+

−

ixs ω )i)

ixs  ω )i)

2n+1 (k)

f (k) (0) f1 (0) 2 I[k, ω, x] + I[k − 1, ω, x], k! (k − 1)!

f (x) − F2 (x)J0 (ωx)

(3.11)

k=0

2n+1 k=0

(k)

f1

(0) k k! x

and F2 (x) = f2 (x) −

2n+1 k=0

(k)

f2

(0) k k! x .

Theorem 3.1. The absolute error of the approximate solutions (3.11) for evaluating the exact solutions (3.2) of the Volterra integral equation of the first kind with a highly oscillatory Bessel kernel x 0

  J0 ω(x − t) y(t)dt = f (x),

x ∈ [0, T ],

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1645

Table 2 Relative errors of Qn [F ] with n = 2. Here f (x) = x sin(x). x

ω=5

ω = 10

ω = 20

ω = 40

ω = 80

1 2 3 4

7.25 × 10−5 8.55 × 10−5 1.18 × 10−3 5.28 × 10−5

1.05 × 10−6 5.89 × 10−6 5.36 × 10−5 8.18 × 10−6

7.68 × 10−8 4.13 × 10−7 8.27 × 10−7 2.19 × 10−7

5.56 × 10−9 1.25 × 10−8 1.07 × 10−8 1.35 × 10−8

1.67 × 10−10 6.66 × 10−10 2.25 × 10−9 6.63 × 10−10

Table 3 Relative errors of Qn [F ] with n = 3. Here f (x) = x sin(x). x

ω=5

ω = 10

ω = 20

ω = 40

ω = 80

1 2 3 4

5.93 × 10−7 6.52 × 10−7 6.13 × 10−6 1.41 × 10−7

1.48 × 10−8 4.34 × 10−10 7.87 × 10−8 9.23 × 10−9

4.38 × 10−11 1.84 × 10−10 4.65 × 10−10 5.88 × 10−11

2.16 × 10−12 1.25 × 10−12 2.46 × 10−12 9.88 × 10−13

1.34 × 10−14 1.99 × 10−14 4.70 × 10−14 1.20 × 10−14

satisfies   2   1 y(x) − yn (x) = O (n!) . (2n)! ω 2n−1/2 x Proof. From Theorem 2.1, we obtain all of the errors of the method (2.25) for approximating 0 (f1 (t) − x 2n+1 f2(k) (0) k  2n+1 f1(k) (0) k (n!)2 1 k=0 k=0 k! t )J0 (ωt)dt and 0 (f2 (t) − k! t ) J0 (ωt)dt satisfying O( (2n)! ω 2n+3/2 ). Thus, there exists     2   (n!)2 1 1 y(x) − yn (x) = ω 2 × O (n!) +O (2n)! ω 2n+3/2 (2n)! ω 2n+3/2       1 1 1 (n!)2 (n!)2 (n!)2 + O = O . =O (2n)! ω 2n−1/2 (2n)! ω 2n+3/2 (2n)! ω 2n−1/2 This completes the proof of Theorem 3.1.

2

Example 3. Let us consider the approximate solution yn (x) for x

  J0 ω(x − t) y(t)dt = x sin(x),

0

whose exact solution can be represented as x y(x) = x

      J0 (ωxt) x(1 − t) sin x(1 − t) ω 2 − 1 + 2 cos x(1 − t) dt

0

(see Tables 2–3 and Fig. 4). 3.2. The approximations of the solutions of Eq. (3.3) According to the following identities ([2, Eqs. 11.1.3–11.1.4] and [14, Eq. 6.511.6]): a

a J0 (t)dt − 2

J2u (t)dt = 0

0

u−1

k=0

J2k+1 (a),

(3.12)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1646

Fig. 4. The absolute error scaled by ω 3.5 of Q2 [f ] (left figure), and scaled by ω 5.5 of Q3 [f ] (right figure), and here we choose f (x) = x sin(x).

a J2u+1 (t)dt = 1 − J0 (a) − 2

u

J2k (a)

(3.13)

k=1

0

and aω  aωπ  J1 (aω)H0 (aω) − J0 (aω)H1 (aω) , J0 (t)dt = aωJ0 (aω) + 2

(3.14)

0

where Hj (z) denotes the StruveH functions [2, pp. 495–498], we have a

 2 u−1 aπ  J1 (aω)H0 (aω) − J0 (aω)H1 (aω) − J2u (ωt)dt = aJ0 (aω) + J2k+1 (aω) 2 ω

(3.15)

k=0

0

and a J2u+1 (ωt)dt =

u 1 1 2

− J0 (aω) − J2k (aω). ω ω ω k=1

0

Substituting (2.25) into (3.4) yields  yn (x) = f (x) + 1 − 

−1

a Jv (ωt)dt 0

2v−k k+3 v

iv+ 2 v! l=1 (2l − 1) 1 k−3 Γ (v + 12 ) k=0 k!(v − k)!2− 3k−1 2 ω− 2

 k−1      2 ixs ixs xs e W− 2v−k , 2v−k 2ω a − i × ws F a− 2 2 ω ω s=1   k−1       n

ixs 2 ixs ixs xs − e W− 2v−k , 2v−k 2ω − i ws − F − 2 2 ω ω ω s=1  n k−1   

 ixs  2  ixs  ixs v+ k−1 2 − (−1) exs W− 2v−k , 2v−k −2ω i ws F 2 2 ω ω ω s=1 n



ixs a− ω

(3.16)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1647

Table 4 Relative errors of Qn [F ]. Here f (x) = cos(x)(sin(x) + ex ). n

ω=5

ω = 10

ω = 20

ω = 40

ω = 80

1 2 3 4

2.42 × 10−2 2.70 × 10−3 2.09 × 10−4 1.16 × 10−5

5.39 × 10−3 6.80 × 10−6 1.11 × 10−7 7.14 × 10−10

4.30 × 10−3 5.00 × 10−7 2.10 × 10−9 4.30 × 10−12

1.98 × 10−3 4.06 × 10−9 5.77 × 10−12 3.89 × 10−15

3.15 × 10−3 8.80 × 10−10 2.45 × 10−13 3.21 × 10−17

 k−1       2 ixs ixs xs e W− 2v−k , 2v−k −2ω a + i − ws F a+ 2 2 ω ω s=1 −1 2n+1  a

f (j) (0)  2j Γ ( j+v+1 ) 2 + 1 − Jv (ωt)dt j+1 Γ ( j−v+1 ) j! ω 2 j=0 0   a (2) (2) + j (j + v − 1)Jv (ωa)sj−1,v−1 (ωa) − Jv−1 (ωa)sj,v (ωa) , ω 

n

where F (x) = f (x) −

2n+1 k=0

ixs a+ ω

(3.17)

f (k) (0) k k! x .

Theorem 3.2. The absolute error of the approximate solutions (3.17) for evaluating the exact solutions (3.4) of the Fredholm integral equation with a highly oscillatory Bessel kernel a y(x) −

Jv (ωt)y(t)dt = f (x), 0

satisfies   2   1 y(x) − yn (x) = O (n!) . (2n)! ω 2n+3/2 Proof. From Theorem 2.1, we can easily obtain the result. 2 Example 4. Let us consider the approximate solution yn (x) for 1 y(x) −

  J5 (ωt)y(t)dt = cos(x) sin(x) + ex ,

0

whose exact solution can be represented as 

y(x) = cos(x) sin(x) + e

 x

 +

−1 1

1 1−

J5 (ωt)dt 0

  cos(t) sin(t) + et J5 (ωt)dt

0

(see Table 4 and Fig. 5). In Figs. 4 and 5 we obtain the scaled errors of the proposed methods (3.11) and (3.17) with the different point numbers: Q2 [F ] with node number n = 2 and Q3 [F ] with node number n = 3. From Figs. 4 and 5, we can see that adding the node number n of Gauss–Laguerre rule can improve the accuracy of the proposed methods. At the same time, we can also see that more accurate approximations can be obtained as ω increases.

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1648

Fig. 5. The absolute error scaled by ω 5.5 of Q2 [f ] (left figure), and scaled by ω 7.5 of Q3 [f ] (right figure), and here we choose f (x) = cos(x)(sin(x) + ex ).

In Tables 2–4 we illustrate the relative errors of the formulas (3.11) and (3.17). We can see that increasing ω can lead to more accurate approximations. 4. Concluding remarks In applied problems, one always encounters the problem of evaluating the integral equation with a highly oscillatory Bessel kernel of the forms (1.4) and (1.6). The standard numerical methods fail for solving these problems. In this paper, we derive an explicit formula for the solutions of the first kind integral equations. Preliminary numerical results also show the efficiency and accuracy of the approximations. In fact, the proposed methods can be generalized to other oscillatory integral equations whose exact solutions can be expressed by highly oscillatory integrals. For example, x

  (x − t)n Jn ω(x − t) y(t)dt = f (x),

0

where f (x) is sufficiently differentiable and satisfies f (0) = f  (0) = · · · = f (2n+1) (0) = 0. The solution is given by [30] 2n+2   2n+1 x   d2 n!(n + 1)! 2 2n+1 2 y(x) = (x − t) J2n+1 ω(x − t) +ω f (t)dt. ω (2n)!(2n + 2)! dt2 0

The solution can be rewritten as n!(n + 1)! y(x) = √ (2n)!(2n + 2)! πΓ (2n + 32 )

 1

x (x − t)

4n+2 −1

0

(1 − τ 2 )2n+1 iω(x−t)τ √ e dτ 1 − τ2



d2 + ω2 dt2

2n+2 f (t)dt.

Then the proposed method can be designed to provide approximations accurately. When a = ∞, the exact solutions of the equation ∞ y(x) −

Jv (ωt)y(t)dt = f (x), 0

(4.1)

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

1649

may be denoted by the form  y(x) = f (x) +

−1 ∞

∞ 1−

Jv (ωt)dt 0

f (t)Jv (ωt)dt.

(4.2)

0

In the future we would investigate these problems. Acknowledgments The author is grateful to the associated editor and referees for their valuable comments and suggestions for great improvement of this paper. References [1] M.J. Ablowitz, A.S. Fokas, Complex Variables: Introduction and Applications, Cambridge University Press, 1997. [2] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Dover Publications, New York, 1965. [3] G. Arfken, Mathematical Methods for Physicists, 3rd ed., Academic Press, Orlando, FL, 1985. [4] G. Bao, W. Sun, A fast algorithm for the electromagnetic scattering form a large cavity, SIAM J. Sci. Comput. 27 (2005) 553–574. [5] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, Cambridge University Press, Cambridge, 2004. [6] H. Brunner, A. Iserles, S.P. Nørsett, 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. [7] R.Y. Chen, X.M. Liang, Asymptotic expansions of Bessel, Anger and Weber transformations, J. Math. Anal. Appl. 372 (2010) 377–389. [8] R.Y. Chen, X.M. Liang, Asymptotic expansions for the Bessel transformations with high oscillations by using the homotopy technique, Int. J. Comput. Math. 88 (13) (2011) 2872–2887. [9] K.C. Chung, G.A. Evans, J.R. Webster, A method to generate generalized quadrature rules for oscillatory integrals, Appl. Numer. Math. 34 (2000) 85–93. [10] 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. [11] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, 2nd ed., Academic Press, 1984. [12] F. de Hoog, R. Weiss, On the solution of Volterra integral equations of the first kind, Numer. Math. 21 (1973) 22–32. [13] G.A. Evans, K.C. Chung, Some theoretical aspects of generalised quadrature methods, J. Complexity 19 (2003) 272–285. [14] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, 5th ed., Academic Press, New York, 1994. [15] G. Gripenberg, S.O. Londen, O. Staffans, Volterra Integral and Functional Equations, Cambridge University Press, Cambridge, 1990. [16] J.H. He, A coupling method of homotopy technique and perturbation technique for nonlinear problem, Int. J. Non-Linear Mech. 35 (1) (2000) 37–43. [17] J.H. He, Homotopy perturbation method: a new nonlinear analytical technique, Appl. Math. Comput. 135 (2003) 73–79. [18] D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (2006) 1026–1048. [19] D. Huybrechs, S. Vandewalle, A sparse discretization for integral equation formulations of high frequency scattering problems, SIAM J. Sci. Comput. 29 (2007) 2305–2328. [20] A. Iserles, S.P. Nørsett, On quadrature methods for highly oscillatory integrals and their implementation, BIT 44 (2004) 755–772. [21] H.C. Kang, S.H. Xiang, On the calculation of highly oscillatory integrals with an algebraic singularity, Appl. Math. Comput. 217 (8) (2010) 3890–3897. [22] D. Levin, Fast integration of rapidly oscillatory functions, J. Comput. Appl. Math. 67 (1996) 95–101. [23] D. Levin, Analysis of a collocation method for integrating rapidly oscillatory functions, J. Comput. Appl. Math. 78 (1997) 131–138. [24] P. Linz, Product integration methods for Volterra integral equations of the first kind, BIT 11 (1971) 413–421. [25] Y.L. Luke, Integrals of Bessel Functions, McGraw–Hill Book Company, New York, 1962. [26] L.G. Mcalevey, Product integration rules for Volterra integral equations of the first kind, BIT 27 (1987) 235–247. [27] G.V. Milovanović, Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures, Comput. Math. Appl. 36 (1998) 19–39. [28] S. Olver, Numerical approximation of vector-valued highly oscillatory integrals, BIT 47 (2007) 637–655. [29] R. Piessens, M. Branders, Modified Clewshaw–Curtis method for the computation of Bessel function integrals, BIT 23 (1983) 370–381.

1650

R. Chen / J. Math. Anal. Appl. 421 (2015) 1635–1650

[30] A.D. Polyanin, A.V. Manzhirov, Handbook of Integral Equations, CRC Press, 1998. [31] W. Sun, J. Wu, The superconvergence of Newton–Cotes rules for the Hadamard finite-part integral on an interval, Numer. Math. 109 (2008) 143–165. [32] W. Sun, N.G. Zamani, Adaptive mesh redistribution for the boundary element method in elastostatics, Comput. Struct. 36 (1990) 1081–1088. [33] H.Y. Wang, S.H. Xiang, Asymptotic expansion and Filon-type methods for a Volterra integral equation with a highly oscillatory kernel, IMA J. Numer. Anal. 31 (2) (2011) 469–490. [34] G.N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge University Press, Cambridge, 1966. [35] L.T. Watson, Numerical linear algebra aspects of globally convergent homotopy methods, SIAM Rev. 28 (1986) 529–545. [36] L.T. Watson, C.Y. Wang, A homotopy method applied to elastica problems, Int. J. Solids Struct. 17 (1981) 29–37. [37] R. Wong, Quadrature formula for oscillating integral transforms, Numer. Math. 39 (1982) 351–360. [38] S.H. Xiang, Numerical analysis of a fast integration method for highly oscillatory functions, BIT Numer. Math. 47 (2007) 469–482. [39] S.H. Xiang, Y. Cho, H.Y. Wang, H. Brunner, Clenshaw–Curtis–Filon-type methods for highly oscillatory Bessel transforms and applications, IMA J. Numer. Anal. (2011), http://dx.doi.org/10.1093/imanum/drq035. [40] S.H. Xiang, W.H. Gui, On generalized quadrature rules for fast oscillatory integrals, Appl. Math. Comput. 197 (2008) 60–75. [41] S.H. Xiang, W.H. Gui, P.H. Mo, Numerical quadrature for Bessel transformations, Appl. Numer. Math. 58 (2008) 1247–1261.