Efficient integration for a class of highly oscillatory integrals

Abstract

This paper presents some quadrature methods for a class of highly oscillatory integrals whose integrands may have singularities at the two endpoints of the interval. One is a Filon-type method based on the asymptotic expansion. The other is a Clenshaw–Curtis–Filon-type method which is based on a special Hermite interpolation polynomial and can be evaluated efficiently in O(N log N) operations, where N + 1 is the number of Clenshaw–Curtis points in the interval of integration. In addition, we derive the corresponding error bound in inverse powers of the frequency ω for the Clenshaw–Curtis–Filon-type method for the class of highly oscillatory integrals. The efficiency and the validity of these methods are testified by both the numerical experiments and the theoretical results.

Introduction

In many areas of science and engineering one often encounters the problem of computing rapidly oscillatory integrals. In most of the cases, such integrals cannot be calculated analytically and one has to resort to numerical methods. However, the numerical evaluation may be very difficult when the parameter ω is large, because in such case the integrand is highly oscillatory. A prohibitively large number of quadrature points is needed if one uses a classic rule such as Gaussian quadrature, or any quadrature method based on (piecewise) polynomial interpolation of the integrand. In this paper we are concerned with the efficient evaluation of highly oscillatory integrals of the type

$$I[f]={\int }_a^b(x-a{)}^{\alpha }(b-x{)}^{\beta }f(x)e^{i\omega x}\mathit{dx}\text{,}$$

where f(x) is a sufficiently smooth function on [a, b], ω is a large parameter, a and b are real and finite, and α > −1, β > −1. If α ⩾ 0, −1 < β < 0 or β ⩾ 0, −1 < α < 0, then the integrand in (1.1) has only one singularity at the endpoint b or a; if −1 < α < 0 and −1 < β < 0, the integrand has singularities at the two endpoints of the interval. Moreover, when α, β ⩾ 0, the integral (1.1) becomes a standard Fourier-type integral of the form

$${\int }_a^{b}F(x)e^{i\omega x}\mathit{dx}\text{,}$$

whose integrand has no singularity, and where F(x) =  (x − a)^α(b − x)^βf(x). A great many methods have been developed for Fourier transformation (1.2), such as Filon method [9], [14], [29], Clenshaw–Curtis-type method [7], Filon-type method [21], [22], asymptotic method [22], Levin method [26], generalized quadrature rule [6], [13], Levin-type method [35], numerical steepest descent method [20], [42], [23] and so on. However, once α and β are non-negative, the function F is non-singular, but it is smooth at the end-points (a condition for the analysis of methods like Filon or Clenshaw–Curtis) only if these two parameters are in addition integer. For the integral (1.1), 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.2) whose integrand does not involve singularity. These singularities mentioned above may impact heavily on its quadrature and its error. But since all these methods proposed in this paper, such as Filon-type method and Clenshaw–Curtis–Filon-type method, are based on the interpolation polynomials of smooth function f, in general, the required accuracy can be obtained by using derivatives at the end-points or adding the number of interior node points. Such integrals (1.1) with the weak singularities [24] are applied to the numerical approximations of solutions to Volterra integral equations of the first kind [5], [4]. In addition, the numerical integration of such integral (1.1) is employed to the solution of the singular integral equation for classical crack problems in plane and antiplane elasticity [18]. Therefore, it is of great importance for the study of the numerical integration of such integrals.

In fact, the integral (1.1) has gained popularity in literatures. In 1955, Erdélyi [11] established asymptotic expansions using neutralizer functions by van der Corput and general integration by parts. In 1956, these types of asymptotic expansions were presented in the treatise [12]. If f(x) is N times continuously differentiable for a ⩽ x ⩽ b, and −1 < α ⩽ 0, −1 < β ⩽ 0, then

$${\int }_a^b(x-a{)}^{\alpha }(b-x{)}^{\beta }f(x)e^{i\omega x}\mathit{dx}=B_N(\omega )-A_N(\omega )+O({\omega }^{-N-1})\text{,}\text{as}\omega \to \infty \text{,}$$

where

$$A_N(\omega )=\sum _{n=1}^N\frac{\Gamma (n+\alpha -1)}{(n-1)!}{e}^{\pi i(n+\alpha -2)/2}{\omega }^{-n-\alpha }{e}^{i\omega a}\frac{d^{n-1}}{{\mathit{dx}}^{n-1}}{\left[(b-x{)}^{\beta }f(x)\right]}_{x=a}\text{,}$$

$$B_N(\omega )=\sum _{n=1}^N\frac{\Gamma (n+\beta -1)}{(n-1)!}{e}^{\pi i(n-\beta -2)/2}{\omega }^{-n-\beta }{e}^{i\omega b}\frac{d^{n-1}}{{\mathit{dx}}^{n-1}}{\left[(x-a{)}^{\alpha }f(x)\right]}_{x=b}\text{.}$$

In 1958, a similar asymptotic result was re-established by Lighthill [27], by generalized function theory. In 1971, in the case that f(x) is analytic in a region containing [a, b], a straightforward proof based on contour integration was given by Lyness [30]. Recently, in 2008 and 2009, Lyness [31], [32] presented asymptotic expansions by theory involving inverse functions. For details one can refer to [31], [32]. However, little research has been done on the numerical computation of the integral (1.1). Based on these relevant background literatures above and the asymptotic expansion (1.3) at hand, in Section 2 of this paper, we will consider both the asymptotic method and the Filon-type method.

In many cases the accuracy of the Filon-type method is significantly higher than those of other methods, such as Filon method and asymptotic method. The computation of the coefficients of interpolation polynomial at the equally-spaced points needs to solve a linear system with O(N³) operations [44]. The main approximative power of the Filon-type method comes from matching function values and derivatives at the end-points (or, with greater generality, at critical points), and the addition of extra nodes is intended merely to decrease the error ‘constant’, rather than improving the asymptotic order. Indeed, just computing one extra derivative at the end-points is more effective asymptotically than any number of interior nodes. Therefore, for a well-behaved function f, it is nonsensical to use a large number of internal quadrature nodes, and this means that ill conditioning of the interpolation matrix will never occur in practice. Thus, in this case the Filon-type method is not prone to the Runge’s phenomenon in practical applications [36]. However, for the Runge’s function $f(x)=\frac{1}{1+25x^2}$ [38], now over the interval [−1, 1], it suffers from Runge’s phenomenon, where certain nonoscillatory functions f have oscillating interpolation polynomials at the equally-spaced points. Moreover, similar to Runge’s phenomenon is the situation when f and its derivatives increase too fast to be accurately approximated at the equally spaced points by polynomials, the use of the Filon-type method may not result in better approximations. Since the Filon-type method is based on interpolation, it is logical that the accuracy of the Filon-type method is directly related to the interpolation accuracy. Thus, in these cases for the complicated f, the accuracy of the Filon-type method is limited. For the complicated f, to avoid the Runge’s phenomenon or the similar situation, the use of Clenshaw–Curtis points, if chosen as interpolation points, is a better choice than that of the equally spaced points [33], [38], [41]. Meanwhile, adding additional node points results in a more accurate approximation [36]. Moreover, the Chebyshev coefficients can be evaluated by using an efficient and numerically stable algorithm in O(N log N) operations, based on FFT [8], [9], [15], [41], which avoids solving a linear system with O(N³) operations. The better choice of interpolation points can be made far more efficient and stable from a computational point of view than the equally-spaced points and it ensures uniform convergence for most of continuous functions [33]. The authors of [44] show that the Clenshaw–Curtis method [7], [3], [37] is much more efficient than the Filon method [14], [22], [43] with the equispaced points or Gauss–Lobatto points on [a, b]. Therefore, based on these aforementioned advantages of Clenshaw–Curtis points, in Section 3 of this paper, our attention focuses on adapting the Clenshaw–Curtis points for the computation of (1.1), which will induce a Clenshaw–Curtis–Filon-type method,

$$Q_s^{C''–''C}[f]={\int }_a^b{H}_{N+2s}(x)(x-a{)}^{\alpha }(b-x{)}^{\beta }{e}^{i\omega x}\mathit{dx}\text{,}$$

where H_N+2s(x) is the special Hermite interpolation polynomial of f(x) at the Clenshaw–Curtis points [2] $c_n=\frac{b+a}{2}+\frac{b-a}{2}\cos \left(\frac{n\pi }{N}\right)$ on [a, b] satisfying

$$H_{N+2s}^{(k)}(a)=f^{(k)}(a)\text{,}{H}_{N+2s}(c_n)=f(c_n)\text{,}{H}_{N+2s}^{(k)}(b)=f^{(k)}(b)\text{,}$$

for n = 1, 2, … , N − 1 and k = 0, 1, … , s; it can be evaluated in O(N log N) operations [44]. The choice of the extreme points as interpolation points for highly oscillatory integrals is not just a technical necessity, still it can improve the accuracy and asymptotic order significantly [22], [35], [43], [44] (or see (2.11), (3.20), (3.33) below). Meanwhile, the Clenshaw–Curtis–Filon-type method, owing to the easy implementation and its good stability and convergence properties in N for a fixed ω or in ω for a fixed N, comes into its own when we need to compute

$${\int }_a^{b}f(x\text{,}y)e^{i\omega x}\mathit{dx}\text{,}$$

say (it is also possible to add singularities at the end-points) for many values of y along a grid, a situation often encountered when solving integral equations, e.g. in high-frequency scattering in [10] and the references therein, which is our main target application.

An outline of this paper is organized as follows. In Section 2, we describe both the asymptotic method and Filon-type method for (1.1) and present explicit expressions of the generalized moments for (1.2). A simply computation of Kummer confluent hypergeometric functions are also discussed there. In Section 3, we give a Clenshaw–Curtis–Filon-type method for (1.1) and present a five-term recurrence relation of the modified moments for (1.1). Finally, both the corresponding error estimate and the uniform convergence for the Clenshaw–Curtis–Filon-type method are illustrated. Meanwhile, in Sections 2 Asymptotic method and Filon-type method, 3 Clenshaw–Curtis–Filon-type method, some preliminary robust numerical examples show the accuracy and effectiveness of this set of ideas. These vigorous results mentioned below indicate that the weak singularity plays a role fairly similar to stationary points in classical asymptotic analysis of highly oscillatory integrals (see (2.10), (2.11), (3.20), (3.33)).