Abstract and Applied Analysis

Volume 2014 (2014), Article ID 984834, 8 pages

http://dx.doi.org/10.1155/2014/984834

Research Article

Fast Computation of Singular Oscillatory Fourier Transforms

Hongchao Kang and Xinping Shao

Department of Mathematics, School of Science, Hangzhou Dianzi University, Hangzhou, Zhejiang 310018, China

Received 1 April 2014; Accepted 3 July 2014; Published 17 July 2014

Academic Editor: Chuanzhi Bai

Copyright © 2014 Hongchao Kang and Xinping Shao. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

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

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 [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 , where is the highest multiplicity of the Hermite interpolation at the endpoints and . 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, , , 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 [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 , , and , 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 ; 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 [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 , 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 [35] 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 [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 . 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 .

Proof. For and , since

or

together with (2), (20), and (22), we have

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, , 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 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 . (b) Absolute errors in -point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral .

(a)  

(b)  

Table 2: (a) Approximate values in -point approximations by the proposed method to the integral , compared with exact values. (b) Approximate values in -point approximations by the Clenshaw-Curtis-Filon-type method [31] to the integral , compared with exact values.

(a)  

Exact values

Approximate values

Exact values

Approximate values

(b)  

Exact values

Approximate values

Exact values

Approximate values

Table 3: Absolute errors in -point approximations by the proposed method to the integral .

Table 4: Approximate values in -point approximations by the proposed method to the integral , compared with exact values.

Exact values

Approximate values

Exact values

Approximate values

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

Table 6: Approximate values in -point approximations by the proposed method to the integral , compared with exact values.

Exact values

Approximate values

Exact values

Approximate values

4. Concluding Remark

In this paper, for an analytic and not increasing fast in , 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).

References

P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, NY, USA, 2nd edition, 1984.

L. N. G. Filon, “On a quadrature formula for trigonometric integrals,” Proceedings of the Royal Society of Edinburgh, vol. 49, pp. 38–47, 1928.

E. A. Flinn, “A modification of Filon's method of numerical integration,” Journal of the Association for Computing Machinery, vol. 7, pp. 181–184, 1960.

Y. L. Luke, “On the computation of oscillatory integrals,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 50, pp. 269–277, 1954.

C. W. Clenshaw and A. R. Curtis, “A method for numerical integration on an automatic computer,” Numerische Mathematik, vol. 2, pp. 197–205, 1960.

S. Xiang, X. Chen, and H. Wang, “Error bounds for approximation in Chebyshev points,” Numerische Mathematik, vol. 116, no. 3, pp. 463–491, 2010.

A. Iserles and S. P. Nørsett, “Efficient quadrature of highly oscillatory integrals using derivatives,” Proceedings of The Royal Society of London A, vol. 461, no. 2057, pp. 1383–1399, 2005.

S. Xiang, “Efficient Filon-type methods for ,” Numerische Mathematik, vol. 105, no. 4, pp. 633–658, 2007.

D. Levin, “Procedures for computing one- and two-dimensional integrals of functions with rapid irregular oscillations,” Mathematics of Computation, vol. 38, no. 158, pp. 531–538, 1982.

D. Levin, “Analysis of a collocation method for integrating rapidly oscillatory functions,” Journal of Computational and Applied Mathematics, vol. 78, no. 1, pp. 131–138, 1997.

K. C. Chung, G. A. Evans, and J. R. Webster, “A method to generate generalized quadrature rules for oscillatory integrals,” Applied Numerical Mathematics, vol. 34, no. 1, pp. 85–93, 2000.

G. A. Evans and K. C. Chung, “Evaluating infinite range oscillatory integrals using generalised quadrature methods,” Applied Numerical Mathematics, vol. 57, no. 1, pp. 73–79, 2007.

S. Olver, “Moment-free numerical integration of highly oscillatory functions,” IMA Journal of Numerical Analysis, vol. 26, no. 2, pp. 213–227, 2006.

N. Bleistein and R. Handelsman, Asymptotic Expansions of Integrals, Holt, Rinehart and Winston, New York, NY, USA, 1975.

A. Erdélyi, Asymptotic Expansions, Dover, New York, NY, USA, 1956.

R. Wong, Asymptotic Approximations of Integrals, SIAM, Philadelphia, PA, USA, 2001.

H. Wang and S. Xiang, “On the evaluation of Cauchy principal value integrals of oscillatory functions,” Journal of Computational and Applied Mathematics, vol. 234, no. 1, pp. 95–100, 2010.

H. Kang and S. Xiang, “On the calculation of highly oscillatory integrals with an algebraic singularity,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 3890–3897, 2010.

H. Brunner, “On the numerical solution of first-kind Volterra integral equations with highly oscillatory kernels,” in Highly Oscillatory Problems: From Theory to Applications, pp. 13–17, Isaac Newton Institute, 2010.

H. Brunner, Open Problems in the Computational Solution of Volterra Functional Equations with Highly Oscillatory Kernels, HOP 2007: Effective Computational Methods for Highly Oscillatory Solutions, Isaac Newton Institute for Mathematical Sciences, Cambridge, UK, 2007.

S. Lu, “A class of oscillatory singular integrals,” International Journal of Applied Mathematical Sciences, vol. 2, no. 1, pp. 47–64, 2005.

D. H. Phong and E. M. Stein, “Singular integrals related to the Radon transform and boundary value problems,” Proceedings of the National Academy of Sciences of the United States of America, vol. 80, no. 24, pp. 7697–7701, 1983.

M. A. Hamed and B. Cummins, “Numerical integration formula for the solution of the singular integral equation for classical crack problems in plane and antiplane elasticity,” Journal of King Saud University, vol. 3, no. 2, pp. 217–231, 1991.

V. Domínguez, I. G. Graham, and V. P. Smyshlyaev, “A hybrid numerical-asymptotic boundary integral method for high-frequency acoustic scattering,” Numerische Mathematik, vol. 106, no. 3, pp. 471–510, 2007.

V. Domínguez, I. G. Graham, and V. P. Smyshlyaev, “Stability and error estimates for Filon-Clenshaw-CURtis rules for highly oscillatory integrals,” IMA Journal of Numerical Analysis, vol. 31, no. 4, pp. 1253–1280, 2011.

A. Erdélyi, “Asymptotic representations of Fourier integrals and the method of stationary phase,” Journal of the Society for Industrial and Applied Mathematics, vol. 3, pp. 17–27, 1955.

M. L. Lighthill, Introduction to Fourier Analisis and Generalized Functions, Cambrige University Press, New York, NY, USA, 1958.

J. N. Lyness, “Adjusted forms of the Fourier coefficient asymptotic expansion and applications in numerical quadrature,” Mathematics of Computation, vol. 25, pp. 87–104, 1971.

J. N. Lyness, “Numerical evaluation of a fixed-amplitude variable-phase integral,” Numerical Algorithms, vol. 49, no. 1–4, pp. 235–249, 2008.

J. N. Lyness and J. W. Lottes, “Asymptotic expansions for oscillatory integrals using inverse functions,” BIT: Numerical Mathematics, vol. 49, no. 2, pp. 397–417, 2009.

H. Kang and S. Xiang, “Efficient integration for a class of highly oscillatory integrals,” Applied Mathematics and Computation, vol. 218, no. 7, pp. 3553–3564, 2011.

H. Kang, S. Xiang, and G. He, “Computation of integrals with oscillatory and singular integrands using Chebyshev expansions,” Journal of Computational and Applied Mathematics, vol. 242, pp. 141–156, 2013.

S. Xiang, Y. Cho, H. Wang, and H. Brunner, “Clenshaw-curtis-filon-type methods for highly oscillatory bessel transforms and applications,” IMA Journal of Numerical Analysis, vol. 31, no. 4, pp. 1281–1314, 2011.

P. Henrici, Applied and Computational Complex Analysis, vol. 1, John Wiley & Sons, New York, NY, USA, 1974.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, USA, 1964.