# Fast Computation of Singular Oscillatory Fourier Transforms

Authors: Hongchao Kang; Xinping Shao
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We consider the problem of the numerical evaluation of singular oscillatory Fourier transforms∫ab‍x-aαb-xβf(x)eiωxdx, where α>-1andβ>-1. Based on substituting the original interval of integration by the paths of steepest descent, if f is analytic in the complex region G containing [a, b], 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.


---

## Body

## 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 ∫ab‍f(x)eiωg(x)dx, 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)I[f]=∫ab‍(x-a)α(b-x)βf(x)eiωxdx,
where f(x) is a sufficiently smooth function in [a,b], ω is a large parameter, a and b 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 b or a; 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 f(x) is analytic in a region containing [a,b], 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 O(ω-s-r-2), where s is the highest multiplicity of the Hermite interpolation at the endpoints and r=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, O((n!Γ(n+r2+1)/(2n)!)(1/ω2n+r1+1)), r1=min⁡{α,β}, and r2=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 (x-a)α(b-x)βf(x)eiωx (see Figure 1). Here, let G, G1, and G2 denote the regions G={z∈C∣a≤R(z)≤b,0≤I(z)≤L}, G1={|z-a|≤ɛ,0≤arg(z)≤π/2}, and G2={|z-b|≤ɛ,π/2≤arg(z)≤π}, respectively, where L is a large number and the positive ɛ is small enough such that G contains G1 and G2.Figure 1
Illustration of the integration paths for the integrand(x-a)α(b-x)βf(x)eiωx.Theorem 1.
Suppose thatf is an analytic function in the region G={z∈C∣a≤R(z)≤b,I(z)≥0}; then
(2)I[f]=∫ab‍(x-a)α(b-x)βf(x)eiωxdx=iα+1eiωaωα+1∫0∞‍(b-a-ipω)βf(a+ipω)pαe-pdp-(-1)β+1iβ+1eiωbωβ+1×∫0∞‍(b-a+ipω)αf(b+ipω)pβe-pdp.Proof.
Sincef(z) is analytic within the complex region G, then the integrand (z-a)α(b-z)βf(z)eiωz is also analytic (has no singularity) in the region G except G1 and G2 (z∈G∖(G1+G2)), enclosed by Γj, j=1,2,…,6, defined by the following parametric forms Γ1:h1(l)=x, a+ɛ≤x≤b-ɛ, Γ2:h2(l)=b+ɛeiθ1, Γ3:h3(l)=b+il, Γ4:h4(l)=x+iL, x∈[a,b], Γ5:h5(l)=a+il, and Γ6:h6(l)=a+ɛeiθ2, where l∈[ɛ,L], θ1∈[π/2,π], θ2∈[0,π/2], and an arbitrary ɛ>0. Moreover, (z-a)α(b-z)βf(z)eiωz is continuous on all contours containing Γj, j=1,2,…,6 (see Figure 1). By the Cauchy theorem [34], we obtain
(3)∫Γ1+Γ2+Γ3+Γ4+Γ5+Γ6‍(z-a)α(b-z)βf(z)eiωzdz=0,
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 existsM≥0, such that
(4)|f(z)|≤M,z∈G.
Since
(5)0≤|∫Γ2‍(z-a)α(b-z)βf(z)eiωzdz|≤ɛβ+1∫π/2π‍|(b-a+ɛeiθ1)α||f(b+ɛeiθ1)|dθ1≤ɛβ+1M∫π/2π‍|(b-a+ɛeiθ1)α|dθ1⟶0,asɛ⟶0,
then
(6)lim⁡ɛ→0∫Γ2‍(z-a)α(b-z)βf(z)eiωzdz=0.
Similarly, it is easy to get(7)lim⁡ɛ→0∫Γ6‍(z-a)α(b-z)βf(z)eiωzdz=0.
Since the integrand(z-a)α(b-z)βf(z) is analytic in the region G except G1 and G2, then (z-a)α(b-z)βf(z)eiωz is also analytic in the region G except G1 and G2. So there must exist a nonnegative number M~ such that
(8)|(z-a)α(b-z)βf(z)|≤M~,z∈G∖(G1+G2).
Therefore,
(9)0≤|∫Γ4‍(z-a)α(b-z)βf(z)eiωzdz|≤e-Lω∫ab‍|(x-a+iL)α(b-x-iL)βf(x+iL)|dx≤(b-a)M~e-Lω⟶0,asL⟶∞;
namely,
(10)lim⁡L→∞∫Γ4‍(z-a)α(b-z)βf(z)eiωzdz=0.
Since the Gamma functionΓ(z) and the incomplete Gamma function Γ(s,z) [35] are defined by
(11)Γ(z)=∫0∞‍tz-1e-tdt,Γ(s,z)=∫s∞‍tz-1e-tdt,
which hold for the equality
(12)Γ(α+1)=Γ(α+1,0),
then
(13)∫0ɛ‍lαe-ωldl=1ωα+1∫0ωɛ‍tαe-tdt=1ωα+1[∫0∞‍tα+1-1e-tdt-∫ωɛ∞‍tα+1-1e-tdt]=1ωα+1[Γ(α+1)-Γ(α+1,ωɛ)]⟶0,asɛ⟶0.
Therefore,
(14)0≤|∫0ɛ‍(b-a-il)βlαf(a+il)e-ωldl|≤M·M1|∫0ɛ‍lαe-ωldl|⟶0,asɛ⟶0,
where there must exist a nonnegative M1 such that
(15)|(b-a-il)β|≤M1,l∈[0,ɛ].
Thus, using (14),
(16)lim⁡ɛ→0∫Γ5‍(z-a)α(b-z)βf(z)eiωzdz=lim⁡ɛ→0(-iα+1eiωa)×∫ɛL‍(b-a-il)βlαf(a+il)e-ωldl=lim⁡ɛ→0(-iα+1eiωa)×[∫0L‍(b-a-il)βlαf(a+il)e-ωldlh.h-∫0ɛ‍(b-a-il)βlαf(a+il)e-ωldl]=-iα+1eiωa∫0L‍(b-a-il)βlαf(a+il)e-ωldl=-iα+1eiωaωα+1×∫0ωL‍(b-a-ipω)αpαf(a+ipω)e-pdp.
In the same way, we have(17)lim⁡ɛ→0∫Γ3‍(z-a)α(b-z)βf(z)eiωzdz=(-1)β+1iβ+1eiωbωβ+1×∫0ωL‍(b-a+ipω)αf(b+ipω)pβe-pdp.
Therefore, combining (3), (6), (9), (10), (16), and (17), we gain
(18)∫ab‍(x-a)α(b-x)βf(x)eiωxdx=lim⁡ɛ→0∫Γ1‍(z-a)α(b-z)βf(z)eiωzdz=-lim⁡ɛ→0[lim⁡L→∞∫Γ2+Γ3+Γ4+Γ5+Γ6‍(z-a)α(b-z)βhhhhhhhhhhhhhhlim⁡L→∞∫Γ2+Γ3+Γ4+Γ5+Γ6‍×f(z)eiωzdz]=iα+1eiωaωα+1∫0∞‍(b-a-ipω)βf(a+ipω)pαe-pdp-(-1)β+1iβ+1eiωbωβ+1×∫0∞‍(b-a+ipω)αf(b+ipω)pβe-pdp.
This completes the proof.The succeeding part is to consider the efficient evaluation of the nethermost formulas of (2).


### 2.2. Calculation ofI[f] 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:
(19)∫0∞‍w(x)S(x)dx≈∑k=1n‍wkS(xk),
where the xk and wk have been determined so that the formula is exact for a polynomial S(x) up to degree 2n-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 functionxλe-x, λ>-1. From [1], the generalized Gauss Laguerre quadrature rule is given, as follows:
(20)∫0∞‍S(x)xλe-xdx=∑k=1n‍wk(λ)S(xk(λ))+n!Γ(n+λ+1)(2n)!S(2n)(ξ),hhhhhhhhhhhh0<ξ<∞,
where the nodes xk(λ) are the zeros of the generalized or associated Laguerre polynomial Ln(λ)(x) [1] and
(21)wk(λ)=Γ(n+λ+1)xk(λ)n![Ln+1(λ)(xk(λ))]2.In (2), the integrands have a singularity of the form pα or pβ, as p→0. Fortunately, this type of singularity can be handled efficiently by the generalized Gauss-Laguerre quadrature rule.


#### 2.2.1. Calculation ofI[f] by the Generalized Gauss Laguerre Quadrature Rule

The generalized Gauss-Laguerre formula (20) with n points xk(λ) and weights wk(λ) is applied to (2). Then we can easily obtain the following n-point approximation formula to (1):
(22)Qn[f]=iα+1eiωaωα+1×∑k=1n‍wk(α)(b-a-ixk(α)ω)βf(a+ixk(α)ω)-(-1)β+1iβ+1eiωbωβ+1×∑k=1n‍wk(β)(b-a+ixk(β)ω)αf(b+ixk(β)ω),
where xk(α) and xk(β) are the zeros of the generalized or associated Laguerre polynomials Ln(α)(x) and Ln(β)(x), respectively. And, wk(α) and wk(β) are given by choosing “λ=α” and “λ=β” in (21), respectively.Theorem 2.
Iff holds for the conditions of Theorem 1, then we have
(23)I[f]-Qn[f]~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),hhhhhhhhhhhhhhhhhhhhhhhhhhasω⟶∞,
where r1=min⁡{α,β} and r2=max⁡{α,β}.Proof.
Forr1=min⁡{α,β} and r2=max⁡{α,β}, since
(24)ω2n+r1+1=ω2n+α+1,ω2n+r1+1<ω2n+β+1,Γ(n+r2+1)=Γ(n+β+1),Γ(n+r2+1)>Γ(n+α+1),(α≤β);
or
(25)ω2n+r1+1=ω2n+β+1,ω2n+r1+1<ω2n+α+1,Γ(n+r2+1)=Γ(n+α+1),Γ(n+r2+1)>Γ(n+β+1),(α>β)
together with (2), (20), and (22), we have
(26)I[f]-Qn[f]=n!Γ(n+λ+1)(2n)!S(2n)(ξ)=iα+1eiωaωα+1n!Γ(n+α+1)(2n)!×d[(b-a-i(q/ω))βf(a+i(q/ω))]dq2n|q=ξ-(-1)β+1iβ+1eiωbωβ+1n!Γ(n+β+1)(2n)!×d[(b-a+i(q/ω))αf(b+i(q/ω))]dq2n|q=ξ~O(n!Γ(n+α+1)(2n)!ω-2n-α-1)-O(n!Γ(n+β+1)(2n)!ω-2n-β-1)~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),asω⟶∞.
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 a and b.In the subsequent section, we will make use of numerical examples to illustrate the absolute error and approximate values.


## 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 (x-a)α(b-x)βf(x)eiωx (see Figure 1). Here, let G, G1, and G2 denote the regions G={z∈C∣a≤R(z)≤b,0≤I(z)≤L}, G1={|z-a|≤ɛ,0≤arg(z)≤π/2}, and G2={|z-b|≤ɛ,π/2≤arg(z)≤π}, respectively, where L is a large number and the positive ɛ is small enough such that G contains G1 and G2.Figure 1
Illustration of the integration paths for the integrand(x-a)α(b-x)βf(x)eiωx.Theorem 1.
Suppose thatf is an analytic function in the region G={z∈C∣a≤R(z)≤b,I(z)≥0}; then
(2)I[f]=∫ab‍(x-a)α(b-x)βf(x)eiωxdx=iα+1eiωaωα+1∫0∞‍(b-a-ipω)βf(a+ipω)pαe-pdp-(-1)β+1iβ+1eiωbωβ+1×∫0∞‍(b-a+ipω)αf(b+ipω)pβe-pdp.Proof.
Sincef(z) is analytic within the complex region G, then the integrand (z-a)α(b-z)βf(z)eiωz is also analytic (has no singularity) in the region G except G1 and G2 (z∈G∖(G1+G2)), enclosed by Γj, j=1,2,…,6, defined by the following parametric forms Γ1:h1(l)=x, a+ɛ≤x≤b-ɛ, Γ2:h2(l)=b+ɛeiθ1, Γ3:h3(l)=b+il, Γ4:h4(l)=x+iL, x∈[a,b], Γ5:h5(l)=a+il, and Γ6:h6(l)=a+ɛeiθ2, where l∈[ɛ,L], θ1∈[π/2,π], θ2∈[0,π/2], and an arbitrary ɛ>0. Moreover, (z-a)α(b-z)βf(z)eiωz is continuous on all contours containing Γj, j=1,2,…,6 (see Figure 1). By the Cauchy theorem [34], we obtain
(3)∫Γ1+Γ2+Γ3+Γ4+Γ5+Γ6‍(z-a)α(b-z)βf(z)eiωzdz=0,
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 existsM≥0, such that
(4)|f(z)|≤M,z∈G.
Since
(5)0≤|∫Γ2‍(z-a)α(b-z)βf(z)eiωzdz|≤ɛβ+1∫π/2π‍|(b-a+ɛeiθ1)α||f(b+ɛeiθ1)|dθ1≤ɛβ+1M∫π/2π‍|(b-a+ɛeiθ1)α|dθ1⟶0,asɛ⟶0,
then
(6)lim⁡ɛ→0∫Γ2‍(z-a)α(b-z)βf(z)eiωzdz=0.
Similarly, it is easy to get(7)lim⁡ɛ→0∫Γ6‍(z-a)α(b-z)βf(z)eiωzdz=0.
Since the integrand(z-a)α(b-z)βf(z) is analytic in the region G except G1 and G2, then (z-a)α(b-z)βf(z)eiωz is also analytic in the region G except G1 and G2. So there must exist a nonnegative number M~ such that
(8)|(z-a)α(b-z)βf(z)|≤M~,z∈G∖(G1+G2).
Therefore,
(9)0≤|∫Γ4‍(z-a)α(b-z)βf(z)eiωzdz|≤e-Lω∫ab‍|(x-a+iL)α(b-x-iL)βf(x+iL)|dx≤(b-a)M~e-Lω⟶0,asL⟶∞;
namely,
(10)lim⁡L→∞∫Γ4‍(z-a)α(b-z)βf(z)eiωzdz=0.
Since the Gamma functionΓ(z) and the incomplete Gamma function Γ(s,z) [35] are defined by
(11)Γ(z)=∫0∞‍tz-1e-tdt,Γ(s,z)=∫s∞‍tz-1e-tdt,
which hold for the equality
(12)Γ(α+1)=Γ(α+1,0),
then
(13)∫0ɛ‍lαe-ωldl=1ωα+1∫0ωɛ‍tαe-tdt=1ωα+1[∫0∞‍tα+1-1e-tdt-∫ωɛ∞‍tα+1-1e-tdt]=1ωα+1[Γ(α+1)-Γ(α+1,ωɛ)]⟶0,asɛ⟶0.
Therefore,
(14)0≤|∫0ɛ‍(b-a-il)βlαf(a+il)e-ωldl|≤M·M1|∫0ɛ‍lαe-ωldl|⟶0,asɛ⟶0,
where there must exist a nonnegative M1 such that
(15)|(b-a-il)β|≤M1,l∈[0,ɛ].
Thus, using (14),
(16)lim⁡ɛ→0∫Γ5‍(z-a)α(b-z)βf(z)eiωzdz=lim⁡ɛ→0(-iα+1eiωa)×∫ɛL‍(b-a-il)βlαf(a+il)e-ωldl=lim⁡ɛ→0(-iα+1eiωa)×[∫0L‍(b-a-il)βlαf(a+il)e-ωldlh.h-∫0ɛ‍(b-a-il)βlαf(a+il)e-ωldl]=-iα+1eiωa∫0L‍(b-a-il)βlαf(a+il)e-ωldl=-iα+1eiωaωα+1×∫0ωL‍(b-a-ipω)αpαf(a+ipω)e-pdp.
In the same way, we have(17)lim⁡ɛ→0∫Γ3‍(z-a)α(b-z)βf(z)eiωzdz=(-1)β+1iβ+1eiωbωβ+1×∫0ωL‍(b-a+ipω)αf(b+ipω)pβe-pdp.
Therefore, combining (3), (6), (9), (10), (16), and (17), we gain
(18)∫ab‍(x-a)α(b-x)βf(x)eiωxdx=lim⁡ɛ→0∫Γ1‍(z-a)α(b-z)βf(z)eiωzdz=-lim⁡ɛ→0[lim⁡L→∞∫Γ2+Γ3+Γ4+Γ5+Γ6‍(z-a)α(b-z)βhhhhhhhhhhhhhhlim⁡L→∞∫Γ2+Γ3+Γ4+Γ5+Γ6‍×f(z)eiωzdz]=iα+1eiωaωα+1∫0∞‍(b-a-ipω)βf(a+ipω)pαe-pdp-(-1)β+1iβ+1eiωbωβ+1×∫0∞‍(b-a+ipω)αf(b+ipω)pβe-pdp.
This completes the proof.The succeeding part is to consider the efficient evaluation of the nethermost formulas of (2).


## 2.2. Calculation ofI[f] 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:
(19)∫0∞‍w(x)S(x)dx≈∑k=1n‍wkS(xk),
where the xk and wk have been determined so that the formula is exact for a polynomial S(x) up to degree 2n-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 functionxλe-x, λ>-1. From [1], the generalized Gauss Laguerre quadrature rule is given, as follows:
(20)∫0∞‍S(x)xλe-xdx=∑k=1n‍wk(λ)S(xk(λ))+n!Γ(n+λ+1)(2n)!S(2n)(ξ),hhhhhhhhhhhh0<ξ<∞,
where the nodes xk(λ) are the zeros of the generalized or associated Laguerre polynomial Ln(λ)(x) [1] and
(21)wk(λ)=Γ(n+λ+1)xk(λ)n![Ln+1(λ)(xk(λ))]2.In (2), the integrands have a singularity of the form pα or pβ, as p→0. Fortunately, this type of singularity can be handled efficiently by the generalized Gauss-Laguerre quadrature rule.


### 2.2.1. Calculation ofI[f] by the Generalized Gauss Laguerre Quadrature Rule

The generalized Gauss-Laguerre formula (20) with n points xk(λ) and weights wk(λ) is applied to (2). Then we can easily obtain the following n-point approximation formula to (1):
(22)Qn[f]=iα+1eiωaωα+1×∑k=1n‍wk(α)(b-a-ixk(α)ω)βf(a+ixk(α)ω)-(-1)β+1iβ+1eiωbωβ+1×∑k=1n‍wk(β)(b-a+ixk(β)ω)αf(b+ixk(β)ω),
where xk(α) and xk(β) are the zeros of the generalized or associated Laguerre polynomials Ln(α)(x) and Ln(β)(x), respectively. And, wk(α) and wk(β) are given by choosing “λ=α” and “λ=β” in (21), respectively.Theorem 2.
Iff holds for the conditions of Theorem 1, then we have
(23)I[f]-Qn[f]~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),hhhhhhhhhhhhhhhhhhhhhhhhhhasω⟶∞,
where r1=min⁡{α,β} and r2=max⁡{α,β}.Proof.
Forr1=min⁡{α,β} and r2=max⁡{α,β}, since
(24)ω2n+r1+1=ω2n+α+1,ω2n+r1+1<ω2n+β+1,Γ(n+r2+1)=Γ(n+β+1),Γ(n+r2+1)>Γ(n+α+1),(α≤β);
or
(25)ω2n+r1+1=ω2n+β+1,ω2n+r1+1<ω2n+α+1,Γ(n+r2+1)=Γ(n+α+1),Γ(n+r2+1)>Γ(n+β+1),(α>β)
together with (2), (20), and (22), we have
(26)I[f]-Qn[f]=n!Γ(n+λ+1)(2n)!S(2n)(ξ)=iα+1eiωaωα+1n!Γ(n+α+1)(2n)!×d[(b-a-i(q/ω))βf(a+i(q/ω))]dq2n|q=ξ-(-1)β+1iβ+1eiωbωβ+1n!Γ(n+β+1)(2n)!×d[(b-a+i(q/ω))αf(b+i(q/ω))]dq2n|q=ξ~O(n!Γ(n+α+1)(2n)!ω-2n-α-1)-O(n!Γ(n+β+1)(2n)!ω-2n-β-1)~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),asω⟶∞.
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 a and b.In the subsequent section, we will make use of numerical examples to illustrate the absolute error and approximate values.


## 2.2.1. Calculation ofI[f] by the Generalized Gauss Laguerre Quadrature Rule

The generalized Gauss-Laguerre formula (20) with n points xk(λ) and weights wk(λ) is applied to (2). Then we can easily obtain the following n-point approximation formula to (1):
(22)Qn[f]=iα+1eiωaωα+1×∑k=1n‍wk(α)(b-a-ixk(α)ω)βf(a+ixk(α)ω)-(-1)β+1iβ+1eiωbωβ+1×∑k=1n‍wk(β)(b-a+ixk(β)ω)αf(b+ixk(β)ω),
where xk(α) and xk(β) are the zeros of the generalized or associated Laguerre polynomials Ln(α)(x) and Ln(β)(x), respectively. And, wk(α) and wk(β) are given by choosing “λ=α” and “λ=β” in (21), respectively.Theorem 2.
Iff holds for the conditions of Theorem 1, then we have
(23)I[f]-Qn[f]~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),hhhhhhhhhhhhhhhhhhhhhhhhhhasω⟶∞,
where r1=min⁡{α,β} and r2=max⁡{α,β}.Proof.
Forr1=min⁡{α,β} and r2=max⁡{α,β}, since
(24)ω2n+r1+1=ω2n+α+1,ω2n+r1+1<ω2n+β+1,Γ(n+r2+1)=Γ(n+β+1),Γ(n+r2+1)>Γ(n+α+1),(α≤β);
or
(25)ω2n+r1+1=ω2n+β+1,ω2n+r1+1<ω2n+α+1,Γ(n+r2+1)=Γ(n+α+1),Γ(n+r2+1)>Γ(n+β+1),(α>β)
together with (2), (20), and (22), we have
(26)I[f]-Qn[f]=n!Γ(n+λ+1)(2n)!S(2n)(ξ)=iα+1eiωaωα+1n!Γ(n+α+1)(2n)!×d[(b-a-i(q/ω))βf(a+i(q/ω))]dq2n|q=ξ-(-1)β+1iβ+1eiωbωβ+1n!Γ(n+β+1)(2n)!×d[(b-a+i(q/ω))αf(b+i(q/ω))]dq2n|q=ξ~O(n!Γ(n+α+1)(2n)!ω-2n-α-1)-O(n!Γ(n+β+1)(2n)!ω-2n-β-1)~O(n!Γ(n+r2+1)(2n)!1ω2n+r1+1),asω⟶∞.
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 a and b.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 xk(λ) and weights wk(λ) 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.Tables1, 2, 3, 4, 5, and 6 demonstrate the absolute errors and approximate values in n-point approximations by the proposed method to the above integrals. Furthermore, they exhibit the fast convergence of the approximations as n increases. Meanwhile, all these tables above show that more and more accurate approximations can be obtained as ω increases and n is fixed. Conversely, as n 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, n=2,3,4, and so forth. Also, for sufficiently large ω and some n, 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 n, 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,
(27)N1(ω)=⌊1.05ω+57.5⌋,
where N1(ω) 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 inn-point approximations by the proposed method to the integral ∫01‍x-(1/2)(x-1)-(1/2)cos⁡xeiωxdx. (b) Absolute errors in n-point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral ∫01‍x-(1/2)(x-1)-(1/2)cos⁡xeiωxdx.
(a)
ω
∖
n
2
3
4
1
0
2.7
×
1
0
-
5
1.0
×
1
0
-
6
3.9
×
1
0
-
8
20
1.1
×
1
0
-
6
9.9
×
1
0
-
9
2.2
×
1
0
-
10
40
6.3
×
1
0
-
8
2.4
×
1
0
-
10
1.6
×
1
0
-
12
80
1.1
×
1
0
-
9
1.3
×
1
0
-
12
6.3
×
1
0
-
15
(b)
ω
∖
n
2
3
4
10
8.6
×
1
0
-
1
7.8
×
1
0
-
1
6.9
×
1
0
-
1
20
7.6
×
1
0
-
1
5.1
×
1
0
-
1
3.4
×
1
0
-
1
40
6.3
×
1
0
-
1
4.6
×
1
0
-
1
3.0
×
1
0
-
1
80
5.9
×
1
0
-
1
4.4
×
1
0
-
1
2.2
×
1
0
-
1Table 2
(a) Approximate values in2-point approximations by the proposed method to the integral ∫01‍x-(1/2)(x-1)-(1/2)cos⁡xeiωxdx, compared with exact values. (b) Approximate values in 2-point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral ∫01‍x-(1/2)(x-1)-(1/2)cos⁡xeiωxdx, compared with exact values.
(a)
ω
1
000
2000
Exact values
0.06937869322675
+
0.04527657820424
i
0.03655011452354
+
0.04767073962282
i
Approximate values
0.06937869322672
+
0.04527657820423
i
0.03655011452354
+
0.04767073962281
i
ω
3000
5000
Exact values
0.01353266279337
+
0.03766001599991
i
0.00974123171200
+
0.00678449969265
i
Approximate values
0.01353266279337
+
0.03766001599991
i
0.00974123171200
+
0.00678449969265
i
(b)
ω
1000
2000
Exact values
0.06937869322675
+
0.04527657820424
i
0.03655011452354
+
0.04767073962282
i
Approximate values
0.06937879634322
+
0.04527685643561
i
0.03655017563845
+
0.0476703845971
i
ω
3000
5000
Exact values
0.01353266279337
+
0.03766001599991
i
0.00974123171200
+
0.00678449969265
i
Approximate values
0.01353266884523
+
0.03766001457824
i
0.00974123131452
+
0.00678449989714
iTable 3
Absolute errors inn-point approximations by the proposed method to the integral ∫23‍(x-2)-(1/2)(3-x)-(1/2)sin⁡⁡xeiωxdx.
ω
∖
n
2
3
4
8
4.5
×
1
0
-
5
2.8
×
1
0
-
6
3.0
×
1
0
-
7
16
2.4
×
1
0
-
6
5.2
×
1
0
-
8
1.8
×
1
0
-
9
32
1.1
×
1
0
-
7
6.2
×
1
0
-
10
6.7
×
1
0
-
12
64
5.2
×
1
0
-
9
7.7
×
1
0
-
12
2.5
×
1
0
-
14Table 4
Approximate values in2-point approximations by the proposed method to the integral ∫23‍(x-2)-(1/2)(3-x)-(1/2)sin⁡⁡xeiωxdx, compared with exact values.
ω
1000
2000
Exact values
-
0.05096624886110
+
0.02697237461921
i
0.00069025638199
-
0.04128898797276
i
Approximate values
-
0.05096624886109
+
0.02697237461920
i
0.00069025638199
-
0.04128898797276
i
ω
3000
5000
Exact values
0.02715481058125
+
0.01444143741392
i
-
0.00930648155343
-
0.01691388809588
i
Approximate values
0.02715481058125
+
0.01444143741392
i
-
0.00930648155343
-
0.01691388809588
iTable 5
Absolute errors inn-point approximations by the proposed method to the integral ∫-11‍(x+1)-(1/2)(1-x)-(1/2)exeiωxdx.
ω
∖
n
2
3
4
5
1.1
×
1
0
-
4
3.5
×
1
0
-
6
2.0
×
1
0
-
7
10
7.3
×
1
0
-
6
5.0
×
1
0
-
8
1.0
×
1
0
-
9
20
3.4
×
1
0
-
7
7.2
×
1
0
-
10
4.2
×
1
0
-
12
40
1.5
×
1
0
-
9
5.3
×
1
0
-
12
9.5
×
1
0
-
15Table 6
Approximate values in2-point approximations by the proposed method to the integral ∫-11‍(x+1)-(1/2)(1-x)-(1/2)exeiωxdx, compared with exact values.
ω
500
1000
Exact values
-
0.16534892514854
+
0.03862546212408
i
0.12015044233980
+
0.01747126126322
i
Approximate values
-
0.16534892514837
+
0.03862546212408
i
0.12015044233979
+
0.01747126126322
i
ω
3000
5000
Exact values
-
0.03778034068446
+
0.04543591233617
i
-
0.03222911938840
-
0.03366226121525
i
Approximate values
-
0.03778034068446
+
0.04543591233617
i
-
0.03222911938840
-
0.03366226121525
i


## 4. Concluding Remark

In this paper, for an analytic and not increasing fastf in G={z∈C∣a≤R(z)≤b,I(z)≥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].



---

*Source: 984834-2014-07-17.xml*