- Email: [email protected]
Evaluation of Cauchy principal value integrals of oscillatory kind
Applied Mathematics and Computation 217 (2010) 2390–2396
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Evaluation of Cauchy principal value integrals of oscillatory kind Jianbing Li *, Xuesong Wang, Tao Wang College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, PR China
a r t i c l e
i n f o
Keywords: Cauchy principal value integral Oscillatory integral Improved-Levin method Chebyshev differentiation matrix Gauss–Laguerre quadrature
a b s t r a c t This paper studies the evaluation of Cauchy principal value (c.p.v.) integrals of oscillatory R b f ðxÞ ixx kind: --a x s e dx, where s 2 (a, b). The singularity of the c.p.v. integral is transferred to an individual oscillatory integral that can be analytically calculated. The remaining integral which is non-singular can be well calculated with a stable and fast quadrature method, namely the improved-Levin method. Numerical examples confirm the advantages of the proposed method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction The rapid and accurate evaluation of oscillatory integrals is one of the major issues arising in many fields, especially in the computational electromagnetics, optics, and quantum mechanics [1]. Efforts have been made in the past several decades on this issue, and several efficient quadrature methods have been developed. These methods include the asymptotic expansion [2], the Filon(-type) method [3,4], the Levin(-type) method [5,6], the (numerical-)steepest descent method [1,7], etc. [8]. The existing methods are mainly concerned with oscillatory integrals with common integrands. However, in many situations of practical interest, the oscillatory integrals may be singular. For example, the Green function in computational electromagnetics is both oscillatory and singular [1]. Therefore, special attention should be paid on these integrals. In this paper, we study the evaluation of a kind of oscillatory integral with a Cauchy type singularity:
Z Ix ðf ; s; a; bÞ :¼ --
b
a
f ðxÞ ixx e dx; xs
s 2 ða; bÞ; x 0;
ð1Þ
where f(x) is a smooth and non-oscillatory function. Some work has been done on this issue. In [9], Okecha calculates the integral (1) by replacing the function f(x) with the corresponding Lagrange interpolating polynomials on the zeros of the Legendre polynomial and on the singularity point s. The algorithm was improved in [10] by introducing the Hermite interpolation polynomial to obtain a higher efficiency. But this algorithm was proved to be not stable enough by Capobianco and Criscuolo in [11], and the new method they proposed is based on an interpolation procedure at the zeros of the orthogonal polynomials with respect to a Jacobi weight. Recently, Haiyong Wang and Shuhuang Xiang proposed a uniform approximations to the integrals in [12]. It replaces the function f(x) by its Lagrange interpolating polynomial based on the Chebyshev– Gauss–Labatto nodes, and transfers the singularity to an integral that can be analytically calculated. In the method proposed by this paper, the original integral is also separated into two oscillatory integrals: one is non-singular and another one is singular. The non-singular one is computed with an improved-Levin method, and very accurate integral result can be obtained. The second integral has an analytical solution, but several special functions (sine integral and cosine integral) have to be computed. We will propose an alternate approach to calculate the second integral, which does not need to compute special functions and could be more efficient.
* Corresponding author. E-mail address: [email protected] (J. Li). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.07.039
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
2391
2. Quadrature of the Cauchy principal integral based on the improved-Levin method 2.1. Transformation of the Cauchy principal integral The original integral can be easily transformed into a uniform type if the coordinate transformation x0 ¼ ba x þ bþa is em2 2 ployed. Therefore, the present work only takes the uniform type into account:
Z1 f ðxÞ Ix ðf ; sÞ :¼ -e 1 x s
ixx
Z
dx ¼ limþ e!0
jxsjPe
f ðxÞ ixx e dx; xs
s 2 ð1; 1Þ:
The singularity of the integral comes from the existence of a pole at x = s. A reasonable approach of removing the singularity is to make the integrand be the form of ‘0/0’ at the pole, so we transform the original integral into the following two integrals:
Z
Ix ðf ; sÞ ¼
1
1
f ðxÞ f ðsÞ e xs
ixx
dx þ f ðsÞ
Z
1
1
eixx dx , Ix ðf ; sÞ þ f ðsÞws;x : xs
ð2Þ
It is well known that the Taylor expansion of f(x) at x = s is
f ðxÞ ¼ f ðsÞ þ
1 X f ðnÞ ðsÞðx sÞn ; n! n¼1
ð3Þ
so the first integral on the right-hand side of (2) can be re-expressed as
Z
Ix ðf ; sÞ ¼
1
1
f ðxÞ f ðsÞ ixx e dx ¼ xs
Z
1
1
1 X f ðnÞ ðsÞðx sÞn1 ixx e dx: n! n¼1
ð4Þ
Obviously, this integral is non-singular. However, the integral in the second integral is singular at x = s:
Z1 ws;x ¼ --
eixx dx: 1 x s
ð5Þ
Fortunately, this integral has an analytical result [13]:
Re½ws;x ¼ cosðxsÞ½Ciðu1 Þ Ciðju2 jÞ sinðxsÞ½Siðu1 Þ þ Siðju2 jÞ; and
Im½ws;x ¼ sinðxsÞ½Ciðu1 Þ Ciðju2 jÞ þ cosðxsÞ½Siðu1 Þ þ Siðju2 jÞ; where u1 = x(1 s), u2 = x(1 + s), and
CiðuÞ ¼ SiðuÞ ¼
Z Z
u
0
0
u
cosðfÞ 1 df þ ln u þ C; f sin f df; f
u > 0;
are respectively the cosine and sine integrals with C being the Euler constant. In the analytical result, one has to compute two special functions: the cosine and sine integrals. These two integrals can be well handled by some mathematical software such as matlab, but the computational complexity is generally not small enough. In this paper, we present an alternate approach to compute the integral ws,x. By choosing a proper integration path in the complex plane, the integral ws,x is transformed into several singular-free and oscillatory-free integrals which can be well calculated by the Gauss–Laguerre quadrature algorithm [14,15]. Because the abscissas and weights of the Gauss–Laguerre quadrature algorithm can be pre-computed, the complexity of the newly proposed approach is generally very small. For more details, please refer to the Appendix A. Because the second term on the right-hand side of (2), ws,x, can be well calculated, we shall pay the main attention on the computation of the first integral Ix ðf ; sÞ. 2.2. Quadrature of Ix ðf ; sÞ by the improved-Levin method ðsÞ Let hðx; sÞ ¼ f ðxÞf . It is immediately derived from the Taylor expansion (or the L’Hôspital rule for the limitation of ‘0/0’ xs type) that
hðs; sÞ ¼ f 0 ðsÞ: The integral Ix ðf ; sÞ is then denoted by
ð6Þ
2392
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Ix ðf ; sÞ ¼
Z
1
hðx; sÞe
ix x
dx;
ð7Þ
1
where h(x, s) is non-singular. We have proposed an improved-Levin quadrature method for oscillatory integrals [16–18]. The present work adopts this quadrature method to calculate the integral (7). In light of the Levin’s theory, the quadrature of integral (7) is reduced to the calculation of the following ordinary differential equation (ODE) without boundary conditions:
p0 ðxÞ þ ixpðxÞ ¼ hðx; sÞ:
ð8Þ
The Chebyshev differentiation matrix, which is a powerful tool in the spectral method [19–22], is used to solve the ODE. An explicit expression of the matrix follows
Dkj ¼
8 > ck ð1Þkþj > > < cj ðxk xj Þ ;
k; j ¼ 0; . . . ; N; k – j; cj ¼
N > P > > Dkn ; :
2 j ¼ 0; N; 1 j ¼ 1; . . . N 1;
ð9Þ
k ¼ j;
n¼0;n–k
where N + 1 is the number of nodes, xj ¼ cos pNj ; j ¼ 0; 1; . . . ; N is called the Chebyshev–Gauss–Labatto nodes. If the function values of f(x) on the Chebyshev–Gauss–Labatto nodes compose a vector f, then the function values of f0 (x) on these nodes can be approximated by f0 = Df with very high accuracy. Therefore, applying the Chebyshev differentiation matrix on ODE (8) yields
ðD þ ixEÞp ¼ H;
ð10Þ
where E is an identity matrix, p and H are the numerical vectors of p(x) and h(x, s) on the Chebyshev–Gauss–Labatto nodes, respectively, i.e.,
p ¼ ½pðx0 Þ; pðx1 Þ; ; pðxN ÞT ;
H ¼ ½hðx0 ; sÞ; hðx1 ; sÞ; ; hðxN ; sÞT :
The determination of H is a little special. For a given node xj, if xj – s, the component h(xj, s) is computed with
hðxj ; sÞ ¼
f ðxj Þ f ðsÞ ; xj s
ð11Þ
otherwise, the component should be computed according to (6). This treatment does not truncate the Taylor series (3), so it could lead to an exact numerical vector H. It should be also noted that derivatives are not needed in closed form because the differentiation matrix D can be applied to compute the derivatives. The behavior and the proper solution method of (10) have been comprehensively studied in [16–18]. From this system of linear equations, one can easily work out the solution p, and then the integral is yielded as
Ix ðf ; sÞ ¼ pð1Þeix pðN þ 1Þeix : The differentiation matrix D is generally pre-computed, so the determination of the system of linear Eq. (10) is with very small computational complexity. 3. Error analysis The error of calculating the oscillatory integrals via the improved-Levin method may comes from two aspects: the determination of the target system of linear equations and the computation of it. 1. For the aspect of determining the system of linear equations: As has been stated, the numerical vector H is exactly obtained from (11) and (6). Therefore, the error of this aspect could be caused by the process of approximating p0 by Dp. However, it is found that the approximation is also very accurate if the number of nodes is relatively large. Here, we give two examples in Fig. 1 to demonstrate this property. Fig. 1(a) presents the maximum relative error of first-order derivative caused by the differentiation matrix, where the number of nodes is specified to be N + 1 = 3, 4, 5, . . . , 30, respectively. Fig. 1(b) shows the relative error on each node while N = 20. It is observed that the derivatives can be well obtained by the differentiation matrix if the number of nodes is relatively large. 2. For the aspect of solving the system of linear equations: In the present study, we assume that the frequency x is relative large. In effect, if x is very small, the traditional quadrature method can also well deal with this kind of integral, then the advantage of rapid solutions (including the present method) is limited. Under this assumption, the target system of linear Eq. (10) is well-conditioned [16–18] and its solution error is almost negligible. It should be noted that even if the system of linear equations is ill-conditioned (this case corresponds to small x or large N), the impact of its solution error on the final integral result will still be very small if the truncated singular value decomposition (TSVD) method is adopted to solve it [16–18]. Conclusively, the error of the present method should be very small; this makes the method to be attractive.
2393
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396 0
−13
10
10
sinh(x) exp(x)
−5
−14
10
Relative error
Maximum relative error
sinh(x) exp(x)
−10
10
−15
10
10
−15
10
−16
5
10 15 20 25 Total number of nodes
10
30
(a) Maximum relative error for different N
5
10 15 Node No.
20
(b) Relative error on each node while N = 20
Fig. 1. The relative errors of derivative caused by the differentiation matrix.
4. Numerical examples This section focuses on demonstrating the applicability and advantages of the present method, and the algorithm proposed in [12] is taken as a comparison. R1 ixx Example 1. The calculation of Ix ðsinhðbxÞ; sÞ ¼ --1 sinhðbxÞ dx; xs e
Ix ðsinhðbxÞ; sÞ ¼
1 < s < 1 in [9,12]. The analytic result of the integral is
1 ðbþ ixÞs e ½Eiððb þ ixÞð1 sÞÞ Eiððb þ ixÞð1 sÞÞ ip eðbþixÞs ½ Eiððb þ ixÞð1 sÞÞ 2 Eiððb þ ixÞð1 sÞÞ ipg;
where Ei(x) is the exponential integral [13]. Refs [9,12] only consider the imaginary part and the function f(x) = sinh(x) varies very slowly. Here, we set b = 5 and take both the real and imaginary parts into account. We use the present method and the method proposed in [12] (hereinafter denoted by ‘‘Wang method” for convenience) to solve the c.p.v. integral at different number of nodes (N = 2, 4, . . ., , 16). Moreover, different x and s are also involved. For the convenience of comparison, we take the relative error into account, which is defined by Er = jInum Ij/jIj with I and Inum being the real and the numerical integral result, respectively. Fig. 2 shows the relative errors of the two algorithms with ‘‘Li” and ‘‘Wang” denoting the present method and the Wang method, respectively.
0
10
−2
10
−4
r
Relative error (E )
10
−6
10
Wang: τ0.69,ω32 Li: τ0.69,ω32 Wang: τ0.99,ω32 Li: τ0.99,ω32 Wang: τ0.69,ω320 Li: τ0.69,ω320 Wang: τ0.99,ω320 Li: τ0.99,ω320
−8
10
−10
10
−12
10
−14
10
2
4
6
8 10 Number of nodes (N)
12
14
16
Fig. 2. Comparison of relative error between the present method and the Wang method [12] in example 1.
2394
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
From Fig. 2, the following phenomena are observed: 1. Both these two methods become more and more accurate (the curves or the symbols descent) as the number of nodes increases. This is caused by the property that more nodes can give a better interpolation of the unknown functions. 2. The present method is more accurate than the Wang method. For example, the curve ‘—’ is below the symbols ‘’. The same phenomena are observed in the other three groups of curves and symbols (the group of ‘’ and ‘ ’, the group of ‘h’ and ‘’, and the group of ‘4’ and ‘’). The difference in relative error between them is about 2–3 orders of magnitude. The reason of this phenomenon is that the present method does not truncate the series (3) while obtaining the vector H, but the series of Chebyshev expansion in the Wang method is truncated. 3. Both the two methods become more and more accurate as the frequency x increases. This can be seen from the fact that the curves ‘’ and ‘’ are below curves ‘ ’ and ‘—’, and the symbols ‘h’ and ‘4’ are below symbols ‘’ and ‘’. The reason of this phenomenon is that both the two quadrature methods have the asymptotic order of Oðx2 Þ at least. In the preceding example, the function f(x) is very simple. In the following one, we would like to study a c.p.v. integral whose f(x) is more complicated. R 1 f ðxÞ ixx cosð10xÞ ffiffiffiffiffiffiffiffiffiffiffi. Example 2. The calculation of Ix ðf ðxÞ; sÞ ¼ --1 x dx; 1 < s < 1 with f ðxÞ ¼ p se 2 x þ0:1
It is seen that the function f(x) is relatively complicated, so a relatively large number of nodes N (or the degree of Chebysehv polynomial) is required to compute the c.p.v. integral properly. The exact solution of the integral is not known explicitly. In order to make the comparison available, we assume that the result obtained by the Wang method with N = 100 is very close to the exact solution. Under this assumption, the comparison of relative error is presented in Fig. 3, where the number of nodes varies as N = 8, 10, 12, ,. . . , 40. From Fig. 3, the similar phenomena as the preceding example can be observed. There is no need to give more details. In [12], the Wang method is proved to be a good approximation method for the c.p.v. integrals, and some other methods [9,11] have been taken as comparisons to show the advantage of the method. But as shown in Figs. 2 and 3, the method proposed in this paper has even better accuracy than the Wang method. Therefore, the present method is more attractive for the evaluation of c.p.v. integrals. 5. Closing remark A method to calculate the Cauchy principal value (c.p.v.) integral is proposed in this paper. The integral is separated into two oscillatory integrals: a non-singular one and a singular but analytically calculable one. The non-singular one is computed by the improved-Levin method with high accuracy and efficiency, and the singular one is computed by the numerical steepest descent method which serves as an alternative to the analytical solution. Also, the advantage of the proposed method is testified by numerical examples, with the Wang method [12] as a comparison.
0
10
Wang: τ0.69,ω200 Li: τ0.69,ω200 Wang: τ0.99,ω200 Li: τ0.99,ω200 Wang: τ0.69,ω1000 Li: τ0.69,ω1000 Wang: τ0.99,ω1000 Li: τ0.99,ω1000
−2
Relative error (Er)
10
−4
10
−6
10
−8
10
−10
10
10
15
20 25 30 Number of nodes (N)
35
40
Fig. 3. Comparison of relative error between the present method and the Wang method [12] for a more complicated c.p.v. integral.
2395
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Acknowledgment The authors would like to express their gratitude to the editors and the anonymous reviewers for their helpful comments and useful suggestions for improvement of this paper. R 1 ixx Appendix A. An alternate approach of evaluating the integral ws;x ¼ --1 exs dx It is known that
Z Z 1 ixx e eixx ws;x ¼ -dx ¼ limþ dx: e!0 jxsjPe x s 1 x s
ð12Þ
We will show that this integral can be well calculated if some proper transformations are adopted in the complex plane. As shown in Fig. 4, we set a very small half-circle Cr around the pole, i.e., Cr = {s + r eihjh:p ? 0}, where r is the radius of the halfcircle. Consequently, the integral (12) can be rewritten as
Z 1 ixx Z e eixx ws;x ¼ -dx ¼ lim ! ! dx; r!0 1 x s AB þ CD x s
x 2 R:
ð13Þ
In Fig. 4, we have constructed a closed integral path that does not contain the pole; therefore, the following conclusion is immediately obtained according to the Cauchy integral theorem [23]:
Z !
!
!
!
!
AB þC r þCD þ DF þ FE þ EA
eixx dx ¼ 0; xs
x 2 C:
ð14Þ
The combination of (13) and (14) shows that
ws;x ¼ lim r!0
Z !
!
!
C r þAE þ EF þ FD
eixx dx: xs
We now pay the main attention on the computation of the following four integrals: R eixx ! dx. xs
R
eixx C r xs
dx;
R
eixx ! AE xs
dx;
R
eixx ! EF xs
dx, and
FD
1. For the first integral: It is well known that n
eiz 1 z i zn1 1 ¼ þ i þ þ þ , þ /ðzÞ: z 2! z z n!
ð15Þ
In (15), /(z) = i z/2!+ +inzn1/n! + is analytical at z = 0 and j/(0)j = jij = 1, so there exist a sufficient small d to make j/ (z)j 6 2 hold for jzj < d. Applying the coordinate transformation z = x(x s) and (15) to the first integral yields
Z C r
eixx dx ¼ eixs xs
Z C 0r
eiz dz ¼ e z
ix s
Z C 0r
dz þ z
Z
! /ðzÞdz ;
C 0r
where C 0r ¼ fxre ih jh : 0 ! pg. It is known that
Fig. 4. Integral path of the integral.
ð16Þ
2396
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Z C 0r
dz ¼ z
Z p ixreih dh ¼ ip: xreih 0
ð17Þ
And the derivation
Z Z Z j/ðzÞjdz 6 2 dl ¼ 2pxr; /ðzÞdz 6 C 0r C 0r C0 r shows
lim
Z
r!0
C 0r
/ðzÞdz ¼ 0:
ð18Þ
Therefore, the combination of (16)–(18) gives
lim
Z
r!0
C r
eixx dx ¼ ieixs p: xs !
2. For the second integral: Its integral path is AE ¼ f1 þ igjg : 0 ! 1g. In this manner, we have
Z !
AE
eixx dx ¼ ie xs
ix
Z
1
0
exg dg: 1 s þ ig
Obviously, the integral is not oscillatory anymore, and it can be well calculated by the Gauss–Laguerre quadrature method. Moreover, the Gauss–Laguerre quadrature method could be of high efficiency because the abscissas and weights of it can be pre-computed [14,15]. R ixx 3. For the third integral: It is known that jeixxj = 1 and limjxj?11/jx sj = 0, so there should be ! exs dx ¼ 0. EF 4. For the fourth integral: It can obtained in the same way as the second one:
Z !
FD
eixx dx ¼ ieix xs
Z 0
1
exg dg: 1 s þ ig
Up to now, the four integrals are well calculated. This completes the derivation of the alternative approach. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
W.C. Chew, Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, New York, 1990. A. Iserles, S.P. Nørsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proc. Roy. Soc. 461 (2005) 1383–1399. L.N.G. Filon, On a quadrature formula for trigonometric integrals, Proc. Roy. Soc. 49 (1928) 38–47. A. Iserles, S.P. Nørsett, On the computation of highly oscillatory multivariate integrals with critical points, BIT Numer. Math. 46 (2006) 549–566. D. Levin, Procedures for computing one and two dimensional integrals of functions with rapid irregular oscillations, Math. Comput. 38 (158) (1982) 531–538. S. Olver, Moment-free numerical integration of highly oscillatory functions, IMA J. Numer. Anal. 26 (2) (2006) 213–227. D. Huybrechs, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (3) (2006) 1026–1048. G.A. Evans, J.R. Webster, A comparison of some methods for the evaluation of highly oscillatory integrals, J. Comput. Appl. Math. 112 (1999) 55–69. G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comput. 49 (1987) 259–268. G.E. Okecha, Hermite interpolation and a method for evaluating Cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91– 98. M.R. Capobianco, G. Criscuolo, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. H. Wang, S. Xiang, Uniform approximations to Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 215 (2009) 1886–1894. I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, sixth ed., Academic Press, New York, 2004. P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, 1984. M. Kzaz, Convergence acceleration of the Gauss–Laguerre quadrature formula, Appl. Numer. Math. 29 (2) (1999) 201–220. J. Li, X. Wang, T. Wang, A universal solution to one-dimensional highly oscillatory integrals, Sci. China 51 (10) (2008) 1614–1622. J. Li, X. Wang, T. Wang, S. Xiao, An improved Levin quadrature method for highly oscillatory integrals, Appl. Numer. Math. 60 (8) (2010) 833–842. J. Li, X. Wang, T. Wang, C. Shen, Delaminating quadrature method for multi-dimensional highly oscillatory integrals, Appl. Math. Comput. 209 (2) (2009) 327–338. L.N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, 2000. J. Shen, T. Tang, Spectral and High-Order Methods with Applications, Science Press, Beijing, 2006. J.P. Boyd, Chebyshev and Fourier Spectral Methods, DOVER Publications, New York, 2000. C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid dynamics, Springer, New York, 1988. S.G. Krantz, Handbook of Complex Variables, Birkhäuser Press, Boston, 1999.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Evaluation of Cauchy principal value integrals of oscillatory kind Jianbing Li *, Xuesong Wang, Tao Wang College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, PR China
a r t i c l e
i n f o
Keywords: Cauchy principal value integral Oscillatory integral Improved-Levin method Chebyshev differentiation matrix Gauss–Laguerre quadrature
a b s t r a c t This paper studies the evaluation of Cauchy principal value (c.p.v.) integrals of oscillatory R b f ðxÞ ixx kind: --a x s e dx, where s 2 (a, b). The singularity of the c.p.v. integral is transferred to an individual oscillatory integral that can be analytically calculated. The remaining integral which is non-singular can be well calculated with a stable and fast quadrature method, namely the improved-Levin method. Numerical examples confirm the advantages of the proposed method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction The rapid and accurate evaluation of oscillatory integrals is one of the major issues arising in many fields, especially in the computational electromagnetics, optics, and quantum mechanics [1]. Efforts have been made in the past several decades on this issue, and several efficient quadrature methods have been developed. These methods include the asymptotic expansion [2], the Filon(-type) method [3,4], the Levin(-type) method [5,6], the (numerical-)steepest descent method [1,7], etc. [8]. The existing methods are mainly concerned with oscillatory integrals with common integrands. However, in many situations of practical interest, the oscillatory integrals may be singular. For example, the Green function in computational electromagnetics is both oscillatory and singular [1]. Therefore, special attention should be paid on these integrals. In this paper, we study the evaluation of a kind of oscillatory integral with a Cauchy type singularity:
Z Ix ðf ; s; a; bÞ :¼ --
b
a
f ðxÞ ixx e dx; xs
s 2 ða; bÞ; x 0;
ð1Þ
where f(x) is a smooth and non-oscillatory function. Some work has been done on this issue. In [9], Okecha calculates the integral (1) by replacing the function f(x) with the corresponding Lagrange interpolating polynomials on the zeros of the Legendre polynomial and on the singularity point s. The algorithm was improved in [10] by introducing the Hermite interpolation polynomial to obtain a higher efficiency. But this algorithm was proved to be not stable enough by Capobianco and Criscuolo in [11], and the new method they proposed is based on an interpolation procedure at the zeros of the orthogonal polynomials with respect to a Jacobi weight. Recently, Haiyong Wang and Shuhuang Xiang proposed a uniform approximations to the integrals in [12]. It replaces the function f(x) by its Lagrange interpolating polynomial based on the Chebyshev– Gauss–Labatto nodes, and transfers the singularity to an integral that can be analytically calculated. In the method proposed by this paper, the original integral is also separated into two oscillatory integrals: one is non-singular and another one is singular. The non-singular one is computed with an improved-Levin method, and very accurate integral result can be obtained. The second integral has an analytical solution, but several special functions (sine integral and cosine integral) have to be computed. We will propose an alternate approach to calculate the second integral, which does not need to compute special functions and could be more efficient.
* Corresponding author. E-mail address: [email protected] (J. Li). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.07.039
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
2391
2. Quadrature of the Cauchy principal integral based on the improved-Levin method 2.1. Transformation of the Cauchy principal integral The original integral can be easily transformed into a uniform type if the coordinate transformation x0 ¼ ba x þ bþa is em2 2 ployed. Therefore, the present work only takes the uniform type into account:
Z1 f ðxÞ Ix ðf ; sÞ :¼ -e 1 x s
ixx
Z
dx ¼ limþ e!0
jxsjPe
f ðxÞ ixx e dx; xs
s 2 ð1; 1Þ:
The singularity of the integral comes from the existence of a pole at x = s. A reasonable approach of removing the singularity is to make the integrand be the form of ‘0/0’ at the pole, so we transform the original integral into the following two integrals:
Z
Ix ðf ; sÞ ¼
1
1
f ðxÞ f ðsÞ e xs
ixx
dx þ f ðsÞ
Z
1
1
eixx dx , Ix ðf ; sÞ þ f ðsÞws;x : xs
ð2Þ
It is well known that the Taylor expansion of f(x) at x = s is
f ðxÞ ¼ f ðsÞ þ
1 X f ðnÞ ðsÞðx sÞn ; n! n¼1
ð3Þ
so the first integral on the right-hand side of (2) can be re-expressed as
Z
Ix ðf ; sÞ ¼
1
1
f ðxÞ f ðsÞ ixx e dx ¼ xs
Z
1
1
1 X f ðnÞ ðsÞðx sÞn1 ixx e dx: n! n¼1
ð4Þ
Obviously, this integral is non-singular. However, the integral in the second integral is singular at x = s:
Z1 ws;x ¼ --
eixx dx: 1 x s
ð5Þ
Fortunately, this integral has an analytical result [13]:
Re½ws;x ¼ cosðxsÞ½Ciðu1 Þ Ciðju2 jÞ sinðxsÞ½Siðu1 Þ þ Siðju2 jÞ; and
Im½ws;x ¼ sinðxsÞ½Ciðu1 Þ Ciðju2 jÞ þ cosðxsÞ½Siðu1 Þ þ Siðju2 jÞ; where u1 = x(1 s), u2 = x(1 + s), and
CiðuÞ ¼ SiðuÞ ¼
Z Z
u
0
0
u
cosðfÞ 1 df þ ln u þ C; f sin f df; f
u > 0;
are respectively the cosine and sine integrals with C being the Euler constant. In the analytical result, one has to compute two special functions: the cosine and sine integrals. These two integrals can be well handled by some mathematical software such as matlab, but the computational complexity is generally not small enough. In this paper, we present an alternate approach to compute the integral ws,x. By choosing a proper integration path in the complex plane, the integral ws,x is transformed into several singular-free and oscillatory-free integrals which can be well calculated by the Gauss–Laguerre quadrature algorithm [14,15]. Because the abscissas and weights of the Gauss–Laguerre quadrature algorithm can be pre-computed, the complexity of the newly proposed approach is generally very small. For more details, please refer to the Appendix A. Because the second term on the right-hand side of (2), ws,x, can be well calculated, we shall pay the main attention on the computation of the first integral Ix ðf ; sÞ. 2.2. Quadrature of Ix ðf ; sÞ by the improved-Levin method ðsÞ Let hðx; sÞ ¼ f ðxÞf . It is immediately derived from the Taylor expansion (or the L’Hôspital rule for the limitation of ‘0/0’ xs type) that
hðs; sÞ ¼ f 0 ðsÞ: The integral Ix ðf ; sÞ is then denoted by
ð6Þ
2392
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Ix ðf ; sÞ ¼
Z
1
hðx; sÞe
ix x
dx;
ð7Þ
1
where h(x, s) is non-singular. We have proposed an improved-Levin quadrature method for oscillatory integrals [16–18]. The present work adopts this quadrature method to calculate the integral (7). In light of the Levin’s theory, the quadrature of integral (7) is reduced to the calculation of the following ordinary differential equation (ODE) without boundary conditions:
p0 ðxÞ þ ixpðxÞ ¼ hðx; sÞ:
ð8Þ
The Chebyshev differentiation matrix, which is a powerful tool in the spectral method [19–22], is used to solve the ODE. An explicit expression of the matrix follows
Dkj ¼
8 > ck ð1Þkþj > > < cj ðxk xj Þ ;
k; j ¼ 0; . . . ; N; k – j; cj ¼
N > P > > Dkn ; :
2 j ¼ 0; N; 1 j ¼ 1; . . . N 1;
ð9Þ
k ¼ j;
n¼0;n–k
where N + 1 is the number of nodes, xj ¼ cos pNj ; j ¼ 0; 1; . . . ; N is called the Chebyshev–Gauss–Labatto nodes. If the function values of f(x) on the Chebyshev–Gauss–Labatto nodes compose a vector f, then the function values of f0 (x) on these nodes can be approximated by f0 = Df with very high accuracy. Therefore, applying the Chebyshev differentiation matrix on ODE (8) yields
ðD þ ixEÞp ¼ H;
ð10Þ
where E is an identity matrix, p and H are the numerical vectors of p(x) and h(x, s) on the Chebyshev–Gauss–Labatto nodes, respectively, i.e.,
p ¼ ½pðx0 Þ; pðx1 Þ; ; pðxN ÞT ;
H ¼ ½hðx0 ; sÞ; hðx1 ; sÞ; ; hðxN ; sÞT :
The determination of H is a little special. For a given node xj, if xj – s, the component h(xj, s) is computed with
hðxj ; sÞ ¼
f ðxj Þ f ðsÞ ; xj s
ð11Þ
otherwise, the component should be computed according to (6). This treatment does not truncate the Taylor series (3), so it could lead to an exact numerical vector H. It should be also noted that derivatives are not needed in closed form because the differentiation matrix D can be applied to compute the derivatives. The behavior and the proper solution method of (10) have been comprehensively studied in [16–18]. From this system of linear equations, one can easily work out the solution p, and then the integral is yielded as
Ix ðf ; sÞ ¼ pð1Þeix pðN þ 1Þeix : The differentiation matrix D is generally pre-computed, so the determination of the system of linear Eq. (10) is with very small computational complexity. 3. Error analysis The error of calculating the oscillatory integrals via the improved-Levin method may comes from two aspects: the determination of the target system of linear equations and the computation of it. 1. For the aspect of determining the system of linear equations: As has been stated, the numerical vector H is exactly obtained from (11) and (6). Therefore, the error of this aspect could be caused by the process of approximating p0 by Dp. However, it is found that the approximation is also very accurate if the number of nodes is relatively large. Here, we give two examples in Fig. 1 to demonstrate this property. Fig. 1(a) presents the maximum relative error of first-order derivative caused by the differentiation matrix, where the number of nodes is specified to be N + 1 = 3, 4, 5, . . . , 30, respectively. Fig. 1(b) shows the relative error on each node while N = 20. It is observed that the derivatives can be well obtained by the differentiation matrix if the number of nodes is relatively large. 2. For the aspect of solving the system of linear equations: In the present study, we assume that the frequency x is relative large. In effect, if x is very small, the traditional quadrature method can also well deal with this kind of integral, then the advantage of rapid solutions (including the present method) is limited. Under this assumption, the target system of linear Eq. (10) is well-conditioned [16–18] and its solution error is almost negligible. It should be noted that even if the system of linear equations is ill-conditioned (this case corresponds to small x or large N), the impact of its solution error on the final integral result will still be very small if the truncated singular value decomposition (TSVD) method is adopted to solve it [16–18]. Conclusively, the error of the present method should be very small; this makes the method to be attractive.
2393
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396 0
−13
10
10
sinh(x) exp(x)
−5
−14
10
Relative error
Maximum relative error
sinh(x) exp(x)
−10
10
−15
10
10
−15
10
−16
5
10 15 20 25 Total number of nodes
10
30
(a) Maximum relative error for different N
5
10 15 Node No.
20
(b) Relative error on each node while N = 20
Fig. 1. The relative errors of derivative caused by the differentiation matrix.
4. Numerical examples This section focuses on demonstrating the applicability and advantages of the present method, and the algorithm proposed in [12] is taken as a comparison. R1 ixx Example 1. The calculation of Ix ðsinhðbxÞ; sÞ ¼ --1 sinhðbxÞ dx; xs e
Ix ðsinhðbxÞ; sÞ ¼
1 < s < 1 in [9,12]. The analytic result of the integral is
1 ðbþ ixÞs e ½Eiððb þ ixÞð1 sÞÞ Eiððb þ ixÞð1 sÞÞ ip eðbþixÞs ½ Eiððb þ ixÞð1 sÞÞ 2 Eiððb þ ixÞð1 sÞÞ ipg;
where Ei(x) is the exponential integral [13]. Refs [9,12] only consider the imaginary part and the function f(x) = sinh(x) varies very slowly. Here, we set b = 5 and take both the real and imaginary parts into account. We use the present method and the method proposed in [12] (hereinafter denoted by ‘‘Wang method” for convenience) to solve the c.p.v. integral at different number of nodes (N = 2, 4, . . ., , 16). Moreover, different x and s are also involved. For the convenience of comparison, we take the relative error into account, which is defined by Er = jInum Ij/jIj with I and Inum being the real and the numerical integral result, respectively. Fig. 2 shows the relative errors of the two algorithms with ‘‘Li” and ‘‘Wang” denoting the present method and the Wang method, respectively.
0
10
−2
10
−4
r
Relative error (E )
10
−6
10
Wang: τ0.69,ω32 Li: τ0.69,ω32 Wang: τ0.99,ω32 Li: τ0.99,ω32 Wang: τ0.69,ω320 Li: τ0.69,ω320 Wang: τ0.99,ω320 Li: τ0.99,ω320
−8
10
−10
10
−12
10
−14
10
2
4
6
8 10 Number of nodes (N)
12
14
16
Fig. 2. Comparison of relative error between the present method and the Wang method [12] in example 1.
2394
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
From Fig. 2, the following phenomena are observed: 1. Both these two methods become more and more accurate (the curves or the symbols descent) as the number of nodes increases. This is caused by the property that more nodes can give a better interpolation of the unknown functions. 2. The present method is more accurate than the Wang method. For example, the curve ‘—’ is below the symbols ‘’. The same phenomena are observed in the other three groups of curves and symbols (the group of ‘’ and ‘ ’, the group of ‘h’ and ‘’, and the group of ‘4’ and ‘’). The difference in relative error between them is about 2–3 orders of magnitude. The reason of this phenomenon is that the present method does not truncate the series (3) while obtaining the vector H, but the series of Chebyshev expansion in the Wang method is truncated. 3. Both the two methods become more and more accurate as the frequency x increases. This can be seen from the fact that the curves ‘’ and ‘’ are below curves ‘ ’ and ‘—’, and the symbols ‘h’ and ‘4’ are below symbols ‘’ and ‘’. The reason of this phenomenon is that both the two quadrature methods have the asymptotic order of Oðx2 Þ at least. In the preceding example, the function f(x) is very simple. In the following one, we would like to study a c.p.v. integral whose f(x) is more complicated. R 1 f ðxÞ ixx cosð10xÞ ffiffiffiffiffiffiffiffiffiffiffi. Example 2. The calculation of Ix ðf ðxÞ; sÞ ¼ --1 x dx; 1 < s < 1 with f ðxÞ ¼ p se 2 x þ0:1
It is seen that the function f(x) is relatively complicated, so a relatively large number of nodes N (or the degree of Chebysehv polynomial) is required to compute the c.p.v. integral properly. The exact solution of the integral is not known explicitly. In order to make the comparison available, we assume that the result obtained by the Wang method with N = 100 is very close to the exact solution. Under this assumption, the comparison of relative error is presented in Fig. 3, where the number of nodes varies as N = 8, 10, 12, ,. . . , 40. From Fig. 3, the similar phenomena as the preceding example can be observed. There is no need to give more details. In [12], the Wang method is proved to be a good approximation method for the c.p.v. integrals, and some other methods [9,11] have been taken as comparisons to show the advantage of the method. But as shown in Figs. 2 and 3, the method proposed in this paper has even better accuracy than the Wang method. Therefore, the present method is more attractive for the evaluation of c.p.v. integrals. 5. Closing remark A method to calculate the Cauchy principal value (c.p.v.) integral is proposed in this paper. The integral is separated into two oscillatory integrals: a non-singular one and a singular but analytically calculable one. The non-singular one is computed by the improved-Levin method with high accuracy and efficiency, and the singular one is computed by the numerical steepest descent method which serves as an alternative to the analytical solution. Also, the advantage of the proposed method is testified by numerical examples, with the Wang method [12] as a comparison.
0
10
Wang: τ0.69,ω200 Li: τ0.69,ω200 Wang: τ0.99,ω200 Li: τ0.99,ω200 Wang: τ0.69,ω1000 Li: τ0.69,ω1000 Wang: τ0.99,ω1000 Li: τ0.99,ω1000
−2
Relative error (Er)
10
−4
10
−6
10
−8
10
−10
10
10
15
20 25 30 Number of nodes (N)
35
40
Fig. 3. Comparison of relative error between the present method and the Wang method [12] for a more complicated c.p.v. integral.
2395
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Acknowledgment The authors would like to express their gratitude to the editors and the anonymous reviewers for their helpful comments and useful suggestions for improvement of this paper. R 1 ixx Appendix A. An alternate approach of evaluating the integral ws;x ¼ --1 exs dx It is known that
Z Z 1 ixx e eixx ws;x ¼ -dx ¼ limþ dx: e!0 jxsjPe x s 1 x s
ð12Þ
We will show that this integral can be well calculated if some proper transformations are adopted in the complex plane. As shown in Fig. 4, we set a very small half-circle Cr around the pole, i.e., Cr = {s + r eihjh:p ? 0}, where r is the radius of the halfcircle. Consequently, the integral (12) can be rewritten as
Z 1 ixx Z e eixx ws;x ¼ -dx ¼ lim ! ! dx; r!0 1 x s AB þ CD x s
x 2 R:
ð13Þ
In Fig. 4, we have constructed a closed integral path that does not contain the pole; therefore, the following conclusion is immediately obtained according to the Cauchy integral theorem [23]:
Z !
!
!
!
!
AB þC r þCD þ DF þ FE þ EA
eixx dx ¼ 0; xs
x 2 C:
ð14Þ
The combination of (13) and (14) shows that
ws;x ¼ lim r!0
Z !
!
!
C r þAE þ EF þ FD
eixx dx: xs
We now pay the main attention on the computation of the following four integrals: R eixx ! dx. xs
R
eixx C r xs
dx;
R
eixx ! AE xs
dx;
R
eixx ! EF xs
dx, and
FD
1. For the first integral: It is well known that n
eiz 1 z i zn1 1 ¼ þ i þ þ þ , þ /ðzÞ: z 2! z z n!
ð15Þ
In (15), /(z) = i z/2!+ +inzn1/n! + is analytical at z = 0 and j/(0)j = jij = 1, so there exist a sufficient small d to make j/ (z)j 6 2 hold for jzj < d. Applying the coordinate transformation z = x(x s) and (15) to the first integral yields
Z C r
eixx dx ¼ eixs xs
Z C 0r
eiz dz ¼ e z
ix s
Z C 0r
dz þ z
Z
! /ðzÞdz ;
C 0r
where C 0r ¼ fxre ih jh : 0 ! pg. It is known that
Fig. 4. Integral path of the integral.
ð16Þ
2396
J. Li et al. / Applied Mathematics and Computation 217 (2010) 2390–2396
Z C 0r
dz ¼ z
Z p ixreih dh ¼ ip: xreih 0
ð17Þ
And the derivation
Z Z Z j/ðzÞjdz 6 2 dl ¼ 2pxr; /ðzÞdz 6 C 0r C 0r C0 r shows
lim
Z
r!0
C 0r
/ðzÞdz ¼ 0:
ð18Þ
Therefore, the combination of (16)–(18) gives
lim
Z
r!0
C r
eixx dx ¼ ieixs p: xs !
2. For the second integral: Its integral path is AE ¼ f1 þ igjg : 0 ! 1g. In this manner, we have
Z !
AE
eixx dx ¼ ie xs
ix
Z
1
0
exg dg: 1 s þ ig
Obviously, the integral is not oscillatory anymore, and it can be well calculated by the Gauss–Laguerre quadrature method. Moreover, the Gauss–Laguerre quadrature method could be of high efficiency because the abscissas and weights of it can be pre-computed [14,15]. R ixx 3. For the third integral: It is known that jeixxj = 1 and limjxj?11/jx sj = 0, so there should be ! exs dx ¼ 0. EF 4. For the fourth integral: It can obtained in the same way as the second one:
Z !
FD
eixx dx ¼ ieix xs
Z 0
1
exg dg: 1 s þ ig
Up to now, the four integrals are well calculated. This completes the derivation of the alternative approach. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
W.C. Chew, Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, New York, 1990. A. Iserles, S.P. Nørsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proc. Roy. Soc. 461 (2005) 1383–1399. L.N.G. Filon, On a quadrature formula for trigonometric integrals, Proc. Roy. Soc. 49 (1928) 38–47. A. Iserles, S.P. Nørsett, On the computation of highly oscillatory multivariate integrals with critical points, BIT Numer. Math. 46 (2006) 549–566. D. Levin, Procedures for computing one and two dimensional integrals of functions with rapid irregular oscillations, Math. Comput. 38 (158) (1982) 531–538. S. Olver, Moment-free numerical integration of highly oscillatory functions, IMA J. Numer. Anal. 26 (2) (2006) 213–227. D. Huybrechs, On the evaluation of highly oscillatory integrals by analytic continuation, SIAM J. Numer. Anal. 44 (3) (2006) 1026–1048. G.A. Evans, J.R. Webster, A comparison of some methods for the evaluation of highly oscillatory integrals, J. Comput. Appl. Math. 112 (1999) 55–69. G.E. Okecha, Quadrature formulae for Cauchy principal value integrals of oscillatory kind, Math. Comput. 49 (1987) 259–268. G.E. Okecha, Hermite interpolation and a method for evaluating Cauchy principal value integrals of oscillatory kind, Kragujevac J. Math. 29 (2006) 91– 98. M.R. Capobianco, G. Criscuolo, On quadrature for Cauchy principal value integrals of oscillatory functions, J. Comput. Appl. Math. 156 (2003) 471–486. H. Wang, S. Xiang, Uniform approximations to Cauchy principal value integrals of oscillatory functions, Appl. Math. Comput. 215 (2009) 1886–1894. I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, sixth ed., Academic Press, New York, 2004. P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, 1984. M. Kzaz, Convergence acceleration of the Gauss–Laguerre quadrature formula, Appl. Numer. Math. 29 (2) (1999) 201–220. J. Li, X. Wang, T. Wang, A universal solution to one-dimensional highly oscillatory integrals, Sci. China 51 (10) (2008) 1614–1622. J. Li, X. Wang, T. Wang, S. Xiao, An improved Levin quadrature method for highly oscillatory integrals, Appl. Numer. Math. 60 (8) (2010) 833–842. J. Li, X. Wang, T. Wang, C. Shen, Delaminating quadrature method for multi-dimensional highly oscillatory integrals, Appl. Math. Comput. 209 (2) (2009) 327–338. L.N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, 2000. J. Shen, T. Tang, Spectral and High-Order Methods with Applications, Science Press, Beijing, 2006. J.P. Boyd, Chebyshev and Fourier Spectral Methods, DOVER Publications, New York, 2000. C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid dynamics, Springer, New York, 1988. S.G. Krantz, Handbook of Complex Variables, Birkhäuser Press, Boston, 1999.