Fast Computation of Singular Oscillatory Fourier Transforms

and Applied Analysis 3 (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 + εe iθ 1 , Γ3 : h3(l) = b + il, Γ4 : h4(l) = x + iL, x ∈ [a, b], Γ5 : h5(l) = a + il, and Γ6 : h6(l) = a + εe iθ 2 , where l ∈ [ε, L], θ1 ∈ [π/2, π], θ2 ∈ [0, π/2], and an arbitrary ε > 0. Moreover, (z − a)(b − z)f(z)e is continuous on all contours containing Γj, j = 1, 2, . . . , 6 (see Figure 1). By the Cauchy theorem [34], we obtain ∫ Γ 1 +Γ 2 +Γ 3 +Γ 4 +Γ 5 +Γ 6 (z − a)(b − z)f (z) edz = 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 existsM ≥ 0, such that 󵄨󵄨󵄨f (z) 󵄨󵄨󵄨 ≤ M, z ∈ G. (4) Since 0 ≤ 󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨 ∫ Γ 2 (z − a)(b − z)f (z) edz 󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨 ≤ ε ∫ π π/2 󵄨󵄨󵄨󵄨󵄨(b − a + εe iθ 1) 󵄨󵄨󵄨󵄨󵄨 󵄨󵄨󵄨󵄨f (b + εe iθ 1) 󵄨󵄨󵄨󵄨 dθ1 ≤ εM∫ π π/2 󵄨󵄨󵄨󵄨󵄨(b − a + εe iθ 1) 󵄨󵄨󵄨󵄨󵄨 dθ1 󳨀→ 0, as ε 󳨀→ 0, (5) then lim ε→0 ∫ Γ 2 (z − a)(b − z)f (z) edz = 0. (6) Similarly, it is easy to get lim ε→0 ∫ Γ 6 (z − a)(b − z)f (z) edz = 0. (7) 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)e is also analytic in the region G except G1 and G2. So there must exist a nonnegative number ?̃? such that 󵄨󵄨󵄨󵄨(z − a) (b − z)f (z) 󵄨󵄨󵄨󵄨 ≤ ?̃?, z ∈ G \ (G1 + G2) . (8) Therefore, 0 ≤ 󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨 ∫ Γ 4 (z − a)(b − z)f (z) edz 󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨


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][2][3][4], Clenshaw-Curtis-type [5,6], Filon-type [7,8], asymptotic [7], Levin [9,10], generalized quadrature rule [11,12], Levintype [13], and complex integration methods [14][15][16][17][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 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 wellknown 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 highfrequency 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 Filontype 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, ((!Γ( +  2 + 1)/(2)!)(1/2+ 1 +1 )),  1 = min{, }, and  2 = max{, }, is nearly twice as high as that of the Filontype 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.(1).In this subsection, we will study the problem of computing the integral (1) in depth.

A Decomposition of the Integral
Here, we provide the illustration of the special integration paths for the integrand (−)  (−)  ()  (see Figure 1).

Theorem 1. Suppose that 𝑓 is an analytic function in the region
Proof.Since () is analytic within the complex region , then the integrand ( − )  1).By the Cauchy theorem [34], we obtain 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      ()     ≤ ,  ∈ . Since then lim Similarly, it is easy to get lim Since the integrand ( − )  ( − )  () is analytic in the region  except  1 and  2 , then ( − )  ( − )  ()  is also analytic in the region  except  1 and  2 .So there must exist a nonnegative number M such that Since the Gamma function Γ() and the incomplete Gamma function Γ(, ) [35] are defined by which hold for the equality then Therefore, where there must exist a nonnegative  1 such that Thus, using (14), In the same way, we have lim Therefore, combining (3), ( 6), ( 9), ( 10), (16), and ( 17), we gain This completes the proof.
The succeeding part is to consider the efficient evaluation of the nethermost formulas of (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 Gausstype: 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: where the nodes  ()  are the zeros of the generalized or associated Laguerre polynomial  ()   () [1] and 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.
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.

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, where  1 () 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 ∫

Table 5 :
Absolute errors in -point approximations by the proposed method to the integral ∫