Fast Computation of Integrals with Fourier-Type Oscillator Involving Stationary Point

: An adaptive splitting algorithm was implemented for numerical evaluation of Fourier-type highly oscillatory integrals involving stationary point. Accordingly, a modiﬁed Levin collocation method was coupled with multi-resolution quadratures in order to tackle the stationary point and irregular oscillations of the integrand caused by ω . Some test problems are included to verify the accuracy of the proposed methods.


Introduction
Accurate and efficient evaluation of highly oscillatory integrals is challenging as the analytical and classical computational methods fail to compute the integrals. The interest of computational scientists to evaluate these integrals quickly and accurately is developed due to the wide range of applications of these integrals in different fields of science and engineering such as optics, acoustics, quantum mechanics, seismology image processing, and electromagnetic [1][2][3][4]. Generally, these integrals can be written as where r, Θ are both smooth and non-oscillatory functions over the interval [−1, 1] , ω 1 and the function Θ(x) has a stationary point of order k at x = x 0 i.e., Θ (x 0 ) = Θ (x 0 ) = ... = Θ (k−1) (x 0 ) = 0 and Θ (k) (x 0 ) = 0. Large values of the frequency ω cause a highly oscillations of the integrand, which make the existing quadratures and softwares like MATHEMATICA in some cases, inaccurate. New algorithms need to be developed for accurate and efficient computation of these integrals.
In the last two decades, a number of accurate and efficient methods have been designed for numerical evaluation of one-dimensional highly oscillatory integrals, which include: the asymptotic method [5][6][7][8], the numerical steepest descent method [9], the Filon(-type) methods [10,11], and the Levin(-type) methods [12][13][14][15][16]. Among which, the Levin method has attracted much attention as it can handle highly oscillatory integrals with complicated phase functions. The asymptotic expansion theory which is considered in [10] is also one of the best methods for evaluation of highly oscillatory integrals. Since then, many authors [11,17,18] have implemented this method. The Filon method [11] gives better approximation when the frequency ω −→ ∞, but the method has its limitations, that is, it is only applicable for linear phase functions. Compared to the Filon method, Levin's method has the advantage of solving oscillatory integrals with complicated phase functions. According to the Levin method, a one-dimensional oscillatory integral is transformed into an ordinary differential equation (ODE) and then the ODE can be solved by the collocation method. Subsequently, a solution of such an integral is obtained.
In [13], the authors presented an improved Levin method with Chebyshev polynomials as a basis function and better accuracy was obtained. In [15], the authors evaluated non-oscillatory, mildly oscillatory, and highly oscillatory multi-dimensional integrals by using a meshless method based on a multi-quadric radial basis function and multi-resolution quadratures. In [14], the authors presented a meshless method based on a multi-quadric radial basis function (MQ-RBF) and multi-resolution quadrature to solve highly oscillatory integrals with and without critical points. In [19], the author split the integration domain into two or more sub-domains in order to separate the critical point and a two-point Gauss-Legendre quadrature was used to evaluate the integrals with critical points. The remaining integrals were computed by the modified Levin collocation method.
In [16], the authors developed an algorithm in which the meshless method with a multi-quadric radial basis function was coupled with the multi-resolution quadrature based on hybrid functions and the Haar wavelet to evaluate highly oscillatory integrals with critical point(s). The multi-resolution quadrature based on hybrid function of order 8 to evaluate integral of the form I where h = b−a 8n . For a = 0, b = 1, and n = 4, the truncation error bound of the hybrid function for some ξ 1 in (a, b) [20]. The formula of the Haar wavelet based quadrature for computing the same integral is given as where h = b−a n and n = 2M. The truncation error of the quadrature rule Q w H [ f ] is given by for some ς 1 ∈ (a, b) [20]. The current work is the extension of the work reported in [16,19]. According to the proposed procedure, the interval of integration is divided into two or more subintervals in order to isolate the stationary point. Integrals defined over a stationary point oriented interval can be computed by the multi-resolution quadratures even for higher frequencies, while the remaining integrals with high oscillations can be approximated by the Levin collocation quadrature. Theoretical error analysis of the individual methods and the proposed algorithm was performed. Numerical results met the theoretical proofs of the proposed algorithm.

Evaluation Procedure
Levin type quadratures fail to compute Fourier-type oscillatory integrals (1) containing a stationary point. We propose a splitting algorithm, which couples two types of quadrature rules to tackle the stationary point. In this algorithm, we subdivide the interval [−1, 1] in such a manner that the stationary point is separated in a very small interval and can be handled by the multi-resolution quadratures. The evaluation procedure is discussed as follows.

Levin Quadrature
The Levin quadrature is applicable to evaluate Fourier-type oscillatory integrals of the form where Θ (x) = 0. To find the approximate solution of the integral (5), an approximate functionŜ(x) is assumed to satisfy the following differential equation where L is defined as On applying the interpolation condition to (6), one can find the approximate solutionŜ(x). Consequently, the desired solution of integral (5) can be obtained as

Chebyshev Differentiation Matrix and Its Approximation
The Chebyshev differentiation matrix D [13] is used to approximate the derivative of a function and can be explicitly defined as where and x j = cos( πj N−1 ), j = 0, 1, 2, ..., N − 1 are the N Chebyshev-Gauss-Lobatto nodes. If an integral is defined over any domain [a, b], then, the Chebyshev-Gauss-Lobatto nodes x j can be obtained by the The procedure is applied to compute the first derivative of the following oscillatory type functions. The results in terms of absolute errors are shown in Figure 1. From Figure 1, it is shown that the accuracy for the derivative of the functions approximated by the Chebyshev differentiation matrix Du at Chebyshev-Gauss-Lobatto nodes improves on increasing N, and decaying for large frequency ω. In the Levin procedure, we are concerned with the derivative of a smooth function S(x), which is independent of ω.

Chebyshev-Levin Quadrature
The Chebyshev differentiation matrix is used to find the approximate solutionŜ(x) of ODE (6). If Θ(x) is any phase function of integral (1) such that Θ (x) = 0 for all x ∈ [a, b], then the discretized form of the ODE (6) at Chebyshev-Gauss-Lobatto nodes is given bŷ where represents Hadamard product of the real valued function Θ and the vectorŜ. IfŜ = DŜ, then the descretized form ODE (10) is given as where ∑ = diag(Θ ) and can be written in Matlab built-in form as (11) can be written as where I is the identity matrix and each ofŜ and R are column matrices of order N × 1. The vectorŝ S, Θ and R are obtained from the real valued functionsŜ(x), Θ(x) and r(x), respectively, at Chebyshev-Gauss-Lobatto nodes aŝ Equation (12) can be written in matrix form as where A is a square matrix of order N × N and both column vectorsŜ and R have length N. Then the system of Equation (13) for the unknown vectorŜ(x) is solved. Subsequently, the desired numerical solution of the oscillatory integral (5) can be obtained by the Levin's procedure as

Adaptive Splitting
The Chebyshev-Levin quadrature Q C−L [r] fails to evaluate the Fourier-type oscillatory integral (1) with a stationary point at x = x 0 . In this section, an adaptive splitting is used to evaluate the integral (1). According to this splitting, the domain interval is bifurcated into two or more sub-intervals in such manner that the stationary point is isolated in a very small interval. For this, a parameter ζ = ( n 10ω ) 1 k , a point in the vicinity of stationary point x 0 , is defined such that for fixed n, ζ → 0 as ω → ∞.
Case-I When the stationary point x 0 = a lies at the left end point of the domain, then integral (1) can be split as The integral As the stationary point occurs anywhere in the interval [a, b] and hence for this, we have to discuss the following case as case-II.
Case-II When a < x 0 < b, then the split form of the integral (1) is given as The integral I 2 [r, ω] contains the stationary point and is computed by the multi-resolution with n quadrature points. The integrals containing no stationary point are computed by Q C−L [r] with N collocation points. The following transformation Equation (15) for integral I 3 [r, ω] is used to descretize the interval in Chebyshev-Gauss-Lobatto nodes.

Error Bounds
In this section, error bounds of the individual methods as well as the proposed Algorithm-I in the inverse power of ω are obtained to ensure the asymptotic convergence rate of the new algorithm.
Lemma 3. Suppose that the integral I 2 (r, ω) of (14) has a unique stationary point at x = x 0 which lies at the left end of the domain. Then the error bound of the Chebyshev-Levin quadrature Q C−L [r] with N collocation points for computing I 2 is given by where k is the order of stationary point.

Proof.
Letp(x) = ∑ N−1 i=0 λ i T i (x) be the approximate Chebyshev interpolating polynomial of the first kind of order N − 1 that satisfies the following ODE Here λ i , i = 0, 1, ...N − 1 are the N unknown coefficients, which can be determined by the following interpolation condition where ψ(x) =p (x) + iωΘ (x)p(x). By applying the recurrence relation T i (x) = 2xT i−1 (x) − T i−2 (x), i = 2, 3, ... [20] with T 0 = 1 and T 1 = x, it follows that the function ψ(x) can be written as ). Using Lemma 2 and [16], the error bound of the integral I 2 (r, ω) of (14) computed by the Q C−L [r] is given as Following the procedure [16] and (17), we can find the error bound of the proposed method r(x)e iωΘ(x) dx with n quadrature points is given by where n = ξ−x 0 8h .

Proof
where a is the stationary point, C is a constant independent of ω, h and k is the order of stationary point. Since Q 8 h [ f ] can compute the integral I 1 (r, ω) having stationary point at x = x 0 . Therefore, the error bound of the hybrid function Q 8 h [ f ] for computing I 1 (r, ω) = ζ x 0 r(x)e iωΘ(x) dx, (where ζ is the nearest point of x 0 ) is given as If x 0 = 0, then (18) becomes where C 1 is independent of h and ω.
Since n = ζ 8h , the inequality in (19) can be written as Using ζ = ( n 10ω ) 1 k , Equation (20) can be written as As the Chebyshev-hybrid quadrature splitting method (ChQ) is obtained as therefore, the error bound of the splitting method ChQ can be obtained as where C 2 and C are independent of h and ω. For each subinterval of size h, Equation (21) can be written as

Numerical Examples and Discussion
In this section, the proposed method was tested on solving some benchmark problems [11,19]. The reference solution was obtained by MAPLE 18. The absolute errors L abs and scaled absolute errors were computed in each test problem. Results of the new methods were compared with methods reported in [11,19]. MATLAB 2009a platform was used for computation of numerical results.  Table 1. Figure 2a shows that the asymptotic order of convergence of the proposed method reached to O(ω −3.5 ), while the asymptotic order of convergence of the methods reported in [11] is O(ω −2 ) as shown in Figure 2b. The proposed method was tested for higher frequencies. The absolute errors are shown in Table 1. It is shown in the table that the proposed method Q C−L [r] improved accuracy on increasing ω.
Integral (22) was computed by the proposed methods Q C−L [r] for fixed frequency and varying nodal points and the results were compared with those obtained by The results are shown in Figure 3a, while in Figure 3b, the nodal points are fixed and the frequencies vary. It is shown in the figure that the proposed method was better than all the other methods. The main advantage of the proposed method Q C−L [r] was that it improved accuracy on increasing ω as well as N. Levin method (Middle and Bottom) [11] for test problem 1. Table 1. L abs produced by the Q C−L [r], (n, N = 10) for test Problem 1.

Example 2. Consider the following integral [19]
Oscillatory behavior of the real part of the integrand (23) is shown in Figure 4a. The integral has a stationary point of order k = 10 at x = 0, which lies at the extreme left position of the domain. The Chebyshev-Levin quadrature Q C−L [r] failed to evaluate the integral due to the existence of stationary point in the domain interval. Multi-resolution quadrature Q 8 nodes to give the desired accuracy for high frequencies, which was impractical. The integral was computed by the new splitting procedure ChQ accurately. Results in the form of scaled absolute errors and relative errors were computed and are shown in Figure 5. In Figure 5a, it has been shown that the asymptotic order of convergence of the splitting method ChQ reached to O(ω −3 ). The new splitting method improved accuracy on increasing nodal points as shown in Figure 5b. The two splitting methods ChQ and Chebyshev-Haar quadrature (CHQ) were implemented for fixed nodes and varying frequencies. The results are shown in Figure 6a. The new methods were compared with the result of the meshless splitting method [16] in Figure 6b. From the figures, it is shown that the new splitting methods were better than the methods reported in [16]. From all figures, it is evident that the new splitting methods were more accurate than the existing methods even for smaller nodes.  Example 3. Consider the following integral [11] The integral in (24) is highly oscillatory. Irregular high oscillations of real part of the integrand are shown in Figure 4b. The integral has a stationary point of order 2 at x = 0 ∈ [−1, 1]. The integral was split according to the case-II of the splitting methods and was computed by the two new splitting methods ChQ and CHQ. Absolute errors and scaled absolute errors are analyzed in Table 2 and Figure 7. Figure 7a indicates that the new splitting method ChQ retained the same asymptotic order of convergence O(ω −3 ) for this problem as well. In Figure 7b, results of the new splitting methods were compared with multi-resolution quadratures for fixed frequency and varying nodal points. The new splitting methods reduced the error up to O(10 −15 ), which was a higher rate of convergence than the method reported in [11].
Method ChQ was tested for higher frequencies. The results are shown in Table 2. It is shown in the table that the new methods improved the results on increasing ω as well as N.
The integral in (25) has a stationary point at x = 0. The integral was computed by the new splitting algorithms ChQ and CHQ. Results are analyzed in Figures 8 and 9. The result in terms of scaled absolute errors is shown in Figure 8a. In the same figure, results of the proposed methods are compared with the results of Filon-type method [11]. It has been shown that the asymptotic order of convergence of the new method ChQ reached to O(ω −3 ), while the method reported in [11] was of order O(ω −2 ).
Comparison of the proposed methods was performed for fixed ω and varying nodal points and vice versa as shown in Figure 9. It has been shown that the proposed splitting algorithms ChQ and CHQ improved the accuracy on increasing ω and N.

Conclusions
An adaptive splitting procedure which couples Chebyshev-Levin quadrature Q C−L [r] with multi-resolution quadratures was implemented to evaluate highly oscillatory integrals with a stationary point. Theoretical error analysis in terms of inverse powers of ω of the proposed procedure was performed and numerically verified by solving a few benchmark problems. The asymptotic order of convergence of both, the Chebyshev-Levin quadrature Q C−L [r] and the splitting algorithm reached to O(ω −3 ). Results of the proposed algorithm were compared with some existing methods reported in the literature. Moreover, the improvement in the accuracy could easily be witnessed from the results of new algorithm.