Research Article | Open Access
Fast Computation of Singular Oscillatory Fourier Transforms
We consider the problem of the numerical evaluation of singular oscillatory Fourier transforms , where . 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.
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 , Levin [9, 10], generalized quadrature rule [11, 12], Levin-type , 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 where is a sufficiently smooth function in , is a large parameter, and are real and finite, and , . In (1), if , , or , , the integrand has a singularity of a simple type at one sided endpoint or ; if , , the integrand has singularities of a complicated type at both endpoints of the interval; if, in particular, , 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 . 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  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  by Erdélyi. Three years later, a similar asymptotic result was reestablished by Lighthill in , 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 . 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  presented a Filon-type method and a Clenshaw-Curtis-Filon-type method for computing the integral (1), the error of which satisfies , where is the highest multiplicity of the Hermite interpolation at the endpoints and . In 2013, the recent literature  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, , , and , 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 .
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 , , and , respectively, where is a large number and the positive is small enough such that contains and .
Theorem 1. Suppose that is an analytic function in the region ; then
Proof. Since is analytic within the complex region , then the integrand is also analytic (has no singularity) in the region except and (), enclosed by , , defined by the following parametric forms , , , , , , , and , where , , , and an arbitrary . Moreover, is continuous on all contours containing , (see Figure 1). By the Cauchy theorem , 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 , such that Since then
Similarly, it is easy to get
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 Therefore, namely,
Since the Gamma function and the incomplete Gamma function  are defined by which hold for the equality then Therefore, where there must exist a nonnegative such that
Thus, using (14),
In the same way, we have
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).
2.2. Calculation of by the Generalized Gauss Laguerre Quadrature Rule
The key point in the interpolatory rule for the infinite interval is, of course, the rule of Gauss-type: where the and have been determined so that the formula is exact for a polynomial up to degree . 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 , . From , the generalized Gauss Laguerre quadrature rule is given, as follows: where the nodes are the zeros of the generalized or associated Laguerre polynomial  and
In (2), the integrands have a singularity of the form or , as . 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): 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 where and .
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 . The nodes and weights of the generalized Gauss-Laguerre quadrature rule are given in . 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, , 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 . In , the Chebyshev expansions method depends on ; for example, 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 . For details one can refer to . 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).