A Gaussian quadrature rule for oscillatory integrals on a bounded interval

We investigate a Gaussian quadrature rule and the corresponding orthogonal polynomials for the oscillatory weight function $e^{i\omega x}$ on the interval $[-1,1]$. We show that such a rule attains high asymptotic order, in the sense that the quadrature error quickly decreases as a function of the frequency $\omega$. However, accuracy is maintained for all values of $\omega$ and in particular the rule elegantly reduces to the classical Gauss-Legendre rule as $\omega \to 0$. The construction of such rules is briefly discussed, and though not all orthogonal polynomials exist, it is demonstrated numerically that rules with an even number of points are always well defined. We show that these rules are optimal both in terms of asymptotic order as well as in terms of polynomial order.

. Prompted by the unexpected succes of modified Magnus expansions for oscillatory differential equations [9], these authors introduced asymptotic and Filon-type methods as a generalisation of the classical Filon method [3]. Other innovations in this field include the numerical method of steepest descent [6], based on numerical evaluation of steepest descent integrals, and Levin-type methods [15,18,19], based on solving an associated differential equation. A common property of all these methods is that high asymptotic order can be attained, by which we mean an error that behaves like O(ω −α ), ω → ∞, for a positive α. Unlike classical asymptotic expansions, whose error behaves similarly as ω → ∞, these methods are in principle not limited in accuracy for a fixed ω.
This work concerns an asymptotically optimal oscillatory quadrature rule that is valid for all frequencies. The discussion will focus on the Fourierintegral where we shall consider ω ranging from 0 to ∞. Of the above mentioned quadrature rules, the numerical method of steepest descent can be considered to be asymptotically optimal; it delivers an error of size O(ω −2n−1 ), when using 2n quadrature points in the complex plane for this integral. Filon-type methods are based on interpolation of the amplitude function f . By choosing interpolation points that scale towards the endpoints like 1/ω, they deliver O(ω −n−1 ) error with the same number of points [10]. In [5] it was pointed out that by choosing the interpolation points for the Filon-type method to be precisely the quadrature points in the complex plane of the numerical steepest descent method, the so-called superinterpolation points, also here an error of O(ω −2n−1 ) can actually be attained with only 2n points. Next, one could ask, how well these methods fare when ω is small. Unfortunately, the superinterpolation points are unbounded in the limit ω → 0. As such, they are not a good choice for small ω. For ω = 0 the Gauss-Legendre points are optimal in the sense that they integrate exactly polynomials of degree 2n − 1 using n points. As for Filon-type methods, the interpolation points can be chosen such that they approach the Legendre points as ω → 0 [10]. However, it is not clear what is an optimal way to do so. In the context of the so-called exponential fitting methods, such strategies have been described with relatively few quadrature points [13], or with an heuristic dependence on ω [14]. However, in both cases, using 2n points still yields an error of size O(ω −n−1 ), which is not optimal.
In a recent work, [1], it was shown that Gaussian quadrature rules can be constructed for certain oscillatory integral transforms, i.e., integrals over unbounded domains. This goes against common wisdom, saying that Gaussian rules can be found for positive weight functions only. A more correct statement is that existence and uniqueness proofs, as well as many construction methods, rely on positive weight functions. For integral transforms of the form it was shown in [1] that applying the Gauss-Laguerre rule along the path of steepest descent yields a Gaussian rule, which is exact for f being a polynomial of degree ≤ 2n − 1. Moreover, a clear connection was shown between polynomial accuracy and asymptotic accuracy for such integrals, where higher polynomial accuracy means higher asymptotic accuracy. Thus, the Gaussian rule for this integral transform attains an error O(ω −2n−1 ). Similar results were shown for other oscillators, all resulting in rules with nodes in the complex plane.
The topic of the current work is a further investigation of Gaussian rules for oscillatory weight functions, but now on a bounded domain. This, we will see, is a more difficult case, since it is no longer possible to remove the dependency on ω by a simple rescaling as in the unbounded case. For integral (1), the oscillatory part e iωx is the weight function, and one thus seeks rules that integrate polynomials up to degree 2n − 1 exactly. This endeavour poses several problems. The rules depend on ω in a non-trivial way, and must be computed numerically. One cannot guarantee existence of the rules, though numerical evidence indicates that all rules with an even number of points exist. On the other hand, the resulting rules reduce elegantly to Gauss-Legendre by construction in the limit ω → 0. Moreover, it will be proved, under mild assumptions, that the quadrature points tend to the superinterpolation points in the high-frequency limit ω → ∞. This implies that the method yields optimal asymptotic order in addition to optimal polynomial order. These observations are the basis of the statement that this Gaussian rule is truly optimal for oscillatory integrals of the form (1) throughout the frequency regime.
This paper is built up as follows. In §1, we recall some preliminary facts regarding Gaussian quadrature and highly oscillatory quadrature methods. In §2, analytic expressions for the cases n = 1 and n = 2, as well as several numerical experiments, shed light on the questions of existence and other properties of Gaussian rules for integrals of the form (1). This section is concluded with a set of conjectures, most of which are proved in §3, on the properties of the orthogonal polynomials, and §4, on the properties of the quadrature rule.

Gaussian quadrature
An n-point Gaussian quadrature rule {x j , w j } n j=1 is a quadrature rule which has optimal polynomial accuracy, i.e., it is exact for all polynomials of degree It is well known that the nodes of a Gaussian rule are precisely the zeros of the n-th orthogonal polynomial with respect to the weight function h [4]. Classical theory of Gaussian quadrature ensures that such rules exist and that they are uniquely defined whenever h is positive. Moreover, the positivity of the weight function also ensures that x j ∈ [a, b], j = 1, . . . , n, and that the weights w j are all positive. The monic orthogonal polynomial p n can be computed in various ways, for example by the recurrence with initial values p −1 = 0 and p 0 = 1. Defining the pairing which is an inner product if h is positive, the recurrence coefficients are given by Alternatively, writing p n (x) = x n + n−1 k=0 a k x k , the coefficients a 0 , . . . , a n−1 satisfy the linear system where the moments µ m are defined as

Gaussian quadrature and the method of steepest descent
In the method of steepest descent, an oscillatory integral of the form (1) is rewritten along the so-called paths of steepest descent. Assuming f is analytic, we may deform the path of integration as follows: (5) Applying the Gauss-Laguerre quadrature on the two resulting integrals yields an error that decays like O(ω −2n−1 ), when using n points for each integral (so 2n points in total) [6]. The Filon-type method, due to Iserles & Nørsett [8], is based on polynomial interpolation of the amplitude function f . In [5] it was shown that using the complex points obtained by applying Gaussian quadrature on the steepest descent integrals in (5) as interpolation points in a Filon-type method, also yields an error decay of O(ω −2n−1 ). Following the terminology of [5], these points are referred to as superinterpolation points: denote the n point Gauss-Laguerre quadrature rule. The corresponding 2n superinterpolation points are defined as .
The corresponding quadrature weights for these points are Note that Filon-type methods have the advantage that accuracy can be increased by adding interpolation points wherever they are needed, while maintaining asymptotic accuracy. This can be used to minimise the effect of eventual singularities in the complex plane. This leads, of course, to different weights. Also note that the superinterpolation points are not well defined in the limit ω → 0. and uniqueness theory. However, assuming existence for the time being, at least for some n and ω, we denote the n-th monic orthogonal polynomial by p ω n . The polynomial p ω n is orthogonal in the sense that, The moments can be computed explicitly, where Γ(z, a) is the incomplete Gamma functions [17], or by the recurrence which is obtained through integration by parts. This recurrence should be used with great care, since it's forward stable only when m < ω. For m > ω, the backward recursion should be used. Assuming existence of these polynomials, they can be computed using the linear system for the coefficients (4). For the cases n = 1 and n = 2 the polynomials can be computed analytically, and this gives some insight into the nature of the higher order rules.

The case n = 1
In the case n = 1 the orthogonal polynomial takes the form p ω Note that the single quadrature point −a 0 remains on the imaginary axis. The rule is undefined when ω is a multiple of π, except in the limit ω → 0 where a 0 → 0. This limit corresponds to the Gauss-Legendre rule.

The case n = 2
In the case n = 2, one obtains for p ω 2 (x) = x 2 + a 1 x + a 0 the coefficients  In the low-frequency limit, for ω → 0, we recover the Gauss-Legendre rule again. Observe that Since the roots They are a perturbation of the Gauss-Legendre nodes as expected.
For ω → ∞, we observe that . These turn out to be the two-point superinterpolation points as defined in Definition 1.1 -recall that 1 is the single root of the Laguerre polynomial of degree 1.
The roots of the polynomials, which are the two quadrature points for the Gaussian rule, are given explicitly by Although this expression is rather complicated, one sees that the points exist for all ω, unlike in the case n = 1. Fig. 1 shows the curves that the two quadrature points trace out in the complex plane as ω increases. The qualitative behaviour seen in this figure also appears for higher n. The points leave the real line and drift into the complex plane orthogonal to the real line. After an excursion into the complex plane the points attract in a rather irregular manner towards the two endpoints of the interval. The curves appear to have cusps for certain values of ω. Regarding the curves as parametric curves, it is clear that a cusp can only happen at a singular point, i.e., where ∂x ± ∂ω vanishes. In order to compute the recurrence coefficients, defined in (3), that would enable us to compute p ω 3 , we need the quantity  Figure 3: The absolute value of (p ω 2 , p ω 2 ) (continuous) and ∂x 1 ∂ω (dashed) as a function of ω.
For certain values of ω, this expression vanishes and the recurrence coefficient α 2 does not exist. The first such value occurs near ω = 2π, and this value can be numerically computed to be ω * 1 = 5.92966 . . .. The values of ω for which (p ω 2 , p ω 2 ) = 0 are difficult to compute explicitly, but one can give some information for large values of ω. Dividing by ω 4 throughout, we obtain that (p ω 2 , p ω 2 ) = 0 is equivalent to 2 cos(ω) ω Hence, to leading order we have the zeros of sin ω as solutions, so ω k = kπ + O(k −1 ), k → ∞. The first correction to this estimation gives This is illustrated in Figure 3, where also | ∂x + ∂ω | is plotted. The figure shows that the zeros of (p ω 2 , p ω 2 ) appear to coincide with the zeros of | ∂x + ∂ω |, and thus these values of ω correspond to the the cusps in Fig. 1. We may conclude from this that there is a breakdown of the recurrence (2) for countably many values of ω, and that these problematic values correspond to cusps in the curves x j (ω). As a result, p ω 3 is undefined for these special values of ω, just like p ω 1 was undefined for values of ω that are exact multiples of π. However, by avoiding the recurrence and by computing with the moments as in eq. (4), for example, one can still compute the orthogonal polynomial of degree 4 for all values of ω. The quadrature method obtained from the superinterpolation points has asymptotic order 3, and this method presumably has the same order. Applying it to the integral the error is shown in Fig. 4. We see that the resulting method indeed appears to have asymptotic order 3. Note that a Filon-type method with two quadrature points is in general only expected to have asymptotic order 2. This can be achieved by using the endpoints ±1, or any two points that move towards ±1 at a rate of 1/ω [10].

Numerical experiments for n > 2
For n > 2, obtaining closed forms expressions for zeros of the orthogonal polynomials is highly impractical, if at all possible. Turning to numerics, the coefficients of the monic polynomial can be computed from the Hankel system (4), and the roots can be found numerically. Note that this system is likely to be ill-conditioned. The computations in this paper have been carried out in high precision in Maple.
The roots for n = 3 behave in a similar manner as in the case n = 2. In Fig. 5, we see that the cusps in this case seem to coincide with the cusps in the case n = 2. However, there is an extra root which is always on the imaginary axis, a consequence of the symmetry of the polynomials  with respect to the imaginary axis. In Fig 7, one observes that this root is indeed undefined for a set of discrete values ω identified as the cusps in the case n = 2.
One can also compute zeros of polynomials of higher order. For n = 16, the zeros behave qualitatively in much the same way as the zeros of p ω 2 . The cusps in the curves correspond to the same critical values of ω regardless of the root, but these values are different from those in the case n = 2. The next experiment concerns the behaviour of the nodes for large ω. Figure 9 shows the difference between the 8 roots of p ω 16 near −1 and the 8 corresponding superinterpolation points. It appears that the difference goes like O(ω −2 ), similar to what was established for the case n = 2. For large values of ω, the roots of the orthogonal polynomial tend to the superinterpolation points.

Conjectures based on observations
From the above discussion we can list several conjectures: 1. Rules with an even number of points exist for all ω.
2. Rules with an odd number of points are not defined for all ω.
3. In particular, p ω 2n+1 is not defined for those critical values of ω that correspond to the cusps of the roots of p ω 2n .
4. The polynomials p ω n are symmetric with respect to the imaginary axis.  5. The roots of p ω 2n approach the 2n superinterpolation points at a rate of O(ω −2 ) in the high-frequency limit. 6. The Gaussian rule based on the zeros of p ω 2n has asymptotic order O(ω −2n−1 ).
In this paper we shall not prove the first conjecture, on the existence of evendegree polynomials, but assume it to be true. We shall prove conjectures 3 (Lemma 3.1) and 4 in §3 and conjectures 5 and 6 in §4.

Properties of the orthogonal polynomials
In this section we set out to describe a number of interesting properties of the orthogonal polynomials that are useful later on to explain the observations that were made in the previous section.

Symmetry
Symmetries of weight functions and intervals lead to symmetries of the corresponding polynomials. The complex exponential weight function has symmetries with respect to the imaginary axis, in the sense that w(z) = w(−z). Note that the point −z is the reflection of z with respect to the imaginary axis. This leads to the following. Lemma 3.1 Let w(z) be a weight function such that w(z) = w(−z) and Γ be a contour that is invariant under reflection with respect to the imaginary axis, i.e., z ∈ Γ ⇒ −z ∈ Γ.
Proof 1 Using the symmetry of Γ and of the weight function respectively, we find Thus, p n (−z) satisfies the same orthogonality conditions as p n (z). Since the latter is unique, we must have that for some constant c. Using that p n is monic yields c = (−1) n .
The weight function e iωx satisfies the required symmetry of Lemma 3.1 and so does the interval [−1, 1]. Therefore the polynomials p ω n satisfy the symmetry (7). This implies among other things that they are either purely real or purely imaginary along the imaginary axis, depending on the parity of n.

Derivatives with respect to ω
We gain useful knowledge about the polynomials by differentiating with respect to ω. Surprisingly, this differentiation yields the orthogonal polynomial of smaller degree.
where the proportionality constant, which depends on ω, is one of the recurrence coefficients (3) of the orthogonal polynomials.
Differentiating with respect to ω yields (11) Since F k (ω) = 0 for k = 0, . . . , n − 1 and F k+1 (ω) = 0 for k = 0, . . . , n − 2, we must have that Because p ω n is monic, ∂p ω n ∂ω is a polynomial of degree at most n − 1. By the relations above, it satisfies the exact same orthogonality conditions as p ω n−1 . Thus, there must be a constant depending on ω, but not on x, such that It remains to determine this constant. Setting k = n − 1 in (11) we find Note from the orthogonality of p ω n that Similarly, we find that Combined with (12) this leads to the result.

Three-term recurrence relation
As it turns out, the coefficients α k and β k of the three-term recurrence relation (2) can be computed themselves by a recursion. The following result follows from taking the derivative with respect to ω of the recurrence relation. Note that α k and β k are always functions of ω. In the following we frequently omit this dependency and we denote the derivative of α k with respect to ω simply by α k and we do so similarly for β k . Theorem 3.3 Let α k and β k be the ω-dependent recurrence coefficients for p ω n , defined by (3). Then we have Proof 3 Differentiating the recurrence (2) with respect to ω yields Using Thm. 3.2 and collecting terms leads to Matching the leading order coefficients on both sides, using that both p ω k and p ω k−1 are monic, we must have −iβ k+1 + ∂α k ∂ω = −iβ k . The recurrence for β k follows from this. Dividing over yields and the recurrence for α k follows by comparing to the regular recurrence for Remark 1 In the theory of orthogonal polynomials, these are sometimes known as the deformation equations for the recurrence coefficients. General expressions of this kind can be obtained whenever a weight function of exponential type is perturbed with a parameter, see for instance [2, Prop. 2.1] and references therein.

Trajectories of the roots
Finally, we intend to show that the trajectories of the roots in the complex plane may have cusps. If so, the norm of the corresponding orthogonal polynomial vanishes, and as a result the orthogonal polynomial of one degree higher does not exist. We start by showing the equivalence between the vanishing derivatives of the roots x j and the vanishing derivatives of the polynomial p ω n . As in the previous section we frequently omit the dependency of the roots on ω in our notation.
Theorem 3.4 Assume that p ω n exists and let x j (ω), j = 1, . . . , n, denote its zeros. If, for a given ω * , we have then (p ω * n , p ω * n ) = 0. If p ω n is uniquely defined, then the converse is also true.
Proof 4 Writing p ω n in terms of its factors, we have Differentiating with respect to ω yields From this and condition (13) it follows that ∂p ω n ∂ω ω=ω * ≡ 0, i.e., the partial derivative of p ω n vanishes identically. By Theorem 3.2 we find that ω * must be a root of β n (ω), from which we conclude in turn that (p ω n , p ω n ) must vanish at ω = ω * .
Recall the functions F k (ω) from (10) and their derivatives (11). Letting k = n − 1, we find from the above and (11) that To prove the converse, we assume (p ω * n , p ω * n ) = 0, which leads using (11) again to This means that ∂p ω n ∂ω ω=ω * satisfies the orthogonality conditions of p ω n . However, it is a polynomial of degree n − 1 and since p ω n exists uniquely, the lower-degree polynomial can only satisfy the above conditions if it vanishes identically.
If (p ω * n , p ω * n ) = 0 for some value of n and ω * , then the monic orthogonal polynomial p ω * n+1 of degree n + 1 does not exist. In fact, the polynomial p ω * n satisfies all of the orthogonality conditions that p ω * n+1 should satisfy, since This explains why the roots of p ω 3 agree with those of p ω 2 in Figure 5 at specific values of ω. The third root, as a function of ω, has poles at these values of ω.

Properties of the quadrature rule
Since the weight function considered here is non-positive, only few of the classical results for orthogonal polynomials apply. We do not fully settle the questions of existence and uniqueness of the polynomials in this paper. Yet, several interesting properties of the quadrature rule can be established and we start with the asymptotic order.

Asymptotic order
The following theorem establishes a connection between polynomial accuracy and asymptotic accuracy, on the assumption that the quadrature points cluster near the endpoints. Note that a Gaussian quadrature rule with an odd number of points in our setting does not qualify here, since due to the symmetry at least one root is always on the imaginary axis. In the notation of the following theorem, a Gaussian rule with an even number of points corresponds to N = 4n and M = 2n, and the theorem thus explains the error behaviour we observe in Fig. 4. x k e iωx dx, k = 0, 1 . . . , N − 1.
and assume that N ≥ 2n, i.e., the rule is at least interpolatory. If the nodes x j can be split in two groups {x 1 j } n j=1 and {x 2 j } n j=1 , such that x 1 j = −1 + O(ω −1 ) and x 2 j = 1 + O(ω −1 ), j = 1, . . . , n, then, for an analytic function f the error has the asymptotic decay Proof 5 From [16] we have that an analytic function f can be written as Applying the method to the j-th Lagrange polynomial l j (x), we have, To conclude that w j = O(ω −1 ), we need that l (k) j (±1) = O(ω k ). Now assume the node x j is the first member of the group 1: x j = x 1 1 . Writing l j (x) in terms of the two groups of nodes we have The second of these products is clearly bounded with its derivatives, and we can, by Leibniz's formula, concentrate on the first. Here, the denominator is of order O(ω −n+1 ). Similarly, the numerator is of order O(ω −n+1 ), when evaluated in x = −1. Differentiating will lower the degree of the numerator, Evaluating in x = −1 gives that the numerator is of order O(ω −n+2 ). In general the k-th derivative of the numerator is of order Similarly one shows that l

4.2
Large ω behaviour of p ω n Next, a rather remarkable fact will be demonstrated, namely that pointwise we have that where L n (x) is the classical Laguerre polynomial of degree n. That is, the orthogonal polynomial becomes the product of two rescaled classical orthogonal polynomials in the limit ω → ∞. This is, in fact, the polynomial that vanishes at the superinterpolation points. To prove this we need the following intermediate result: Lemma 4.2 Consider a vector of n > 1 points in C, x 1 , . . . , x n , and a sequence of integers α 1 < α 2 < . . . < α n . We construct the generalised where L(x 1 , x 2 , . . . , x n ) is a polynomial in x 1 , . . . , x n .
By expanding the determinant along any row, it is apparent that H(x 1 , . . . , x n ) is actually a polynomial. Now H(x j , x 2 , . . . , x n ) = 0, j = 2, 3, . . . , n, since this corresponds to a matrix with duplicate rows. This implies that x 1 − x j , j = 2, . . . , n is a factor of H. Similarly, if one considers H in terms of the remaining arguments, one sees that x i − x j , i = j, are all factors of H. The result follows from this.
The key to showing something like (15), is to look at p ω n evaluated in the superinterpolation points. Applying the interpolatory quadrature rule based on interpolation at the superinterpolation points, along with the assumption on the boundedness of p ω 2n , gives [5, Th.
Here the weights are given in terms of the Gauss-Laguerre weights η j , where V = {x j−1 i } 2n i,j=1 . By Cramer's rule we have Entries of Adj(V ) are computed in terms of determinants: where V j,i is V with row j and column i deleted. Using this we can find the asymptotics of V −1 elementwise. First, using the well known formula for the Vandermonde determinant, where C = 0, which follows from the fact that the points belong two groups of n points which are ω −1 close to either −1 or 1. Similarly, if we delete one row and one column of V , one of the groups will be missing a point, and V j,i will be a generalised Vandermonde matrix. The form of the determinant is given by Lemma 4.2, and from this it follows that det(V j,i ) = O(ω −(n−1) 2 ), 0 < i, j ≤ 2n.
Finally we can state the correspondence between roots of p ω 2n and superinterpolation points, ω → ∞.
From this expression, it is clear that for each j = 1, . . . , 2n we must have y j ∼ ±1, which gives O(ω −n ). The only way O(ω −n−1 ) can be attained is if for some l.
Corollary 1 Assume p ω 2n exists, and that it is bounded with all its derivatives in ω. The 2n point Gaussian rule, with nodes being the zeros of p ω 2n , is of asymptotic order n + 1. Note that in the last two results we assumed that the polynomial p ω 2n is bounded in ω, along with all its derivatives with respect to x. This, along with the assumption that the polynomials exist for all ω, is not proved in this paper.

Conclusions and outlook
We have presented a Gaussian quadrature rule for integrals on [−1, 1] with weight function e iωx , ω ≥ 0. The associated family of orthogonal polynomials p ω n is non standard because of the oscillatory nature of the weight function. It is shown that for some critical values of ω, the orthogonal polynomials of odd degree fail to exist, and that phenomenon is related to the fact that the bilinear form defined with respect to e iωx is not positive definite. However, based on numerical experiments we conjecture that all polynomials of even degree exist for any value of ω. We give some results on this family of orthogonal polynomials, but their global behaviour in terms of ω is bound to be very complicated. We note that they should bridge between the Legendre case (when ω = 0) and a product of two rotated and scaled Laguerre polynomials in the complex plane, as ω → ∞.
We also present results connecting the concept of standard polynomial accuracy of Gaussian quadrature rules and that of asymptotic order in terms of ω, which is of great interest in the context of highly oscillatory quadrature. These two ideas are related under the assumption that the zeros of p ω n cluster near the endpoints, something that is observed numerically in the paper.