Efﬁcient Computation of Highly Oscillatory Fourier Transforms with Nearly Singular Amplitudes over Rectangle Domains

: In this paper, we consider fast and high-order algorithms for calculation of highly oscillatory and nearly singular integrals. Based on operators with regard to Chebyshev polynomials, we propose a class of spectral efﬁcient Levin quadrature for oscillatory integrals over rectangle domains, and give detailed convergence analysis. Furthermore, with the help of adaptive mesh reﬁnement, we are able to develop an efﬁcient algorithm to compute highly oscillatory and nearly singular integrals. In contrast to existing methods, approximations derived from the new approach do not suffer from high oscillatory and singularity. Finally, several numerical experiments are included to illustrate the performance of given quadrature rules.


Introduction
Highly oscillatory integrals frequently arise in acoustic scattering [1], computational physical optics [2], computational electromagnetics [3], and related fields. Generally, dramatically changing integrands make classical approximations perform poor. Therefore, studies on numerical calculation of highly oscillatory integrals have attracted much attention during the past few decades, and a variety of contributions has been made, for example, Filon-type quadrature [4,5], numerical steepest descent method [6], Levin method [7], and so on.
When the phase is nonlinear, researchers usually resort to Levin-type quadrature, which originates from David Levin's pioneering work in [7]. By transforming the oscillatory integration problem into a special ordinary differential equation, one could get an efficient approximation to the generalized Fourier transform with the help of collocation methods. Afterwards, Levin analyzed the convergence rate of the innovative approach in [8]. Analogous to Filon-type quadrature, the Levin-type method based on Hermite interpolation was developed by Olver in [9]. Application of Hermite interpolation definitely increased the convergence rate of the numerical method with respect to the frequency. In [10], Li et al. proposed a stable and high-order Levin quadrature rule by employing the spectral Chebyshev collocation method and truncated singular value decomposition technique. Multiquadric radial basis functions were applied to Levin's equation and an innovative composite Levin method was presented in [11]. Numerical tests manifested that such kind of algorithms was able to deal with stationary problems. Sparse solvers for Levin's equation in one-dimension were constructed by employing recurrence of Chebyshev polynomials in [12,13]. Meanwhile, a class of preconditioners was proposed by the second author to deal with the ill-conditioned linear system in [13]. Molabahrami studied the Galerkin method for Levin's equation and developed the Galerkin-Levin method for oscillatory integrals in [14]. F(x, y) (x − a) 2 + (x − b) 2 + 2 e iω(G 1 (x)+G 2 (y)) dxdy. (1) In this paper, we are concerned with efficient computation of Integral (1), and partly fill in the gap in this field. Moreover, we suppose that F(x, y) is analytic with respect to both variables, G 1 (x) and G 2 (y) are sufficiently smooth functions without stationary points, and the frequency parameter ω 1, (a, b) ∈ R 2 , and | | 1. A large frequency parameter ω implies that integrands of Integral (1) are highly oscillatory, and classical quadrature rules suffer from the computational cost. In Table 1, we list numerical results computed by the classical delaminating quadrature rule coupled with Clenshaw-Curtis quadrature (CCQ), where the quadrature nodes are fixed 16. Referenced values are computed by CHEBFUN toolbox (see [21]). CHEBFUN, which approximate functions by Chebyshev interplant, was firstly developed in 2004. Due to the fast and high-order approximation to the integrand, numerical integration methods in CHEBFUN usually provide efficient numerical approaches for univariate and multivariate integrals. Hence, the 2D quadrature method in CHEBFUN is chosen as a benchmark. It can be seen from Table 1 that, as ω goes to infinity, CCQ diverges from the referenced values when we do not add quadrature nodes. In contrast to oscillatory integrals arising in existing studies, when the point (a, b) in Integral (1) is close to or falls in the integration domain and is particularly small, the integrand attains its peak value around (a, b) and decays dramatically away from such a critical point. In general, the point (a, b) is called the nearly singular point. Plenty of additional quadrature nodes have to be used if we want to make the numerical formula retain a tolerance error.
There also exist several contributions to tackle the nearly singular problem. The sinh transformation is deemed one of the most important tools. For nearly singular moments arising in Laplace's equation, Johnston et al. proposed the sinh transformation in [22]. Occorsio and Serafini considered two kinds of cubature rules for nearly singular and highly oscillatory integrals in [23]. With the help of 2D-dilation technique, Occorsio and Serafini were able to relax the fast changing integrand and applied Gauss-Jacobi quadrature to the transformed integral. Numerical experiments verified that such an approximation procedure greatly increased the numerical performance of Gauss quadrature.
The remaining parts are organized as follows. In the second section, we review some results with regard to the calculation of Chebyshev series and present the convergence property of Chebyshev interplant and series. In Section 3, we first extend the idea in [13] to two-dimensional oscillatory integrals. Compared with existing Levin quadrature, the new approach has an advantage in computational time. Then, noting that there is little convergence analysis of 2D Levin quadrature rules, we try to fill the gap through examining the modified Levin equation. Finally, we present an innovative composite Levin quadrature rule for solving nearly singular and highly oscillatory problems. In contrast to existing numerical integration methods, the proposed composite method does not suffer from high oscillation and nearly singular amplitudes. Numerical tests included in Section 4 are conducted to verify the efficiency of the proposed approach, and some remarks are concluded in Section 5.

Auxillary Tools
In this section, we first revisit auxillary operators with regard to Chebyshev series, which help develop numerical algorithms for computation of two-dimensional oscillatory integrals. Then, error bounds for coefficients of Chebyshev series and Clenshaw-Curtis interplant are introduced.
When the given function f (x) is analytic in a sufficiently large domain containing [−1, 1], one can compute its Chebyshev series by (see [24]) Here, T n (x) denotes the first-kind Chebyshev polynomial of order n. Noting the relation between the first-and second-kind Chebyshev polynomials T n (x), U n (x) (see [24]) and we are able to compute which implies Chebyshev coefficients of the derivative can be represented by Secondly, suppose that there exists a sufficiently smooth function a(x) = ∞ ∑ n=0 a n T n (x).
Noting the identity (see [24]) where a = [a 0 , a 1 , · · · ] T , and coefficients {c n } ∞ n=0 are defined by Operators D, M[a] have been verified to be efficient tools for discretizing Levin's equation. For more details, one can refer to [13]. On the other hand, when f (x) is analytic with | f (x)| ≤ M in the region bounded by the Bernstein ellipse with the radius ρ > 1, we have for every n = 0 (see [25]) Noting that Furthermore, according to [ and Here, p N (x) denotes the interplant of f (x) at Clenshaw-Curtis nodes or the truncated Chebyshev series of f (x).

Main Results
This section is devoted to investigating fast algorithms for calculation of Integral (1). To begin with, let us consider the computation of oscillatory integral without nearly singular integrands, that is, Here, F(x, y), G 1 (x), G 2 (y) are smooth functions with sufficiently large analytic regions, and G 1 (x), G 2 (y) do not have stationary points in the complex plane. Consider the inner integral For fixed y j = cos j N π, j = 0, 1, · · · , N, we are restricted to finding a function P j (x) satisfying Noting that G 1 (x) never vanishes over the interval [−1, 1], we can get the modified Levin equation, Let With the help of operators D, M[a], we rewrite modified Levin's Equation (11) as where Solving Equation (12) by the truncation method [26] gives the unknown coefficients p j n , n = 0, 1, · · · , N, and we can get approximations to P j (±1) by Clenshaw algorithm, where p j,N k denotes the approximation to p j k . Hence, the inner integral (9) is computed by Since H N (y j ), the approximation to H(y) at Clenshaw-Curtis nodes, has been obtained, we are able to construct the polynomial H N (y) by where L j (y) denotes Lagrange basis with respect to Clenshaw-Curtis nodes and h n can be computed by fast Fourier transform. Letting where we are able to approximate q 0 , · · · , q N by the truncation method again. Computing a 0 (±1), a 2 (±1) by where q N k denotes the approximation to q k , we arrive at 2D spectral coefficient Levin quadrature for Integral (8) In [27], Xiang established the relation between Filon and Levin quadrature rules in the case of the phase g(x) = 1 and analyzed the convergence property of Levin quadrature. Instead of resorting to Filon quadrature, we consider the convergence rate of the above spectral coefficient Levin method with respect to quadrature nodes and frequency in the case of nonlinear oscillators through examining the decaying rate of the coefficients.
For any fixed y ∈ [−1, 1], F(x, y) turns to the univariate function with regard to x.
we summarize the convergence analysis in the following theorem.
Then, for sufficiently large ω, we have where the constant C does not depend on ω, N.
calculation implies the quadrature error can be decomposed intô Here, N (y)e iωG 2 (y) dy, In the remaining work, we give estimates for E 1 , E 2 , E 3 with respect to the increasing truncation term N and frequency ω. For E 1 , note that H(y) is bounded within its Bernstein ellipse by As a result, we have according to integration by parts where we obtain Furthermore, letting we have A direct calculation as is done in the estimation procedure for E 1 results in Then, let us consider the decaying rate of coefficients of Chebyshev expansions of Err 2,j (x), which helps to analyze Err 2,j (x) ∞ and Err 2,j (x) ∞ . In fact, the truncation technique implies with and p j,N k denotes the approximation to p j k in Equation (12). For sufficiently large ω > N, it follows that 1 Denoting the maximum of The Chebyshev coefficients of Err 2,j (x) can be computed by Noting the construction technique in the modified spectral Levin coefficient method, we get c 2,j n = 0 for n = 0, 1, · · · , N. On the other hand, for n ≥ N + 1, it follows that It is noted that the first term in the right-hand side of the above equation is the coefficient of F(x, y j ) G 1 (x) and the second term is that of P j (x) , where we denote coefficients to be c FG n , c PG n , respectively. Recalling the product operator in Equation (5), we have Therefore, it follows |c 2,j Here, C := 2 max{6M G 1 M F , 12M G 1 M F S N }. Hence, Err 2,j (x) ∞ and Err 2,j (x) ∞ can be bounded by and As a result, it follows that Now, we arrive at the fact where C 2 := 64C π .
The estimation procedure for E 3 is similar to that of E 2,j . We ignore details and give the conclusion directly where C 3 does not depend on N and ω.
To sum up, we arrive at the following error bound by combining Equations (17), (22), and (23), with C = max{C 1 , C 2 , C 3 }. It is easily seen that the constant C does not depend on N and ω. This completes the proof.
Finally, let us turn to the construction of the composite quadrature rule for calculation of Integral (1). It is observed in the above theorem that, when the radiuses ρ F , ρ H are close to 1, the error bound would expand dramatically. Therefore, an efficient quadrature rule has to guarantee the fact that the integrand has a relatively large analytic radius over the integration domain. To make this judgment be satisfied, we choose a non-uniform grid instead of partitioning the integration region uniformly.
To begin with, the singular point z * is projected into the plane containing the integration region and we get the projection point z. In the case of the projected point z falling into the integration domain (Case I), the first box is determined by the distance between z * and z. We construct a square with its center being z and its side length being 2 z * − z . Then, the side length of level-2 box's with the center z is set to be 2 2 z * − z . To devise the composite quadrature rule, we first select level-1 box as a subdomain. Noting that the remaining domain is not a rectangle, we partition it into four subdomains, that is, Box21, Box22, Box23, and Box24 (see Figure 1). In general, the side length of level-l box's with the center z is set to be 2 l z * − z , and the integration subdomain is constructed similarly, which finally results in a nonuniform grid (see Figure 2).
When the projected point falls out of the integration domain, for example, it is around the side (Case II) or vertex (Case III), we implement a similar partition procedure like that in Case I. The final partition grid is shown in Figures 3 and 4. Applying 2D spectral coefficient Levin quadrature rule in the subdomain leads to a class of composite 2D spectral coefficient Levin quadrature. It is noted that such kind of partition techniques guarantee the fact that the distance between the singular and the integration interval is no less than 2 when we map the integration domain into [−1, 1] × [−1, 1].

Numerical Experiments
This section is devoted to illustrating the numerical performance of 2D spectral coefficient Levin quadrature (2DSC-Levin) and composite 2D spectral coefficient Levin quadrature (C2DSC-Levin) given in Section 3. cos(x + y)e iω(x+y) dxdy.

Example 1. Let us consider the computation of the oscillatory integral
The phase x + y has no stationary points within the domain [−1, 1] × [−1, 1], and the amplitude cos(x + y) is analytic with respect to both variables. Therefore, it is expected that the new approach has an exponential convergence rate.
It is noted that approximation results derived from classical cubature usually do not make sense when ω 1. We first employ CCQ and give the computational results in Table 2, where both quantity of quadrature nodes N and oscillation parameter ω are variables. Although plenty of quadrature nodes have been used in the above example, computed results are not satisfactory especially in the case of high oscillation. Now, we list approximated results of 2DSC-Levin in Table 3, where the referenced exact value is computed by the CHEBFUN toolbox again. It can be seen from Table 3 that absolute errors do not increase as the frequency ω enlarges, which implies that the new method is robust to high oscillation. On the other hand, when we raise the truncation term of Chebyshev series, the absolute error decays fast. In Figure 5, we give tendencies of absolute errors with respect to increasing frequencies and compare the computational time of 2DSC-Levin and referenced algorithm in CHEBFUN. It can be found that the consumed time of 2DSC-Levin does not vary as the frequency ω enlarges, whereas that of CHEBFUN's 2D quadrature (2D-Cheb) increases dramatically. Since 2D-Cheb is a class of self-adaptive algorithms, it has to increase quadrature nodes when the frequency goes to infinity to retain a tolerance error, which results in the dramatically growing curve. However, approximations derived from 2DSC-Levin do not suffer from high oscillation according to Theorem 1. Hence, there is no need to raise quadrature nodes of 2DSC-Levin in high oscillation, and we do not witness an obvious change in Figure 5.  We also employ 2D-Cheb-Levin quadrature given in [18] to give a comparison. Computed results are shown in Table 4.  Tables 3 and 4 illustrates that the accuracy of 2DSC-Levin and 2D-Cheb-Levin quadrature is similar. However, 2DSC-Levin does a little better than 2D-Cheb-Levin quadrature when CPU time is considered. Since both approaches consist of approximations to a series of one-dimensional integrals, we show the consuming time of both approaches for computing the final highly oscillatory integrals in Table 5, where it can be seen that 2DSC-Levin is slightly faster than 2D-Cheb-Levin quadrature, which is partly due to the sparse structure of the discretizatized modified Levin equation. Table 5. Comparison of CPU time for 2DSC-Levin and 2D-Cheb-Levin quadrature for fixed ω = 1000. x 2 + y 2 + 15 e iω(x 2 +x+y 2 +y) dxdy.

2DSC-Levin
The amplitude 1 is no longer an entire function, and the inverse of the phase function x 2 + x + y 2 + y can not be calculated directly.
We show absolute errors and CPU time of 2DSC-Levin in Table 6 and Figure 6, respectively. Due to the fact that the amplitude in Example 2 has a limited analytic radius, 2DSC-Levin converges to the machine precision much more slowly than that in Example 1. However, noting the decaying curve in the left part of Figure 6 manifests that 2DSC-Levin has the property that the higher the oscillation, the better the approximation, which also coincides with the theoretical estimate in Theorem 1. Hence, 2DSC-Levin is feasible for calculation of highly oscillatory integrals over rectangle regions when the oscillator g(x, y) is nonlinear. In addition, it is interesting that the curve of CPU time of 2DSC-Levin has a jump at about ω = 5500 in the right part of Figure 6. Such a phenomenon may originate from the fact that the Levin equation can be solved more efficiently as the frequency becomes larger. However, this is still a conjecture and we need more theoretical investigation in the future work.  To illustrate the effectiveness of the composite 2D spectral coefficient Levin quadrature rule (C2DSC Levin), we give a comparison among the new approach, 2D sinh transformation (JJE) in [22], and 2D dilation quadrature (2D-d) in [23].
For Integral (1), the sinh transformation is defined by Since the transformed integrand is no longer nearly singular, a direct 2D Gauss cubature can be applied in practical computation.However, it should be noted that, although JJE can efficiently deal with nearly singular problems, it generally suffers from the highly oscillatory integrands.
In [23], Occorsio and Serafini proposed 2D-d for the integral where D : Letting Properly choosing d ∈ R + and S = 2ω 1 d ∈ N results in Here, D i,j : Employing the transformed Gauss-Jacobi quadrature to the moment integral gives 2D-dilation quadrature.

Example 3. Consider the computation of
It is noted that the integrand will reach its peak value at (−0.5, 0.5), and dramatically decrease away from such a critical point.
We list computed absolute errors of JJE, 2D-d and C2DSC-Levin in Tables 7-9. As the number of quadrature nodes N increases, all of algorithms converge fast to the referenced value. When the integrand does not change rapidly, JJE provides the best numerical approximation. However, as the frequency ω enlarges, JJE and 2D-d suffer from high oscillation, while C2DSC-Levin is still able to maintain a relatively high-order approximation. It should be noted that the dilation parameter d in 2D-d is restricted to make the number of quadrature nodes coincide with the other two methods, and a slightly modified choice as is done in [23] may make 2D-d be able to deal with some highly oscillatory problems. The corresponding results are shown in Figure 7. Numerical results in this figure indicate that the absolute error derived from JJE increases dramatically when ω goes beyond 500, while the error of 2D-d rises slowly. It is also found that both absolute errors and computational time of the new approach do not suffer from the varying frequency ω. Hence, C2DSC-Levin is the most effective tool for computing oscillatory and nearly singular integrals.    Although JJE is efficient for solving some nearly singular problems, it may fail when b = 0 in the integrand.  Tables 10 and 11, respectively. It can be found that 2D-d provides more accurate approximation than that of C2DSC-Levin in the relatively low frequency. Nevertheless, the absolute error of C2DSC-Levin never increases in the high frequency while its 2D-d counterpart enlarges. Furthermore, although 2D-d will provide a more accurate approximation if we employ the choice of the dilation parameter considered in [23], it still cannot beat C2DSC-Levin in the high frequency when both absolute errors and computational time are considered (see Figure 8).

Conclusions
In this paper, we have presented the modified spectral coefficient Levin quadrature for calculation of highly oscillatory integrals over rectangle regions and established its convergence rate with respect to the truncation term and oscillation parameter. Furthermore, by considering numerical calculation of moments over a non-uniform mesh, we derive the composite Levin quadrature. Numerical experiments indicate that the non-uniform partition technique greatly reduces the nearly singular problem. Recently, sharp bounds for coefficients of multivariate Gegenbauer expansion of analytic functions have been studied in [28], which definitely opens a door for our ongoing work about convergence analysis of Levin quadrature in high dimensional hypercube.
On the other hand, studies on the asymptotic and oscillatory behavior of solutions to highly oscillatory integral and differential equations have attracted much attention during the past decades [29][30][31][32]. It is noted that computation and numerical analysis of oscillatory integrals provide efficient tools for such kinds of studies and investigation of application of the proposed approaches to oscillatory equations is also necessary in the future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: