## Abstract

In this paper, we explore quadrature methods for highly oscillatory integrals. Generalizing the method of stationary phase, we expand such integrals into asymptotic series in inverse powers of the frequency. The outcome is two families of methods, one based on a truncation of the asymptotic series and the other extending an approach implicit in the work of Filon (Filon 1928 *Proc. R. Soc. Edinb.* **49**, 38–47). Both kinds of methods approximate the integral as a linear combination of function values and derivatives, with coefficients that may depend on frequency. We determine asymptotic properties of these methods, proving, perhaps counterintuitively, that their performance drastically improves as frequency grows. The paper is accompanied by numerical results that demonstrate the potential of this set of ideas.

## 1. Introduction

The quadrature of highly oscillating integrals is a computational problem of overarching importance in a wide range of applications, e.g. quantum chemistry, image analysis, electrodynamics, computerized tomography and fluid mechanics. It is widely perceived as a ‘difficult’ problem, best overcome by somehow eliminating oscillation, e.g. by choosing exceedingly small subintervals. The contention of this paper is that, once appropriate discretization methods are designed and properly understood, the problem of highly oscillatory quadrature becomes relatively simple and that the precision of the calculation actually *increases* as the frequency of oscillation grows. This phenomenon has already been identified in Iserles (2003*a*,*b*; see also Levin 1996 for an earlier and different computational approach). Here, we extend the approach of Iserles (2003*a*,*b*) and show that the rate of decay of the error, once frequency grows, can be increased arbitrarily by the inclusion of higher derivatives.

Amazingly enough, the origin of the thread that will lead us to a new generation of highly effective methods is a paper that appeared more than seventy years ago. The standard numerical approach toward the discretization of the integral(1.1)where *f* and *g* are suitably smooth functions, is *Gauss–Christoffel quadrature*, where we interpolate the integrand at distinct *nodes c*_{1}<*c*_{2}<⋯<*c _{ν}* in [0, 1] by a polynomial

*p*of degree

*ν*−1, and approximate(Davis & Rabinowitz 1980). Unfortunately, if

*ω*≫1, this results in a completely useless estimation of

*I*[

*f*], bearing an error of O(1) (Iserles 2003

*a*). An alternative approach builds upon an idea first advanced in a special case by Louis Napoleon George Filon (1928) and phrased in more modern terminology by Flinn (1960). Thus, instead of interpolating the integrand

*f*(

*x*)e

^{iωg(x)}, we interpolate at

*c*

_{1},…,

*c*the function

_{ν}*f*(

*x*) by a polynomial and let(1.2)where and ℓ

_{k}is the

*k*th cardinal polynomial of Lagrangian interpolation. Note that the construction of requires the availability of the first few moments in an explicit form. This is tacitly assumed in the remainder of the paper. It was proved in Iserles (2003

*a*) that as long as

*g*′≠0 in [0, 1], it is true that(1.3)

Thus, increased frequency *ω* results in a smaller error and this behaviour is enhanced once the endpoints are included among quadrature nodes. If *g*′ vanishes at *ξ*_{1},…,*ξ*_{s}∈(0, 1), say, and nowhere else in [0, 1], and if *g*″(*ξ*_{k})≠0, *k*=1,…,*r*, then once both the endpoints and the stationary points *ξ*_{1},…,*ξ*_{r} are nodes, we obtain The theory can be extended to more general stationary points, when (Iserles 2003*b*).

The current paper pursues the idea of letting the interpolating polynomial depends on derivatives of *f*, thereby improving upon the asymptotic estimate (1.3). In effect, in place of Lagrange interpolation in the Filon method (1.3), we allow Hermite interpolation. This requires substantive extension of the asymptotic analysis from Iserles (2003*a*,*b*).

In §2, we consider the case when *g*′≠0, while §3 is devoted to the case when *g* is allowed stationary points in [0, 1]. Our main basic tool is the method of *stationary phase* (Stein 1993), except that to derive high-order asymptotic estimates in a convenient form, we need to resort to a new, albeit relatively straightforward, method of analysis.

The work of §§2 and 3, and the numerical results therein, lead to two new families of efficient methods for the quadrature of highly oscillatory integrals: a generalization of Filon's method (1.2), as well as a straightforward truncation of an asymptotic expansion, which we term as an *asymptotic method.* Early numerical results indicate that although both approaches have merit, a generalized Filon method is usually better. More detailed analysis of the two approaches and means for a practical estimation of numerical error are a matter for current research and will feature in a forthcoming publication.

Our results can be easily and transparently extended from [0, 1] to any bounded interval by linear mapping, while the situation for unbounded intervals is, if at all, substantially easier. As will be clear from §§2 and 3, endpoints contribute to the asymptotic expansion and, once an interval is infinite, the integral (1.1) exists and the moments *m*≥0, are bounded, this contribution disappears.

## 2. The case when *g*′(*x*)≠0, 0≤*x*≤1

We commence from the integral (1.1), assuming that the function *g* is strictly monotone in the interval [0, 1]. Both here and in the sequel, we assume that both *f* and *g* are smooth. Generalization to *C*^{r}[0, 1] functions for sufficiently large *r*≥1 is mathematically straightforward but notationally clumsy, which obscures the concepts at the basis of our work. Although the leading term in asymptotic series of the integral as *ω*→∞ can be derived from the van der Corput lemma (Stein 1993), we require more terms in an explicit form. They can be obtained from the following lemma.

*Let*

*Then, for ω*→∞,(2.1)

We prove by induction on *s*≥0 the identity(2.2)

It is certainly true for *s*=0, coinciding with equation (1.1). For *s*≥1, we integrate on the right by parts,

This proves equation (2.2) and by letting *s*→∞, the proof of the lemma is complete.▪

In the simplest, yet most important, case *g*(*x*)=*x*, the expansion (2.1) has a particularly simple form(2.3)which has already been proved in Håvie (1973) and was recently used for the computation of highly oscillatory integrals by Degani & Schiff (2003).

For general monotone *g* we have, by direct calculation,and it is trivial to prove that for every *k*≥0 and *j*=0,1,…,*k*, there exist functions *σ*_{k,j}, which depend upon *g* and its derivatives, such that(2.4)

An immediate consequence of lemma 2.1 is that the *asymptotic method*,(2.5)a truncation of the asymptotic expansion (2.1), represents an efficient approximation to *I*[*f*] once *ω* is large.

*For every smooth f and g it is true that*(2.6)

This follows at once from equation (2.2) since, replacing by *σ _{x}*[] in equation (2.1), we observe that ▪

As an alternative to , we choose an integer *s*≥1, nodes *c*_{1}=0<*c*_{2}<⋯<*c*_{ν−1}<*c*_{v}=1, and integers *θ*_{1}, *θ*_{2},…,*θ _{ν}*≥1. We let be a Hermite polynomial approximation to

*f*at the points

*c*

_{1},…,

*c*of multiplicities

_{ν}*θ*

_{1},…,

*θ*, respectively. In other words, is a polynomial of degree

_{ν}*n*−1, where such that

*j*=0, 1,…,

*θ*−1,

_{l}*l*=1,2,…,

*ν*. It is well known thatwhere each

*α*

_{l,j}is itself a polynomial of degree

*n*−1 such that for all

*i*=0, 1,…,

*θ*

_{l}−1 and

*m*=1, 2,…,

*ν*, except that (Lorenz

*et al.*1983). Consistent with equation (1.2), we define the

*generalized Filon method*as(2.7)where

*θ*

_{1},

*θ*≥

_{ν}*s*and

*Suppose that θ*_{1}, *θ _{ν}*≥

*s. For every smooth f and g, it is true that*(2.8)

We replace by in the asymptotic formula (2.1). Since *j*=0, 1,…,*s*−1, it follows from equation (2.4) that *k*=0, 1,…,*s*−1. This implies the asymptotic formula (2.8). ▪

As an example, we take *g*(*x*)=*x*, *s*=2 and, in Filon's method, *ν*=2, *c*_{1}=0, *c*_{2}=1 and *θ*_{1}=*θ*_{2}=2. The asymptotic and generalized Filon methods are, respectively,

At first glance, the asymptotic method has an obvious advantage in the simplicity of its coefficients. On the other hand, while an asymptotic method blows up for *small* frequency, this is not the case with generalized Filon. In fact, it was proved in Iserles (2003*a*) that for *g*(*x*)=*x* and 0<*ω*≪1, it is true that for all polynomials *f* of degree *p*−1, where *p* is the order of the usual Gauss–Christoffel quadrature at the nodes *c*_{1}, *c*_{2},…,*c _{ν}*. Although this result for general function

*g*is no longer valid (Iserles 2003

*b*), it is nevertheless true from our use of interpolatory quadrature weights that the degree of polynomial recovery in general is at least

*n*−1.

However, the most striking advantage of generalized Filon emerges in numerical experiments. Thus, in figure 1 we display the scaled errors and for *f*(*x*)=cos *x* and *g*(*x*)=*x*. While both methods display the expected error decay for large , performs significantly better. This is confirmed by a raft of other computer experiments. Moreover, although the rate of decay in the error of generalized Filon cannot be improved by taking *ν*≥3, thereby adding internal nodes, this procedure seems to decrease the leading error constant drastically. This can be clearly observed when comparing figures 1 and 2.

In figure 3, we consider the error for *f*(*x*)=e^{x}, *g*(*x*)=(1+*x*)^{2} and compare errors in four methods — the asymptotic method and three generalized Filon methods — with (i) ** c**=[0, 1],

*=[2, 2], (ii)*

**θ***=[2, 1, 2] and (iii)*

**θ***=[2, 2, 2]. To display all the four errors with disparate orders of magnitude in a single figure, we have computed the quantities log |*

**θ***Q*[

*f*]−

*I*[

*f*]|+3 log

*ω*, where

*Q*[

*f*] stands for the underlying quadrature formulae. The rapid improvement in generalized Filon for increasing

*n*is clear.

## 3. Generalized Filon integration in the presence of stationary points

Once *g*′ vanishes at one or more points in [0, 1], the asymptotics of *I*[*f*] are given by the classical method of stationary phase. This statement is, however, somewhat misleading. In principle, the method of stationary phase applies in two situations.

Instead of integrating in [0, 1], we integrate in (−∞,∞) with the proviso that both

*f*and*g*are in*L*[]∩*C*^{∞}[] (Olver 1974); orWe integrate in [0, 1] but require that

*f*has compact support sufficiently near the stationary points of*g*(Stein 1993).

However, the situation can be readily remedied by partitioning *f* into a sum of bump functions and a remainder, as explained in Iserles (2003*b*). The outcome is an important *superposition principle*, namely that the asymptotic form of equation (1.1) is determined by the local behaviour at each stationary point, which can be investigated by the method of stationary phase, *as well as* the contribution of the endpoints, as given by lemma 2.1.

In fact, we require more than the standard method of stationary phase can deliver, namely the explicit form of all the expansion coefficients. This will be provided in the next lemma, which, in a sense, does to the method of stationary phase what lemma 2.1 did to the van der Corput lemma.

We denote the *generalized moments* of the functional *I* byand note that for every *ξ*∈[0, 1] and every pair of smooth functions *f* and *g*, such that *g*′(*ξ*)=0, *g*″(*ξ*)≠0 and *g*′(*x*)≠0 in *x*∈(0, 1)\{*ξ*}, it is true that(3.1)

*Suppose that g is a smooth function and that g*′(*ξ*)=0, *g*″(*ξ*)≠0, *for some ξ*∈(0, 1) *and g*′(*x*)≠0, *x*∈(0, 1)\{*ξ*}. Then for every smooth *f*, it is true that(3.2)where

Similar to the proof of lemma 2.1, we commence by proving by induction that(3.3)

This is certainly true for *s*=0. Moreover, by equation (3.1),and integration by parts, followed by substitution in equation (3.3), completes the inductive proof. Letting *s*→∞ yields the asymptotic estimate (3.2).▪

The asymptotic expansion (3.2) can be readily extended to the case when *g*^{(r+1)}(*ξ*)≠0 for any integer *r*≥1. In place of equation (3.1), we thus have

Consequently, lettingwe have, following a proof identical to that of lemma 3.1, the asymptotic estimate(3.4)for every smooth *f*.

Once we can derive an asymptotic expansion of the integral (1.1) for a single stationary point in (0, 1), we can do so for any finite number of stationary points. Thus, if *ξ*_{1}, *ξ*_{2},…,*ξ*_{q}∈(0, 1) are all the stationary points (perhaps of different orders) of *g*, while *g*′(0), *g*′(1)≠0, we can partition the interval [0, 1] into *q* subintervals _{1}, _{2},…,_{q} such that a single stationary point resides inside each _{j}. Therefore, equation (3.4) applies (with trivial modification to cater for different endpoints) in every _{j}. A general asymptotic formula can easily be written down but such generality probably obscures, rather than illuminates, the issue.

We note, in passing, that the value of *I*[*f*] is trivially *independent* of the partition. In other words, the contribution of endpoints, except for 0 and 1, necessarily cancels, at least in principle, an asymptotic expression for *I*[*f*] samples *f* and its derivatives at the endpoints 0 and 1 and the stationary points *ξ*_{1}, *ξ*_{2},…,*ξ*_{p}. However, this is not as helpful as it sounds since the choice of endpoints of _{j} expresses itself in the asymptotic formula *inter alia* in the definition of generalized moments. Once we attempt to discard the contribution of ‘internal’ endpoints, we are rapidly led to fairly complicated and opaque expressions. On the other hand, we pre-empt our discussion of generalized Filon methods in this setting by noting that, as a consequence of our discussion, they can be formed by interpolating to requisite order just at 0,*ξ*_{1}, *ξ*_{2},…,*ξ*_{p}, 1.

Of course, to appreciate the rate of decay of *I*[*f*] as *ω*≫1 grows, we must determine the behaviour of the generalized moments *μ*_{j}(*ω*;*ξ*) for *j*=0, 1,…,*r*−1. This information follows from the method of the stationary phase, at least to the extent that if *ξ*∈(0, 1), *g*′(*ξ*)=⋯=*g*(*r*)(*ξ*)=0, *g*^{(r+1)}(*ξ*)≠0 and *g*′(*x*)≠0 elsewhere in [0, 1], then there exist linear differential functionals *a*_{k}, *k*=0, 1,…, such that for every smooth function *h*, compactly supported near *ξ*,(3.5)where *a*_{m}[*h*] are functions of *ω*, bounded for *ω*≫1 (Stein 1993). The term *a*_{0}[*h*] is explicitly known and is a multiple of *h*(*ξ*). In addition, if *h*(*x*)=O((*x*−*ξ*)^{j}) for some *ξ*∈(0, 1) then it is easy to see that *a*_{k}[*h*]=0 for *k*≤*j*. Therefore,

Our first conclusion is that the asymptotic method (2.5) can be extended to cater for stationary points. The following result follows readily from our discussion.

*Assume that g*′(*ξ*)=*g*″(*ξ*)=⋯=*g*^{(r)}(*ξ*)=0, *g*^{(r+1)}(*ξ*)≠0 and *g*′≠0 *in* [0, 1]\{*ξ*}. *Then, for every s*≥0 *the method*(3.6)*carries for every smooth function f the asymptotic error*

Likewise, we can extend the generalized Filon method (2.7). We need first, however, to examine in greater detail the dependence of the operators *ρ*_{k}[*f*] on the derivatives of *f*. It is easy to see that for *x*≠*ξ* each *ρk*[*f*](*x*) is a linear combination of *f*(*x*), *f*′(*x*),…,*f*^{(k)}(*x*) and this follows similarly to our treatment of operators *σ*_{k} in §2. However, things are crucially different at *x*=*ξ*, a value that, according to equation (3.3), we must sample. For example, for *r*=1 we have by direct computationand so on. It is easy to prove that each *ρ*_{k}[*f*](*ξ*) for *k*≥1 depends on *f*^{(j)}(*ξ*), *j*=1, 2,…,2*k*.

The situation is more messy for general *r*≥1. Let

Then, after straightforward algebra,

Therefore, *ρ*_{1}[*f*](*ξ*) depends linearly on *f*^{(r)}(*ξ*) and *f*^{(r+1)}(*ξ*). More generally, lettingwe observe that each is a linear combination of *f*_{r}, *f*_{r+1},…,*f*_{r+m+1}. Therefore, applying the same argument to *ρ*_{1}[*f*] as to *f*, we deduce that *ρ*_{2}[*f*](*ξ*) is a linear combination of *f*^{(j)}(*ξ*), *j*=*r*, *r*+1,…,2*r*+2. In general, an inductive argument shows that *ρ*_{k}[*f*](*ξ*) is a linear combination of *f*^{(j)}(*ξ*) for *j*=*r*, *r*+1,…,*k*(*r*+1) for every *k*≥1.

We are now in a position to define and analyse the generalized Filon method in the present setting. The important observation, following Iserles (2003*b*) in the simplest case *s*=1, is that stationary points play a similar role to endpoints and, to improve asymptotic decay of the error, we must sample function values and derivatives there. Specifically, we must calculate sufficient derivatives to be able to compute sufficient values of *ρ*_{k}[*f*](*ξ*).

We choose *ν*≥3 and nodes *c*_{1}<*c*_{2}<⋯<*c _{ν}* such that

*c*

_{1}=0,

*c*=1, and there exists

_{ν}*q*∈{2, 3,…,

*ν*−1} such that

*c*

_{q}=

*ξ*. Given integers

*θ*

_{1},

*θ*

_{2},…,

*θ*≥1, we let be a polynomial of degree , such that(3.7)

_{ν}*Let θ*_{1}, *θ*_{ν},≥*s and θ*_{q}≥*s*(*r*+1). *Then*(3.8)

The proof is virtually identical to that of theorem 2.3. We substitute in place of *f* in equation (3.5), and note that Hermite interpolation conditions at the endpoints and the stationary points, together with the asymptotic formula (3.5), render zero all the leading terms of the asymptotic expansion, whereby confirming (3.8).▪

The practical application of either or requires knowledge of the first few generalized moments of *g* in an explicit form. In the important case *g*(*x*)=(*x*−*ξ*)^{p}, where *ξ*∈(0, 1) and *p*≥1 is an integer, the moments are easily calculated:where, changing a variable,where *Γ*(*z*, *α*) is the *incomplete Gamma function* (Abramowitz & Stegun 1964). In general, however, moments are often unknown. In addition to the need to formulate the underlying problem in the form (1.1), i.e. as a product of a slowly varying function and a fast oscillator, this places obvious restrictions on the applicability of asymptotic and generalized Filon methods alike.

As an example of asymptotic and generalized Filon methods, consider the integral

Hence, and *r*=1. The asymptotic formula (3.4) (actually, (3.2) suffices in this setting) yields, after some algebra,where erf *z* is the familiar error function.

For , we take the least values consistent with a O(*ω*^{−5/2}) asymptotic decay: *ν*=3, *θ*_{1}=*θ*_{3}=2 and *θ*_{2}=3. In other words, we interpolate to the function values *f*(0), *f*′(0), *f*(1) and *f*′(1) with a sixth-degree polynomial. It is illustrative to present the coefficients explicitly. Thus,where

In figure 4 we display the error, scaled by *ω*^{5/2}, for both methods and *f*(*x*)=e^{x}. Note that, although the asymptotic decay is of the order O(*ω*^{−5/2}), as predicted by theorems 3.2 and 3.3, the error constant is smaller by four degrees of magnitude in the case of generalized Filon method.

Although one can present many other numerical results, this course of action is likely to result in little further enlightenment. Instead, we complete our discussion of stationary points by considering the case when *g*′ vanishes at an endpoint. Without loss of generality, we assume that *g*′(0)=0 and *g*′≠0 in (0, 1], bearing in mind that the case when *g*′(0)=*g*′(1)=0 can be resolved by subdividing [0, 1].

We employ a technique similar to that in the proof of lemma 3.1. For simplicity, we assume that *g*″(0)≠0 but the scope of our analysis can be easily extended. Because letting *ξ*=0, *I*[*f*] satisfies equation (3.1), integration by parts and a single application of the l'Hôpital rule demonstrate that

We deduce at once, revisiting the inductive proof of lemma 3.1, that as *ω*→∞(3.9)

Here, *ρ*_{k} are the operators that we have defined earlier in the section, except that *ξ*=0. Let

Then,and, in general, each *ϕ*_{k} is a linear combination of *f*^{(j)}(0), *j*=1, 2,…,2*k*+1. Likewise, each *ρk*[*f*](0) and *ρ*_{k}[*f*](1) is a linear combination of *f*^{(j)}(0), *j*=0, 1,…,2*k*, and *f*^{(j)}(1), *j*=0, 1,…,*k*.

Commencing from equation (3.9), we conclude at once that the asymptotic methodcarries an asymptotic error of O(ω^{−s−(1/2)}). Moreover, by our analysis, it combines the values of *f*^{(j)}(0), *j*=0, 1,…,2*s*−1, as well as *f*^{(j)}(1), *j*=0, 1,…,*s*−1.

Likewise, the generalized Filon method can be extended to this setting, as long as *c*_{1}=0, *c*_{ν}=1, *θ*_{1}≥2*s* and *θ*_{ν}≥*s*.

## Acknowledgements

The authors wish to thank Brad Baxter (Birkbeck College, London), Hermann Brunner (Memorial University of Newfoundland), Fernando Casas (Universitat Jaume I, Castellón) and Alexei Shadrin (DAMTP, University of Cambridge) for their advice. We also wish to thank the anonymous referees for their helpful comments. The work of the second author was performed while a Visiting Fellow of Clare Hall, Cambridge, during sabbatical leave from Norwegian University of Science and Technology.

## Footnotes

- Received January 27, 2004.
- Accepted September 15, 2004.

- © 2005 The Royal Society