Fast Computation of Singular Oscillatory Fourier Transforms

Abstract

We consider the problem of the numerical evaluation of singular oscillatory Fourier transforms ∫^𝑏_𝑎(𝑥−𝑎)^𝛼(𝑏−𝑥)^𝛽𝑓(𝑥)𝑒^𝑖𝜔𝑥𝑑𝑥, where 𝛼>−1and𝛽>−1. Based on substituting the original interval of integration by the paths of steepest descent, if 𝑓 is analytic in the complex region 𝐺 containing [𝑎, 𝑏], the computation of integrals can be transformed into the problems of integrating two integrals on [0, ∞) with the integrand that does not oscillate and decays exponentially fast, which can be efficiently computed by using the generalized Gauss Laguerre quadrature rule. The efficiency and the validity of the method are demonstrated by both numerical experiments and theoretical results. More importantly, the presented method in this paper is also a great improvement of a Filon-type method and a Clenshaw-Curtis-Filon-type method shown in Kang and Xiang (2011) and the Chebyshev expansions method proposed in Kang et al. (2013), for computing the above integrals.

1. Introduction

In many areas of science and engineering one often encounters the problem of computing rapidly oscillatory integrals due to their frequent occurrences in wide fields ranging from quantum chemistry, image analysis, electrodynamics, and computerized tomography to fluid mechanics. The numerical evaluation can be difficult when the parameter 𝜔 is large, because in that case the integrand is highly oscillatory. A prohibitively large number of quadrature points are needed if one uses a classic rule such as Gaussian quadrature or any quadrature method based on (piecewise) polynomial interpolation of the integrand. In most of the cases, such integrals cannot be calculated analytically and one has to resort to numerical methods. In the past nearly hundred years ago, great many methods have been developed for generalized Fourier transformation ∫^𝑏_𝑎𝑓(𝑥)𝑒^{𝑖𝜔𝑔(𝑥)}𝑑𝑥, such as the Filon [1–4], Clenshaw-Curtis-type [5, 6], Filon-type [7, 8], asymptotic [7], Levin [9, 10], generalized quadrature rule [11, 12], Levin-type [13], and complex integration methods [14–18].

In the present paper, based on special contours and the generalized Gauss Laguerre quadrature rule, we will be concerned with the computation for oscillatory Fourier transform of the form 𝐼𝑓=^𝑏_𝑎(𝑥−𝑎)^𝛼(𝑏−𝑥)^𝛽𝑓(𝑥)𝑒^𝑖𝜔𝑥𝑑𝑥, (1) where 𝑓(𝑥) is a sufficiently smooth function in [𝑎,𝑏], 𝜔 is a large parameter, 𝑎 and 𝑏 are real and finite, and 𝛼>−1, 𝛽>−1. In (1), if 𝛼≥0, −1<𝛽<0, or 𝛽≥0, −1<𝛼<0, the integrand has a singularity of a simple type at one sided endpoint 𝑏 or 𝑎; if −1<𝛼<0, −1<𝛽<0, the integrand has singularities of a complicated type at both endpoints of the interval; if, in particular, 𝛼,𝛽≥0, the integrand has no singularity and zero points or stationary points. When the integrand containing algebraic singularities becomes highly oscillatory, it presents more serious difficulties in obtaining numerical convergence of the integration than the computation of (1) whose integrand does not involve singularity. Such integrals (1) with the weak singularities are applied to the numerical approximations of solutions to Volterra integral equations of the first kind [19, 20]. In addition, it is well-known that the Radon transform, which is an important role in the CT, PET, and SPECT technology of medical sciences and widely applicable to tomography, the creation of an image from the scattering data associated to cross-sectional scans of an object, is also closely related to this form of oscillatory singular integrals [21, 22]. These singularities are also called singularities of the Radon transform in medical tomography. Further, the numerical integration of such integrals (1) is used to the solution of the singular integral equation for classical crack problems in plane and antiplane elasticity [23]. Moreover, they can be taken as model integrals of those appearing in solving integral equations, such as, in high-frequency acoustic scattering, for example, high-frequency Helmholtz equation in two dimensions, where those kernels have algebraic or logarithmic singularities ([24, 25] and references therein), which is also our main target application. Because of such a wide range of applications, it is of great importance for the study of the numerical integration of such integrals. In fact, the integral (1) has gained popularity in literatures. In 1955, Erdélyi [26] established asymptotic expansions using neutralizer functions by van der Corput and general integration by parts. One year later, this sort of asymptotic expansions was listed in the treatise [15] by Erdélyi. Three years later, a similar asymptotic result was reestablished by Lighthill in [27], by generalized function theory. In 1971, in the case where 𝑓(𝑥) is analytic in a region containing [𝑎,𝑏], a straightforward proof based on contour integration was published by Lyness in [28]. In 2008 and 2009, Lyness [29, 30] presented asymptotic expansions by theory involving inverse functions. For details one can refer to [29, 30]. In 2011, the authors of   [31] presented a Filon-type method and a Clenshaw-Curtis-Filon-type method for computing the integral (1), the error of which satisfies 𝑂(𝜔^{−𝑠−𝑟−2}), where 𝑠 is the highest multiplicity of the Hermite interpolation at the endpoints and 𝑟=min{𝛼,𝛽}. In 2013, the recent literature [32] gave a widely used Chebyshev expansions method depending on the frequency 𝜔 for computing many types of singular oscillatory integrals, one of which is such integral (1). Based on these relevant background literatures above, in Section 2 of this paper, thanks to analytic continuation, special contours, and generalized Gauss-Laguerre quadrature rule, we devise efficient method to compute the class of integrals (1). Its asymptotic order, 𝑂((𝑛!Γ(𝑛+𝑟₂+1)/(2𝑛)!)(1/𝜔^{2𝑛+𝑟1+1})), 𝑟₁=min{𝛼,𝛽}, and 𝑟₂=max{𝛼,𝛽}, is nearly twice as high as that of the Filon-type method and the Clenshaw-Curtis-Filon-type method [31, 33] outlined above, while using the same number of function evaluations under the given conditions. The results differ from those in previous research in the sense that the constructed rules are asymptotically optimal; that is, among all known methods for oscillatory integrals they deliver the highest possible asymptotic order of convergence, relative to the required number of evaluations of the integrand. It can combine a fixed computational cost and very high asymptotic order with numerical convergence. Of course, the presented method is also much more efficient than the Chebyshev expansions method proposed in [32].

An outline of this paper is as illustrated below. In the next section, we will illustrate some of the main ideas with the theoretical analysis by choosing a special path and demonstrate a decomposition of the integral (1) in details. Meanwhile, the construction and error of a quadrature rule will be also established by the generalized Gauss-Laguerre quadrature rule. In Section 3, some vigorous and robust numerical results will show the accuracy and efficiency of the proposed approach.

2. Main Results

2.1. A Decomposition of the Integral (1)

In this subsection, we will study the problem of computing the integral (1) in depth. Here, we provide the illustration of the special integration paths for the integrand (𝑥−𝑎)^𝛼(𝑏−𝑥)^𝛽𝑓(𝑥)𝑒^𝑖𝜔𝑥 (see Figure 1). Here, let 𝐺, 𝐺₁, and 𝐺₂ denote the regions 𝐺={𝑧∈ℂ∣𝑎≤ℜ(𝑧)≤𝑏,0≤ℑ(𝑧)≤𝐿}, 𝐺₁={|𝑧−𝑎|≤𝜀,0≤arg(𝑧)≤𝜋/2}, and 𝐺₂={|𝑧−𝑏|≤𝜀,𝜋/2≤arg(𝑧)≤𝜋}, respectively, where 𝐿 is a large number and the positive 𝜀 is small enough such that 𝐺 contains 𝐺₁ and 𝐺₂.

Figure 1 

Illustration of the integration paths for the integrand (𝑥−𝑎)𝛼(𝑏−𝑥)𝛽𝑓(𝑥)𝑒𝑖𝜔𝑥.

Theorem 1. Suppose that 𝑓 is an analytic function in the region 𝐺={𝑧∈ℂ∣𝑎≤ℜ(𝑧)≤𝑏,ℑ(𝑧)≥0}; then 𝐼𝑓=^𝑏_𝑎(𝑥−𝑎)^𝛼(𝑏−𝑥)^𝛽𝑓(𝑥)𝑒^𝑖𝜔𝑥=𝑖𝑑𝑥_𝛼+1𝑒^𝑖𝜔𝑎𝜔^𝛼+1∞₀𝑝𝑏−𝑎−𝑖𝜔^𝛽𝑓𝑝𝑎+𝑖𝜔𝑝^𝛼𝑒^−𝑝−()𝑑𝑝−1^𝛽+1𝑖^𝛽+1𝑒^𝑖𝜔𝑏𝜔^𝛽+1×^∞₀𝑝𝑏−𝑎+𝑖𝜔^𝛼𝑓𝑝𝑏+𝑖𝜔𝑝^𝛽𝑒^−𝑝𝑑𝑝. (2)

Proof. Since 𝑓(𝑧) is analytic within the complex region 𝐺, then the integrand (𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧 is also analytic (has no singularity) in the region 𝐺 except 𝐺₁ and 𝐺₂ (𝑧∈𝐺⧵(𝐺₁+𝐺₂)), enclosed by Γ_𝑗, 𝑗=1,2,…,6, defined by the following parametric forms Γ₁∶ℎ₁()=𝑥, 𝑎+𝜀≤𝑥≤𝑏−𝜀, Γ₂∶ℎ₂()=𝑏+𝜀𝑒^𝑖𝜃1, Γ₃∶ℎ₃()=𝑏+𝑖, Γ₄∶ℎ₄()=𝑥+𝑖𝐿, 𝑥∈[𝑎,𝑏], Γ₅∶ℎ₅()=𝑎+𝑖, and Γ₆∶ℎ₆()=𝑎+𝜀𝑒^𝑖𝜃2, where ∈[𝜀,𝐿], 𝜃₁∈[𝜋/2,𝜋], 𝜃₂∈[0,𝜋/2], and an arbitrary 𝜀>0. Moreover, (𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧 is continuous on all contours containing Γ_𝑗, 𝑗=1,2,…,6 (see Figure 1). By the Cauchy theorem [34], we obtain _{Γ1+Γ2+Γ3+Γ4+Γ5+Γ6}(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=0, (3) where all the contours choose the counterclockwise direction as positive direction. In the sequel, we parameterize each of these integrals in a specific way.
First, there exists 𝑀≥0, such that 𝑓(𝑧)≤𝑀,𝑧∈𝐺. (4) Since 0≤_Γ2(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧≤𝜀^𝛽+1𝜋_𝜋/2𝑏−𝑎+𝜀𝑒^{𝑖𝜃1𝛼}𝑓𝑏+𝜀𝑒^𝑖𝜃1𝑑𝜃₁≤𝜀^𝛽+1𝑀^𝜋_𝜋/2𝑏−𝑎+𝜀𝑒^{𝑖𝜃1𝛼}𝑑𝜃₁⟶0,as𝜀⟶0, (5) then lim_𝜀→0Γ2(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=0. (6)
Similarly, it is easy to get lim_𝜀→0Γ6(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=0. (7)
Since the integrand (𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧) is analytic in the region 𝐺 except 𝐺₁ and 𝐺₂, then (𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧 is also analytic in the region 𝐺 except 𝐺₁ and 𝐺₂. So there must exist a nonnegative number 𝑀 such that (𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)≤𝑀,𝑧∈𝐺⧵𝐺₁+𝐺₂. (8) Therefore, 0≤_Γ4(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧≤𝑒^{−𝐿𝜔𝑏}_𝑎(𝑥−𝑎+𝑖𝐿)^𝛼(𝑏−𝑥−𝑖𝐿)^𝛽𝑓(≤()𝑥+𝑖𝐿)𝑑𝑥𝑏−𝑎𝑀𝑒^−𝐿𝜔⟶0,as𝐿⟶∞; (9) namely, lim_{𝐿→∞Γ4}(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=0. (10)
Since the Gamma function Γ(𝑧) and the incomplete Gamma function Γ(𝑠,𝑧) [35] are defined by Γ(𝑧)=^∞₀𝑡^𝑧−1𝑒^−𝑡𝑑𝑡,Γ(𝑠,𝑧)=^∞_𝑠𝑡^𝑧−1𝑒^−𝑡𝑑𝑡, (11) which hold for the equality Γ(𝛼+1)=Γ(𝛼+1,0), (12) then ^𝜀₀^𝛼𝑒^−𝜔=1𝑑𝜔^𝛼+1₀^𝜔𝜀𝑡^𝛼𝑒^−𝑡=1𝑑𝑡𝜔^𝛼+1∞₀𝑡^𝛼+1−1𝑒^−𝑡𝑑𝑡−^∞_𝜔𝜀𝑡^𝛼+1−1𝑒^−𝑡=1𝑑𝑡𝜔^𝛼+1[Γ(𝛼+1)−Γ(𝛼+1,𝜔𝜀)]⟶0,as𝜀⟶0. (13) Therefore, 0≤^𝜀₀(𝑏−𝑎−𝑖)^𝛽𝛼𝑓(𝑎+𝑖)𝑒^−𝜔𝑑≤𝑀⋅𝑀₁^𝜀₀^𝛼𝑒^−𝜔𝑑⟶0,as𝜀⟶0, (14) where there must exist a nonnegative 𝑀₁ such that (𝑏−𝑎−𝑖)^𝛽≤𝑀₁,∈[0,𝜀]. (15)
Thus, using (14), lim_𝜀→0Γ5(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=lim_𝜀→0−𝑖^𝛼+1𝑒^𝑖𝜔𝑎×^𝐿_𝜀(𝑏−𝑎−𝑖)^𝛽𝛼𝑓(𝑎+𝑖)𝑒^−𝜔𝑑=lim_𝜀→0−𝑖^𝛼+1𝑒^𝑖𝜔𝑎×^𝐿₀(𝑏−𝑎−𝑖)^𝛽𝛼𝑓(𝑎+𝑖)𝑒^−𝜔−_𝜀0𝑑(𝑏−𝑎−𝑖)^𝛽𝛼𝑓(𝑎+𝑖)𝑒^−𝜔𝑑=−𝑖^𝛼+1𝑒^{𝑖𝜔𝑎𝐿}₀(𝑏−𝑎−𝑖)^𝛽𝛼𝑓(𝑎+𝑖)𝑒^−𝜔𝑖𝑑=−^𝛼+1𝑒^𝑖𝜔𝑎𝜔^𝛼+1×₀^𝜔𝐿𝑝𝑏−𝑎−𝑖𝜔^𝛼𝑝^𝛼𝑓𝑝𝑎+𝑖𝜔𝑒^−𝑝𝑑𝑝. (16)
In the same way, we have lim_𝜀→0Γ3(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧=()𝑑𝑧−1^𝛽+1𝑖^𝛽+1𝑒^𝑖𝜔𝑏𝜔^𝛽+1×₀^𝜔𝐿𝑝𝑏−𝑎+𝑖𝜔^𝛼𝑓𝑝𝑏+𝑖𝜔𝑝^𝛽𝑒^−𝑝𝑑𝑝. (17)
Therefore, combining (3), (6), (9), (10), (16), and (17), we gain ^𝑏_𝑎(𝑥−𝑎)^𝛼(𝑏−𝑥)^𝛽𝑓(𝑥)𝑒^𝑖𝜔𝑥𝑑𝑥=lim_𝜀→0Γ1(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽𝑓(𝑧)𝑒^𝑖𝜔𝑧𝑑𝑧=−lim_𝜀→0lim_{𝐿→∞Γ2+Γ3+Γ4+Γ5+Γ6}(𝑧−𝑎)^𝛼(𝑏−𝑧)^𝛽×𝑓(𝑧)𝑒^𝑖𝜔𝑧=𝑖𝑑𝑧_𝛼+1𝑒^𝑖𝜔𝑎𝜔^𝛼+1∞₀𝑝𝑏−𝑎−𝑖𝜔^𝛽𝑓𝑝𝑎+𝑖𝜔𝑝^𝛼𝑒^−𝑝−()𝑑𝑝−1^𝛽+1𝑖^𝛽+1𝑒^𝑖𝜔𝑏𝜔^𝛽+1×^∞₀𝑝𝑏−𝑎+𝑖𝜔^𝛼𝑓𝑝𝑏+𝑖𝜔𝑝^𝛽𝑒^−𝑝𝑑𝑝. (18)
This completes the proof.

The succeeding part is to consider the efficient evaluation of the nethermost formulas of (2).

2.2. Calculation of 𝐼[𝑓] by the Generalized Gauss Laguerre Quadrature Rule

The key point in the interpolatory rule for the infinite interval [0,∞) is, of course, the rule of Gauss-type: ^∞₀𝑤(𝑥)𝑆(𝑥)𝑑𝑥≈^𝑛_𝑘=1𝑤_𝑘𝑆𝑥_𝑘, (19) where the 𝑥_𝑘 and 𝑤_𝑘 have been determined so that the formula is exact for a polynomial 𝑆(𝑥) up to degree 2𝑛−1. But the most widely employed efficient approach for the infinite integrals is the Gauss-Laguerre quadrature rule or generalized Gauss Laguerre quadrature rule.

Generalized Laguerre polynomials are orthogonal with respect to the more general weight function 𝑥^𝜆𝑒^−𝑥, 𝜆>−1. From [1], the generalized Gauss Laguerre quadrature rule is given, as follows: ^∞₀𝑆(𝑥)𝑥^𝜆𝑒^−𝑥𝑑𝑥=^𝑛_𝑘=1𝑤_𝑘^(𝜆)𝑆𝑥_𝑘^(𝜆)+𝑛!Γ(𝑛+𝜆+1)()!𝑆2𝑛^(2𝑛)(𝜉),0<𝜉<∞, (20) where the nodes 𝑥_𝑘^(𝜆) are the zeros of the generalized or associated Laguerre polynomial 𝐿_𝑛^(𝜆)(𝑥) [1] and 𝑤_𝑘^(𝜆)=Γ(𝑛+𝜆+1)𝑥_𝑘^(𝜆)𝑛!𝐿^(𝜆)_𝑛+1𝑥_𝑘^(𝜆)2. (21)

In (2), the integrands have a singularity of the form 𝑝^𝛼 or 𝑝^𝛽, as 𝑝→0. Fortunately, this type of singularity can be handled efficiently by the generalized Gauss-Laguerre quadrature rule.

2.2.1. Calculation of 𝐼[𝑓] by the Generalized Gauss Laguerre Quadrature Rule

The generalized Gauss-Laguerre formula (20) with 𝑛 points 𝑥_𝑘^(𝜆) and weights 𝑤_𝑘^(𝜆) is applied to (2). Then we can easily obtain the following 𝑛-point approximation formula to (1): 𝑄_𝑛𝑓=𝑖^𝛼+1𝑒^𝑖𝜔𝑎𝜔^𝛼+1×^𝑛_𝑘=1𝑤_𝑘^(𝛼)𝑥𝑏−𝑎−𝑖^𝑘(𝛼)𝜔^𝛽𝑓𝑥𝑎+𝑖^𝑘(𝛼)𝜔−(−1)^𝛽+1𝑖^𝛽+1𝑒^𝑖𝜔𝑏𝜔^𝛽+1×^𝑛_𝑘=1𝑤_𝑘^(𝛽)𝑥𝑏−𝑎+𝑖^𝑘(𝛽)𝜔^𝛼𝑓𝑥𝑏+𝑖^𝑘(𝛽)𝜔, (22) where 𝑥_𝑘^(𝛼) and 𝑥_𝑘^(𝛽) are the zeros of the generalized or associated Laguerre polynomials 𝐿_𝑛^(𝛼)(𝑥) and 𝐿_𝑛^(𝛽)(𝑥), respectively. And, 𝑤_𝑘^(𝛼) and 𝑤_𝑘^(𝛽) are given by choosing “𝜆=𝛼” and “𝜆=𝛽” in (21), respectively.

Theorem 2. If 𝑓 holds for the conditions of Theorem 1, then we have 𝐼𝑓−𝑄_𝑛𝑓∼𝑂𝑛!Γ𝑛+𝑟₂+1()!12𝑛𝜔^{2𝑛+𝑟1+1},as𝜔⟶∞, (23) where 𝑟₁=min{𝛼,𝛽} and 𝑟₂=max{𝛼,𝛽}.

Proof. For 𝑟₁=min{𝛼,𝛽} and 𝑟₂=max{𝛼,𝛽}, since 𝜔^{2𝑛+𝑟1+1}=𝜔^2𝑛+𝛼+1,𝜔^{2𝑛+𝑟1+1}<𝜔^2𝑛+𝛽+1,Γ𝑛+𝑟₂,Γ+1=Γ𝑛+𝛽+1𝑛+𝑟₂+1>Γ(𝑛+𝛼+1),𝛼≤𝛽; (24) or 𝜔^{2𝑛+𝑟1+1}=𝜔^2𝑛+𝛽+1,𝜔^{2𝑛+𝑟1+1}<𝜔^2𝑛+𝛼+1,Γ𝑛+𝑟₂(),Γ+1=Γ𝑛+𝛼+1𝑛+𝑟₂+1>Γ𝑛+𝛽+1,𝛼>𝛽 (25) together with (2), (20), and (22), we have 𝐼𝑓−𝑄_𝑛𝑓=𝑛!Γ(𝑛+𝜆+1)()!𝑆2𝑛^(2𝑛)(𝜉)=𝑖^𝛼+1𝑒^𝑖𝜔𝑎𝜔^𝛼+1𝑛!Γ(𝑛+𝛼+1)()!×𝑑2𝑛𝑏−𝑎−𝑖𝑞/𝜔^𝛽𝑓𝑎+𝑖𝑞/𝜔𝑑𝑞^2𝑛_𝑞=𝜉−(−1)^𝛽+1𝑖^𝛽+1𝑒^𝑖𝜔𝑏𝜔^𝛽+1𝑛!Γ𝑛+𝛽+1()!×𝑑2𝑛𝑏−𝑎+𝑖𝑞/𝜔^𝛼𝑓𝑏+𝑖𝑞/𝜔𝑑𝑞^2𝑛_𝑞=𝜉()∼𝑂𝑛!Γ𝑛+𝛼+1()!𝜔2𝑛^{−2𝑛−𝛼−1}−𝑂𝑛!Γ𝑛+𝛽+1()!𝜔2𝑛^{−2𝑛−𝛽−1}∼𝑂𝑛!Γ𝑛+𝑟₂+1()!12𝑛𝜔^{2𝑛+𝑟1+1},as𝜔⟶∞. (26) This completes this proof.

From the error formula (23) it can be seen that more accurate approximations can be obtained for the case of the fixed number of nodes and increasing frequency 𝜔, and for the case of the fixed frequency 𝜔 and increasing the number of nodes. In this case, the integral (1) depends asymptotically on the behavior of the integrand near the endpoints 𝑎 and 𝑏.

In the subsequent section, we will make use of numerical examples to illustrate the absolute error and approximate values.

3. Numerical Examples

In this section, we present the results of numerical experiments obtained by using the proposed method. Our algorithm is compared with the Clenshaw-Curtis-Filon-type method presented in [31]. The nodes 𝑥_𝑘^(𝜆) and weights 𝑤_𝑘^(𝜆) of the generalized Gauss-Laguerre quadrature rule are given in [35]. In order to conduct the experiments, we require the knowledge on the exact values of the integrals of the form (1). The values that we assume to be accurate are computed in Maple 13.02 using the 32 decimal digits precision arithmetic. The compared algorithms (the Clenshaw-Curtis-Filon-type method and the presented one) are implemented in Matlab 7.0.1, taking advantage of its fast vectorised arithmetic operations. The experiments are performed on the desktop computer with AMD Athlon(tm) 64 X2 Dual Core Processor 4000 + (2100 Mhz) and 1 GB of RAM.

Tables 1, 2, 3, 4, 5, and 6 demonstrate the absolute errors and approximate values in 𝑛-point approximations by the proposed method to the above integrals. Furthermore, they exhibit the fast convergence of the approximations as 𝑛 increases. Meanwhile, all these tables above show that more and more accurate approximations can be obtained as 𝜔 increases and 𝑛 is fixed. Conversely, as 𝑛 increases and 𝜔 is fixed, higher accuracy can be also achieved. Moreover, they exhibit that the approach requires very small number of function evaluations to produce approximations in higher accuracy, where we only choose several nodes, 𝑛=2,3,4, and so forth. Also, for sufficiently large 𝜔 and some 𝑛, it is very easy to reach machine precision in Matlab. For the case of low or moderate frequency 𝜔, by adding the number of the node points 𝑛, we can get very accurate approximations. Moreover, Tables 1(a), 1(b), 2(a), and 2(b) show the efficiency and accuracy of the proposed method, compared to the Clenshaw-Curtis-Filon-type method [31]. In [32], the Chebyshev expansions method depends on 𝜔; for example, 𝑁₁(𝜔)=⌊1.05𝜔+57.5⌋, (27) where 𝑁₁(𝜔) is the number of the required truncated terms. For the case of moderate and high frequency 𝜔, the Chebyshev expansions method requires longer time consuming in attaining high precision [32]. For details one can refer to [32]. Furthermore, both numerical analysis in Theorem 2 and numerical experiments above show that the presented method in this paper shares the advantageous property that its accuracy improves greatly when 𝜔 increases. Therefore, it is obvious that the presented method in this paper is much more efficient than the Chebyshev expansions method for computing the integral (1).

Table 1 

(a) Absolute errors in 𝑛-point approximations by the proposed method to the integral ∫¹₀𝑥^−(1/2)(𝑥−1)^−(1/2)cos𝑥𝑒^𝑖𝜔𝑥𝑑𝑥. (b) Absolute errors in 𝑛-point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral ∫¹₀𝑥^−(1/2)(𝑥−1)^−(1/2)cos𝑥𝑒^𝑖𝜔𝑥𝑑𝑥.

Table 2 

(a) Approximate values in 2-point approximations by the proposed method to the integral ∫¹₀𝑥^−(1/2)(𝑥−1)^−(1/2)cos𝑥𝑒^𝑖𝜔𝑥𝑑𝑥, compared with exact values. (b) Approximate values in 2-point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral ∫¹₀𝑥^−(1/2)(𝑥−1)^−(1/2)cos𝑥𝑒^𝑖𝜔𝑥𝑑𝑥, compared with exact values.

Table 3 

Absolute errors in 𝑛-point approximations by the proposed method to the integral ∫³₂(𝑥−2)^−(1/2)(3−𝑥)^−(1/2)sin𝑥𝑒^𝑖𝜔𝑥𝑑𝑥.

Table 4 

Approximate values in 2-point approximations by the proposed method to the integral ∫³₂(𝑥−2)^−(1/2)(3−𝑥)^−(1/2)sin𝑥𝑒^𝑖𝜔𝑥𝑑𝑥, compared with exact values.

Table 5 

Absolute errors in 𝑛-point approximations by the proposed method to the integral ∫¹₋₁(𝑥+1)^−(1/2)(1−𝑥)^−(1/2)𝑒^𝑥𝑒^𝑖𝜔𝑥𝑑𝑥.

Table 6 

Approximate values in 2-point approximations by the proposed method to the integral ∫¹₋₁(𝑥+1)^−(1/2)(1−𝑥)^−(1/2)𝑒^𝑥𝑒^𝑖𝜔𝑥𝑑𝑥, compared with exact values.

4. Concluding Remark

In this paper, for an analytic and not increasing fast 𝑓 in 𝐺={𝑧∈ℂ∣𝑎≤ℜ(𝑧)≤𝑏,ℑ(𝑧)≥0}, we present an efficient method for handling the integral (1) based on analytical continuation and special contours. For the case of moderate, high, and even low frequency 𝜔, this approach is easy to implement by using the generalized Gauss-Laguerre quadrature rule without using high precision arithmetic and gives very accurate results. Moreover, the proposed method requires shorter time consuming in high precision in Matlab and its convergence is much faster than that of those methods presented in [31, 32]. In other words, it can combine lower computational cost and higher asymptotic order with numerical convergence than that of the methods proposed in [31, 32]. Of course, there is no denying that the methods proposed in [31, 32] have their advantages, for which one can refer to [31, 32].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 11301125), Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ14A010003), National Natural Science Foundation of China, Tian Yuan Special Foundation (Grant no. 11326231), and Research Foundation of Hangzhou Dianzi University (Grant no. KYS075613017).