- Email: [email protected]
On the evaluation of Cauchy principal value integrals of oscillatory functions
Journal of Computational and Applied Mathematics 234 (2010) 95–100
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the evaluation of Cauchy principal value integrals of oscillatory functionsI Haiyong Wang ∗ , Shuhuang Xiang Department of Mathematics, Central South University, Changsha, Hunan, 410083, PR China
article
abstract
info
Article history: Received 7 April 2008 Received in revised form 31 December 2008
The problem of the numerical evaluation of Cauchy principal value integrals of oscillatory R1 f (x) functions −1 eiωx x−τ dx, where −1 < τ < 1, has been discussed. Based on analytic continuation, if f is analytic in a sufficiently large complex region G containing [−1, 1], the integrals can be transformed into the problems of integrating two integrals on [0, +∞) with the integrand that does not oscillate, and that decays exponentially fast, which can be efficiently computed by using the Gauss–Laguerre quadrature rule. The validity of the method has been demonstrated in the provision of two numerical experiments and their results. © 2009 Elsevier B.V. All rights reserved.
MSC: 65D30 30E20 Keywords: Cauchy principal value Complex integration method Steepest descent method
1. Introduction In this paper we are concerned with the numerical evaluation of the Cauchy principal value (c.p.v.) integrals of oscillatory functions Iω (f ; τ ) :=
Z
1
ei ω x
−1
f (x) x−τ
Z
ei ω x
= lim
ε→0+
dx
|x−τ |≥ε
f ( x) x−τ
dx,
−1 < τ < 1,
(1.1)
where f is analytic in a sufficiently large region of the complex plane containing [−1, 1]. From [1] we know that the integrals (1.1) exist if the function f satisfies Hölder’s condition on the interval [−1, 1]. There are two practical difficulties in the process of evaluating the integrals, namely the integrand is oscillatory and has a singularity of Cauchy type. Okecha [2] and Capobianco and Criscuolo [3] extensively studied the evaluation of (1.1). In [2] the author proposed to compute numerically (1.1) by replacing the function f (x) with the corresponding Lagrange interpolating polynomial on the zeros of Legendre polynomial. In [3] the authors proposed a numerically stable procedure, which is based on an interpolatory procedure at the zeros of the orthogonal polynomials with respect to a Jacobi weight, and the fact that orthogonal polynomials satisfy a three-term recurrence relation to compute the integrals. For more details one can refer to [3,2]. Recently, Okecha proposed an alternative approach in [4] to deal with the c.p.v. integrals of oscillatory functions by using the Hermite interpolation. For details one can see [4].
I The project was sponsored partly by the NSF of China (No.10771218) and partly by the Program for New Century Excellent Talents in University, State Education Ministry, China. ∗ Corresponding author. E-mail address: [email protected] (H. Wang).
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.12.007
96
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
Fig. 1. Integration path for (1.1).
Complex integration methods [5] and steepest descent methods [6] for the finite Fourier integrals of the form
Z
1
g (x)eiωx dx,
(1.2)
−1
have been developed. If g is analytic in a sufficiently large complex region G that contains [−1, 1], the integrals (1.2) can be transformed into the problems of integrating two integrals on [0, +∞), where the integrands are not oscillatory and decay exponentially fast [6,5]. In the present paper, we extend the complex integration method [5] and steepest descent method [6] to the evaluation of the c.p.v. integrals (1.1). Here we use a semicircular contour in the neighborhood of the singularity τ (see Fig. 1) and obtain a new formula for (1.1), which can be efficiently computed by using the Gauss–Laguerre quadrature rule. We also show that the requirement of f can be relaxed comparing with [5] (see (2.1)). Preliminary numerical results show that the precision of the proposed method greatly increases as ω grows, compared with Okecha’s method [2]. 2. Main results Let G denote the region G = {z ∈ C| − 1 ≤ <(z ) ≤ 1, 0 ≤ =(z ) ≤ R} , and G0 = {|z − τ | ≤ r , 0 ≤ arg(z ) ≤ π }, where R is a large number and r is small enough such that G contains G0 (see Fig. 1). Theorem 2.1. Suppose that f is an analytic function in the half-strip of the complex plane, −1 ≤ <(z ) ≤ 1, =(z ) ≥ 0. If there are two constants M and ω0 such that for 0 ≤ ω0 < ω
Z
1
|f (x + iR)|dx ≤ MReω0 R ,
(2.1)
−1
then the c.p.v. integrals of oscillatory functions Iω (f ; τ ) can be transformed into the following form iωτ
Iω (f ; τ ) = ie
π f (τ ) +
ie−iω
ω
it f −1 + ω ieiω dt − ω −1 − τ + ωit
∞
Z
e 0
−t
e 0
it f 1+ ω dt . it 1−τ + ω
∞
Z
−t
(2.2)
f (z )
Proof. Since the integrand eiωz z −τ is analytic in the region G except G0 , by the Cauchy Residue Theorem, we have
Z
ei ω z Γ1 +Γ2 +Γ3 +Γ4 +Γ5 +Γ6
f (z ) z−τ
= 0,
(2.3)
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
97
with all the contours taken in the counterclockwise direction as depicted in Fig. 1. Next, we parameterize each of these integrals in a specific way. Observe
Z
eiωz Γ1
f (z ) z−τ
R
Z
eiω(1+ip)
dz = 0
e−ωp 0
=
ieiω
ωR
Z
f (1 + ip) 1 − τ + ip
dp
it f 1+ ω dt , it 1−τ + ω
e
ω
idp
1 − τ + ip
R
Z
= ieiω
f (1 + ip)
−t
0
(2.4)
with t = ωp in the last expression. Similarly, we have
Z
eiωz Γ3
f (z ) z−τ
R
Z dz = −
eiω(−1+ip)
0 R
Z
= −ie−iω ie−iω
−1 − τ + ip
idp
f (−1 + ip) dp −1 − t + ip it f −1 + ω −t
e−ωp
0
=−
f (−1 + ip)
ωR
Z
e
ω
−1 − τ +
0
it
dt .
(2.5)
ω
From Eq. (2.1) it follows that
Z Z 1 iωz f (z ) iω(x+iR) f (x + iR) = − e dz e dx z−τ x − τ + iR Γ2 −1 Z 1 f (x + iR) dx ≤ e−ωR x − τ + iR −1 Z 1 e−ωR ≤ |f (x + iR)|dx R
−1
≤ Me(ω0 −ω)R → 0,
as R → ∞.
Since the path Γ5 is a half circle enclosing z = τ , we evaluate z − τ = reiθ ,
(2.6)
R
Γ5
f (z )
eiωz z −τ dz by letting
0 ≤ θ ≤ π,
in which case
Z
eiωz Γ5
f (z ) z−τ
π
Z dz = − 0
= −ieiωτ
f (τ + reiθ ) iθ reiθ idθ eiω(τ +re ) τ − τ + reiθ π
Z
iθ
eiωre f (τ + reiθ )dθ .
(2.7)
0
Note that f (z ) is analytic in G. Then f (z ) is continuous at the point τ , namely
|f (z ) − f (τ )| → 0, whenever |z − τ | → 0. Therefore
Z
0
π
e
iωreiθ
Z π iωreiθ iθ f (τ + re )dθ − π f (τ ) = (e f (τ + re ) − f (τ ))dθ Z0 π iωreiθ iθ iωreiθ iωreiθ = (e f (τ + re ) − e f (τ ) + e f (τ ) − f (τ ))dθ 0 Z π Z π iθ ≤ e−ωr sin θ |f (τ + reiθ ) − f (τ )|dθ + |f (τ )| |eiωre − 1|dθ iθ
0
Z ≤ 0
0
π
|f (τ + reiθ ) − f (τ )|dθ + |f (τ )|
π
Z 0
iθ
|eiωre − 1|dθ .
(2.8)
98
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100 iθ
It is evident that eiωre → 1 as r → 0, and it is true that
Z
π
e
iωreiθ
0
f (τ + re )dθ − π f (τ ) → 0, iθ
r → 0.
(2.9)
Recalling the definition of the c.p.v. integrals and by means of (2.3), when R → ∞ and r → 0, we obtain
Z
Iω (f ; τ ) = lim
ε→0+
x−τ
|x−τ |≥ε
Z =
f ( x)
ei ω x
eiωz
lim
r →0,R→∞
Γ4 +Γ6
dx
f (z ) z−τ
Z =−
ei ω z
lim
r →0,R→∞
dz
Γ1 +Γ2 +Γ3 +Γ5
ie−iω
= ieiωτ π f (τ ) +
ω
dz
it f −1 + ω ieiω dt − it ω −1 − τ + ω
∞
Z
f (z ) z−τ
e− t 0
where −1 < τ < 1. This completes the proof.
Z
it f 1+ ω dt , it 1−τ + ω
∞
e− t
0
(2.10)
In analogy to [6,5], the obtained integrals in Theorem 2.1 can be evaluated by using the Gauss–Laguerre quadrature rule. Let xk and wk denote the nodes and weights of the n-point Gauss–Laguerre quadrature rule; then the c.p.v. integrals of oscillatory functions can be approximated by n ie−iω X
Iω (f ; τ ) ≈ Qn (f , τ ) = ieiωτ π f (τ ) +
ω
wk
k=1
f
ixk
−1 +
ω
−1 − τ +
ixk
n ieiω X
−
ω
ω
wk
k=1
f
ix
1 + ωk
1 − τ + ωk ix
.
(2.11)
We can also approximate the following two kinds of integrals
Z
1
cos(ωx) −1
Z
1
sin(ωx) −1
f (x) x−τ f (x)
eiωx
−1 1
Z
x−τ
1
Z dx = < dx = =
e
iω x
−1
f (x) x−τ f ( x)
x−τ
dx ≈ <(Qn (f ; τ )), (2.12) dx ≈ =(Qn (f ; τ )).
Suppose that f satisfies R ∞the conditions of Theorem 2.1. The error formula for the n-point Gauss–Laguerre quadrature rule applied to the integrals 0 e−x f (x)dx is given in [1]
(n!)2 (2n) f (ξ ), (2n)!
0 < ξ < ∞.
Using this formula and from (2.2) and (2.11), we obtain Iω (f ; τ ) − Qn (f ; τ ) =
=
ie−iω
Z
∞
ix it n f −1 + ωk X f −1 + ω −t dt − wk e −1 − τ + ωit −1 − τ + ixωk k=1
ω 0 ix it n Z ∞ f 1 + ωk X f 1 + ieiω ω − e− t dt − wk it ix ω 0 1−τ + ω 1−τ + k k=1 ie−iω (n!)2
d
2n
it f −1 + ω
it −1−τ + ω
ω (2n)!
dt 2n
|t =ξ1 −
d
ieiω (n!)2
2n
ω
it f 1+ ω
it 1−τ + ω
ω (2n)!
dt 2n
|t =ξ2 ,
(2.13)
with ξ1 , ξ2 ∈ C. If ω 1, from [6] we know that the error decays asymptotically as O(ω−2n−1 ). Hence more accurate approximations can be obtained as ω increases. In the next section, we present two examples to illustrate the order of the error. Remark 1. If f has v singularities zk (k = 1, 2, . . . , v) inside the region G. We can take r sufficiently small so that all Pv f (z ) singularities lie inside the region G except G0 . Let a = 2π i k=1 ak , where ak is the residue of eiωz z −τ at z = zk . In a similar way, Theorem 2.1 is also true by replacing Eqs. (2.2) and (2.3) with Iω (f ; τ ) = a + ie
iωτ
π f (τ ) +
ie−iω
ω
Z
e 0
it f −1 + ω ieiω dt − it ω −1 − τ + ω
∞ −t
e 0
it f 1+ ω dt , it 1−τ + ω
∞
Z
−t
(2.14)
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
99
Table 1 x
R1
Absolute errors in n-point approximations by the Gauss–Laguerre quadrature to the integral −1 sin(ωx) ex dx. The last row shows the value of log2 (e32 /e64 ).
ω
n
8 16 32 64 Rate
3
4
5
9.55 × 10−6 8.02 × 10−8 1.89 × 10−9 1.14 × 10−11 7.4
8.83 × 10−7 3.88 × 10−9 2.18 × 10−11 5.26 × 10−14 8.7
2.25 × 10−7 6.04 × 10−10 1.95 × 10−13 3.50 × 10−16 9.1
and
Z
ei ω z Γ1 +Γ2 +Γ3 +Γ4 +Γ5 +Γ6
f (z ) z−τ
= a,
(2.15)
respectively. Remark 2. It is sometimes advantageous to express the error of the Gauss–Laguerre rule in the form
(n!)2 f [ξ , x1 , x1 , x2 , x2 , . . . , xn , xn ]
(2.16)
where 0 < ξ < ∞ and f [ξ , x1 , x1 , x2 , x2 , . . . , xn , xn ] is the 2n-th divided difference of f (x), relative to the abscissas x1 , x1 , . . . , xn , xn and ξ . The corresponding error term in (2.13) can hence be written as ie−iω
Iω (f ; τ ) − Qn (f ; τ ) =
ω
(n!)2 f [ξ1 , x1 , x2 , . . . , xn , xn ] −
ieiω
ω
(n!)2 f [ξ2 , x1 , x2 , . . . , xn , xn ]
(2.17)
where ξ1 , ξ2 ∈ C. 3. Numerical examples In order to illustrate the quality of the approximation provided by (2.11), we test two numerical examples. The nodes xk and weights wk of the Gauss–Laguerre quadrature rule are given in [7, pp. 923]. All computations in this section have been performed by using Maple 11 with 32-digit arithmetic. x
R1
Example 1 ([2,4]). The computation of =(Iω (ex , 0)) = −1 sin(ωx) ex dx: from [2], the analytic expression of the integral can be obtained by taking the series expansion of the exponential function and using the results in [8, pp. 226]
Z
1
sin(ωx) −1
ex x
dx = 2Si(ω) − 2
∞ X k=1
X 2k − 1 1 jπ j! cos ω + , (2k)! j=0 j ω j +1 2 1
2k−1
(3.1)
where Si is the sine integral. Notice that for this integral τ = 0, f (x) = ex . It is evident that f is analytic in the complex plane and 1
Z
|f (x + iR)|dx =
Z
−1
1
|ex+iR |dx = e − e−1 < Reω0 R ,
R 1.
−1
Therefore the proposed method (2.11) can be applied. Table 1 demonstrates the absolute errors in n-point approximations to the integral. It exhibits the fast convergence of the approximations as n increases. The table also shows that more accurate approximations can be obtained as ω increases. Following the example given in [2,4], we take ω = 12 for example. To achieve the accuracy 10−13 , the Gauss–Laguerre rule requires 10 points and Okecha’s method requires 15 points. Example 2 ([2]). The computation of =(Iω (sinh(x), −.13)) = expression of the integral
Z
1
−1
eiωx
sinh(x) x−τ
dx =
R1 −1
sinh(x)
sin(ωx) x+.13 dx: from [2], we obtain the analytic
1 (1+iω)τ e (Ei((1 + iω)(1 − τ )) − Ei((1 + iω)(−1 − τ )) − iπ ) 2
− e(−1+iω)τ (Ei((−1 + iω)(1 − τ )) − Ei((−1 + iω)(−1 − τ )) − iπ ) , where Ei is the exponential integral. Note that for this integral τ = −0.13 and f (x) = sinh(x). It is evident that f is analytic in the complex plane and
Z
1
|f (x + iR)|dx = −1
Z
1
−1
x+iR Z e − e−x−iR dx ≤ 2
1
−1
ex + e− x 2
dx = e − e−1 < Reω0 R ,
Therefore the integral can be handled by the proposed method (2.11).
R 1.
(3.2)
100
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
Table 2
R1
sinh(x)
Absolute errors in n-point approximations by the Gauss–Laguerre quadrature to the integral −1 sin(ωx) x+.13 dx. The last row shows the value of log2 (e64 /e128 ).
ω 16 32 64 128 Rate
n 3
4
5
1.04 × 10−8 3.64 × 10−10 2.36 × 10−12 1.33 × 10−14 7.5
1.23 × 10−9 4.59 × 10−12 1.35 × 10−14 1.21 × 10−17 10.1
1.53 × 10−10 2.05 × 10−14 1.11 × 10−16 1.33 × 10−20 13.0
In Table 2 we present the absolute errors in n-point Gauss–Laguerre approximations. As we can see, the accuracy improves as ω increases, and the convergence is fast and we can take relatively small n in order to get a satisfactory result. 4. Conclusion The numerical evaluation of the c.p.v. integrals of oscillatory functions was considered by using the technique of interpolatory type in [3,2,4]. In this paper, we have presented a method for handling (1.2) based on the analyticity of the integrand. This approach is easy to be implemented by using the Gauss–Laguerre quadrature rule. It may be observed from Tables 1 and 2 that we can obtain quite rapid convergence for small, moderate as well as large values of ω. Furthermore, the accuracy of the proposed method increases both for the case of a fixed number of points and increasing ω, and for the case of an increasing number of points and fixed ω. Acknowledgements The authors are deeply grateful to the anonymous referees for their valuable comments and suggestions for the improvement of this paper. The first author acknowledges the financial support from the Hunan Provincial Innovation Foundation For Postgraduate. References [1] [2] [3] [4] [5] [6] [7] [8]
P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edition, Academic Press, 1984. G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comp. 49 (1987) 259–268. M.R. Capobianco, G. Criscuolo, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. G.E. Okecha, Hermite interpolation and a method for evaluating Cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91–98. G.V. Milovanović, Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures, Comput. Math. Appl. 36 (1998) 19–39. D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (2006) 1026–1048. M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series, and Products, 6th ed., Academic Press, San Diego, CA, 2000.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the evaluation of Cauchy principal value integrals of oscillatory functionsI Haiyong Wang ∗ , Shuhuang Xiang Department of Mathematics, Central South University, Changsha, Hunan, 410083, PR China
article
abstract
info
Article history: Received 7 April 2008 Received in revised form 31 December 2008
The problem of the numerical evaluation of Cauchy principal value integrals of oscillatory R1 f (x) functions −1 eiωx x−τ dx, where −1 < τ < 1, has been discussed. Based on analytic continuation, if f is analytic in a sufficiently large complex region G containing [−1, 1], the integrals can be transformed into the problems of integrating two integrals on [0, +∞) with the integrand that does not oscillate, and that decays exponentially fast, which can be efficiently computed by using the Gauss–Laguerre quadrature rule. The validity of the method has been demonstrated in the provision of two numerical experiments and their results. © 2009 Elsevier B.V. All rights reserved.
MSC: 65D30 30E20 Keywords: Cauchy principal value Complex integration method Steepest descent method
1. Introduction In this paper we are concerned with the numerical evaluation of the Cauchy principal value (c.p.v.) integrals of oscillatory functions Iω (f ; τ ) :=
Z
1
ei ω x
−1
f (x) x−τ
Z
ei ω x
= lim
ε→0+
dx
|x−τ |≥ε
f ( x) x−τ
dx,
−1 < τ < 1,
(1.1)
where f is analytic in a sufficiently large region of the complex plane containing [−1, 1]. From [1] we know that the integrals (1.1) exist if the function f satisfies Hölder’s condition on the interval [−1, 1]. There are two practical difficulties in the process of evaluating the integrals, namely the integrand is oscillatory and has a singularity of Cauchy type. Okecha [2] and Capobianco and Criscuolo [3] extensively studied the evaluation of (1.1). In [2] the author proposed to compute numerically (1.1) by replacing the function f (x) with the corresponding Lagrange interpolating polynomial on the zeros of Legendre polynomial. In [3] the authors proposed a numerically stable procedure, which is based on an interpolatory procedure at the zeros of the orthogonal polynomials with respect to a Jacobi weight, and the fact that orthogonal polynomials satisfy a three-term recurrence relation to compute the integrals. For more details one can refer to [3,2]. Recently, Okecha proposed an alternative approach in [4] to deal with the c.p.v. integrals of oscillatory functions by using the Hermite interpolation. For details one can see [4].
I The project was sponsored partly by the NSF of China (No.10771218) and partly by the Program for New Century Excellent Talents in University, State Education Ministry, China. ∗ Corresponding author. E-mail address: [email protected] (H. Wang).
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.12.007
96
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
Fig. 1. Integration path for (1.1).
Complex integration methods [5] and steepest descent methods [6] for the finite Fourier integrals of the form
Z
1
g (x)eiωx dx,
(1.2)
−1
have been developed. If g is analytic in a sufficiently large complex region G that contains [−1, 1], the integrals (1.2) can be transformed into the problems of integrating two integrals on [0, +∞), where the integrands are not oscillatory and decay exponentially fast [6,5]. In the present paper, we extend the complex integration method [5] and steepest descent method [6] to the evaluation of the c.p.v. integrals (1.1). Here we use a semicircular contour in the neighborhood of the singularity τ (see Fig. 1) and obtain a new formula for (1.1), which can be efficiently computed by using the Gauss–Laguerre quadrature rule. We also show that the requirement of f can be relaxed comparing with [5] (see (2.1)). Preliminary numerical results show that the precision of the proposed method greatly increases as ω grows, compared with Okecha’s method [2]. 2. Main results Let G denote the region G = {z ∈ C| − 1 ≤ <(z ) ≤ 1, 0 ≤ =(z ) ≤ R} , and G0 = {|z − τ | ≤ r , 0 ≤ arg(z ) ≤ π }, where R is a large number and r is small enough such that G contains G0 (see Fig. 1). Theorem 2.1. Suppose that f is an analytic function in the half-strip of the complex plane, −1 ≤ <(z ) ≤ 1, =(z ) ≥ 0. If there are two constants M and ω0 such that for 0 ≤ ω0 < ω
Z
1
|f (x + iR)|dx ≤ MReω0 R ,
(2.1)
−1
then the c.p.v. integrals of oscillatory functions Iω (f ; τ ) can be transformed into the following form iωτ
Iω (f ; τ ) = ie
π f (τ ) +
ie−iω
ω
it f −1 + ω ieiω dt − ω −1 − τ + ωit
∞
Z
e 0
−t
e 0
it f 1+ ω dt . it 1−τ + ω
∞
Z
−t
(2.2)
f (z )
Proof. Since the integrand eiωz z −τ is analytic in the region G except G0 , by the Cauchy Residue Theorem, we have
Z
ei ω z Γ1 +Γ2 +Γ3 +Γ4 +Γ5 +Γ6
f (z ) z−τ
= 0,
(2.3)
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
97
with all the contours taken in the counterclockwise direction as depicted in Fig. 1. Next, we parameterize each of these integrals in a specific way. Observe
Z
eiωz Γ1
f (z ) z−τ
R
Z
eiω(1+ip)
dz = 0
e−ωp 0
=
ieiω
ωR
Z
f (1 + ip) 1 − τ + ip
dp
it f 1+ ω dt , it 1−τ + ω
e
ω
idp
1 − τ + ip
R
Z
= ieiω
f (1 + ip)
−t
0
(2.4)
with t = ωp in the last expression. Similarly, we have
Z
eiωz Γ3
f (z ) z−τ
R
Z dz = −
eiω(−1+ip)
0 R
Z
= −ie−iω ie−iω
−1 − τ + ip
idp
f (−1 + ip) dp −1 − t + ip it f −1 + ω −t
e−ωp
0
=−
f (−1 + ip)
ωR
Z
e
ω
−1 − τ +
0
it
dt .
(2.5)
ω
From Eq. (2.1) it follows that
Z Z 1 iωz f (z ) iω(x+iR) f (x + iR) = − e dz e dx z−τ x − τ + iR Γ2 −1 Z 1 f (x + iR) dx ≤ e−ωR x − τ + iR −1 Z 1 e−ωR ≤ |f (x + iR)|dx R
−1
≤ Me(ω0 −ω)R → 0,
as R → ∞.
Since the path Γ5 is a half circle enclosing z = τ , we evaluate z − τ = reiθ ,
(2.6)
R
Γ5
f (z )
eiωz z −τ dz by letting
0 ≤ θ ≤ π,
in which case
Z
eiωz Γ5
f (z ) z−τ
π
Z dz = − 0
= −ieiωτ
f (τ + reiθ ) iθ reiθ idθ eiω(τ +re ) τ − τ + reiθ π
Z
iθ
eiωre f (τ + reiθ )dθ .
(2.7)
0
Note that f (z ) is analytic in G. Then f (z ) is continuous at the point τ , namely
|f (z ) − f (τ )| → 0, whenever |z − τ | → 0. Therefore
Z
0
π
e
iωreiθ
Z π iωreiθ iθ f (τ + re )dθ − π f (τ ) = (e f (τ + re ) − f (τ ))dθ Z0 π iωreiθ iθ iωreiθ iωreiθ = (e f (τ + re ) − e f (τ ) + e f (τ ) − f (τ ))dθ 0 Z π Z π iθ ≤ e−ωr sin θ |f (τ + reiθ ) − f (τ )|dθ + |f (τ )| |eiωre − 1|dθ iθ
0
Z ≤ 0
0
π
|f (τ + reiθ ) − f (τ )|dθ + |f (τ )|
π
Z 0
iθ
|eiωre − 1|dθ .
(2.8)
98
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100 iθ
It is evident that eiωre → 1 as r → 0, and it is true that
Z
π
e
iωreiθ
0
f (τ + re )dθ − π f (τ ) → 0, iθ
r → 0.
(2.9)
Recalling the definition of the c.p.v. integrals and by means of (2.3), when R → ∞ and r → 0, we obtain
Z
Iω (f ; τ ) = lim
ε→0+
x−τ
|x−τ |≥ε
Z =
f ( x)
ei ω x
eiωz
lim
r →0,R→∞
Γ4 +Γ6
dx
f (z ) z−τ
Z =−
ei ω z
lim
r →0,R→∞
dz
Γ1 +Γ2 +Γ3 +Γ5
ie−iω
= ieiωτ π f (τ ) +
ω
dz
it f −1 + ω ieiω dt − it ω −1 − τ + ω
∞
Z
f (z ) z−τ
e− t 0
where −1 < τ < 1. This completes the proof.
Z
it f 1+ ω dt , it 1−τ + ω
∞
e− t
0
(2.10)
In analogy to [6,5], the obtained integrals in Theorem 2.1 can be evaluated by using the Gauss–Laguerre quadrature rule. Let xk and wk denote the nodes and weights of the n-point Gauss–Laguerre quadrature rule; then the c.p.v. integrals of oscillatory functions can be approximated by n ie−iω X
Iω (f ; τ ) ≈ Qn (f , τ ) = ieiωτ π f (τ ) +
ω
wk
k=1
f
ixk
−1 +
ω
−1 − τ +
ixk
n ieiω X
−
ω
ω
wk
k=1
f
ix
1 + ωk
1 − τ + ωk ix
.
(2.11)
We can also approximate the following two kinds of integrals
Z
1
cos(ωx) −1
Z
1
sin(ωx) −1
f (x) x−τ f (x)
eiωx
−1 1
Z
x−τ
1
Z dx = < dx = =
e
iω x
−1
f (x) x−τ f ( x)
x−τ
dx ≈ <(Qn (f ; τ )), (2.12) dx ≈ =(Qn (f ; τ )).
Suppose that f satisfies R ∞the conditions of Theorem 2.1. The error formula for the n-point Gauss–Laguerre quadrature rule applied to the integrals 0 e−x f (x)dx is given in [1]
(n!)2 (2n) f (ξ ), (2n)!
0 < ξ < ∞.
Using this formula and from (2.2) and (2.11), we obtain Iω (f ; τ ) − Qn (f ; τ ) =
=
ie−iω
Z
∞
ix it n f −1 + ωk X f −1 + ω −t dt − wk e −1 − τ + ωit −1 − τ + ixωk k=1
ω 0 ix it n Z ∞ f 1 + ωk X f 1 + ieiω ω − e− t dt − wk it ix ω 0 1−τ + ω 1−τ + k k=1 ie−iω (n!)2
d
2n
it f −1 + ω
it −1−τ + ω
ω (2n)!
dt 2n
|t =ξ1 −
d
ieiω (n!)2
2n
ω
it f 1+ ω
it 1−τ + ω
ω (2n)!
dt 2n
|t =ξ2 ,
(2.13)
with ξ1 , ξ2 ∈ C. If ω 1, from [6] we know that the error decays asymptotically as O(ω−2n−1 ). Hence more accurate approximations can be obtained as ω increases. In the next section, we present two examples to illustrate the order of the error. Remark 1. If f has v singularities zk (k = 1, 2, . . . , v) inside the region G. We can take r sufficiently small so that all Pv f (z ) singularities lie inside the region G except G0 . Let a = 2π i k=1 ak , where ak is the residue of eiωz z −τ at z = zk . In a similar way, Theorem 2.1 is also true by replacing Eqs. (2.2) and (2.3) with Iω (f ; τ ) = a + ie
iωτ
π f (τ ) +
ie−iω
ω
Z
e 0
it f −1 + ω ieiω dt − it ω −1 − τ + ω
∞ −t
e 0
it f 1+ ω dt , it 1−τ + ω
∞
Z
−t
(2.14)
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
99
Table 1 x
R1
Absolute errors in n-point approximations by the Gauss–Laguerre quadrature to the integral −1 sin(ωx) ex dx. The last row shows the value of log2 (e32 /e64 ).
ω
n
8 16 32 64 Rate
3
4
5
9.55 × 10−6 8.02 × 10−8 1.89 × 10−9 1.14 × 10−11 7.4
8.83 × 10−7 3.88 × 10−9 2.18 × 10−11 5.26 × 10−14 8.7
2.25 × 10−7 6.04 × 10−10 1.95 × 10−13 3.50 × 10−16 9.1
and
Z
ei ω z Γ1 +Γ2 +Γ3 +Γ4 +Γ5 +Γ6
f (z ) z−τ
= a,
(2.15)
respectively. Remark 2. It is sometimes advantageous to express the error of the Gauss–Laguerre rule in the form
(n!)2 f [ξ , x1 , x1 , x2 , x2 , . . . , xn , xn ]
(2.16)
where 0 < ξ < ∞ and f [ξ , x1 , x1 , x2 , x2 , . . . , xn , xn ] is the 2n-th divided difference of f (x), relative to the abscissas x1 , x1 , . . . , xn , xn and ξ . The corresponding error term in (2.13) can hence be written as ie−iω
Iω (f ; τ ) − Qn (f ; τ ) =
ω
(n!)2 f [ξ1 , x1 , x2 , . . . , xn , xn ] −
ieiω
ω
(n!)2 f [ξ2 , x1 , x2 , . . . , xn , xn ]
(2.17)
where ξ1 , ξ2 ∈ C. 3. Numerical examples In order to illustrate the quality of the approximation provided by (2.11), we test two numerical examples. The nodes xk and weights wk of the Gauss–Laguerre quadrature rule are given in [7, pp. 923]. All computations in this section have been performed by using Maple 11 with 32-digit arithmetic. x
R1
Example 1 ([2,4]). The computation of =(Iω (ex , 0)) = −1 sin(ωx) ex dx: from [2], the analytic expression of the integral can be obtained by taking the series expansion of the exponential function and using the results in [8, pp. 226]
Z
1
sin(ωx) −1
ex x
dx = 2Si(ω) − 2
∞ X k=1
X 2k − 1 1 jπ j! cos ω + , (2k)! j=0 j ω j +1 2 1
2k−1
(3.1)
where Si is the sine integral. Notice that for this integral τ = 0, f (x) = ex . It is evident that f is analytic in the complex plane and 1
Z
|f (x + iR)|dx =
Z
−1
1
|ex+iR |dx = e − e−1 < Reω0 R ,
R 1.
−1
Therefore the proposed method (2.11) can be applied. Table 1 demonstrates the absolute errors in n-point approximations to the integral. It exhibits the fast convergence of the approximations as n increases. The table also shows that more accurate approximations can be obtained as ω increases. Following the example given in [2,4], we take ω = 12 for example. To achieve the accuracy 10−13 , the Gauss–Laguerre rule requires 10 points and Okecha’s method requires 15 points. Example 2 ([2]). The computation of =(Iω (sinh(x), −.13)) = expression of the integral
Z
1
−1
eiωx
sinh(x) x−τ
dx =
R1 −1
sinh(x)
sin(ωx) x+.13 dx: from [2], we obtain the analytic
1 (1+iω)τ e (Ei((1 + iω)(1 − τ )) − Ei((1 + iω)(−1 − τ )) − iπ ) 2
− e(−1+iω)τ (Ei((−1 + iω)(1 − τ )) − Ei((−1 + iω)(−1 − τ )) − iπ ) , where Ei is the exponential integral. Note that for this integral τ = −0.13 and f (x) = sinh(x). It is evident that f is analytic in the complex plane and
Z
1
|f (x + iR)|dx = −1
Z
1
−1
x+iR Z e − e−x−iR dx ≤ 2
1
−1
ex + e− x 2
dx = e − e−1 < Reω0 R ,
Therefore the integral can be handled by the proposed method (2.11).
R 1.
(3.2)
100
H. Wang, S. Xiang / Journal of Computational and Applied Mathematics 234 (2010) 95–100
Table 2
R1
sinh(x)
Absolute errors in n-point approximations by the Gauss–Laguerre quadrature to the integral −1 sin(ωx) x+.13 dx. The last row shows the value of log2 (e64 /e128 ).
ω 16 32 64 128 Rate
n 3
4
5
1.04 × 10−8 3.64 × 10−10 2.36 × 10−12 1.33 × 10−14 7.5
1.23 × 10−9 4.59 × 10−12 1.35 × 10−14 1.21 × 10−17 10.1
1.53 × 10−10 2.05 × 10−14 1.11 × 10−16 1.33 × 10−20 13.0
In Table 2 we present the absolute errors in n-point Gauss–Laguerre approximations. As we can see, the accuracy improves as ω increases, and the convergence is fast and we can take relatively small n in order to get a satisfactory result. 4. Conclusion The numerical evaluation of the c.p.v. integrals of oscillatory functions was considered by using the technique of interpolatory type in [3,2,4]. In this paper, we have presented a method for handling (1.2) based on the analyticity of the integrand. This approach is easy to be implemented by using the Gauss–Laguerre quadrature rule. It may be observed from Tables 1 and 2 that we can obtain quite rapid convergence for small, moderate as well as large values of ω. Furthermore, the accuracy of the proposed method increases both for the case of a fixed number of points and increasing ω, and for the case of an increasing number of points and fixed ω. Acknowledgements The authors are deeply grateful to the anonymous referees for their valuable comments and suggestions for the improvement of this paper. The first author acknowledges the financial support from the Hunan Provincial Innovation Foundation For Postgraduate. References [1] [2] [3] [4] [5] [6] [7] [8]
P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edition, Academic Press, 1984. G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comp. 49 (1987) 259–268. M.R. Capobianco, G. Criscuolo, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. G.E. Okecha, Hermite interpolation and a method for evaluating Cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91–98. G.V. Milovanović, Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures, Comput. Math. Appl. 36 (1998) 19–39. D. Huybrechs, S. Vandewalle, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (2006) 1026–1048. M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1964. I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series, and Products, 6th ed., Academic Press, San Diego, CA, 2000.