Highly Oscillatory Quadrature: The Story soFar

TL;DR: Iserles et al. as discussed by the authors presented an analysis of the relationship between the two types of models and showed that Olver's model is more similar to the one proposed by Norsett and Serles.

Abstract: 1 Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom, A.Iserles@damtp.cam.ac.uk 2 Department of Mathematical Sciences, Norwegian University of Science and Technology, 7491 Trondheim, Norway, S.P.Norsett@math.ntnu.no 3 Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom, S.Olver@damtp.cam.ac.uk

Summary

Introduction

• The last few years have witnessed substantive developments in the computation of highly oscillatory integrals in one or more dimensions.
• All these methods share the surprising property that their accuracy increases with growing oscillation.

1 The challenge of high oscillation

• Rapid oscillation is ubiquitous in applications and is, by common consent, considered a ‘difficult’ problem.
• If the integrand oscillates rapidly, and unless the authors use an astronomical number of function evaluations, polynomial interpolation is useless!.
• The purpose of the final section is the sketch gaps in the theory and comment on ongoing challenges and developments.
• Moreover, the authors describe there briefly the recent method of Huybrechs & Vandewalle (2005), as well as the work in progress in Cambridge and Trondheim.
• It is thus of interest to mention that the availability of efficient highly oscillatory quadrature is critical to a number of contemporary methods for ordinary differential equations that exhibit rapid oscillation (Degani & Schiff 2003, Iserles 2002, Iserles 2004, Lorenz, Jahnke & Lubich 2005).

2 Asymptotic expansion in the absence of critical points

• The authors restrict their analysis to R2, directing the reader to (Iserles & Nørsett 2006) for the general case.
• Because of the nonresonance condition, the gradient of the oscillator does not vanish in any of these simplices and the authors can apply (3) therein: this expresses I[f,S] as an asymptotic expansion over (n − 2)-dimensional simplices.
• The authors continue with this procedure until they reach 0-dimensional simplices: the n+1 vertices of the original simplex.
• A moment’s reflection clarifies that only the original vertices of Ω may influence the expansion: the internal vertices are arbitrary, since there is an infinity of simplicial complexes consistent with the nonresonance condition.

3.1 Asymptotic methods

• The simplest and most natural means of approximating (1) consists of a truncation of the asymptotic expansion (5) (replacing S by a polytope Ω).
• Asymptotic quadrature is particularly straightforward in a single dimension, since then its coefficients are readily provided explicitly by an affine mapping of (4) from (0, 1) to an arbitrary bounded real interval.
• Note that all the coefficients are well defined, because of the nonresonance condition.
• Another important shortcoming of an asymptotic method is that, given ω and the number of derivatives that the authors may use, its accuracy, although high, is predetermined.

3.2 Filon-type methods

• A more sophisticated use of the asymptotic expansion rapidly leads to far superior, accurate and versatile quadrature schemes.
• The authors will return to this restriction upon the applicability of (8) in the sequel.
• It is important to observe that in the ‘minimalist’ case, when ϕ interpolates only at the vertices of Ω, (7) and (8) use exactly the same information.
• The difference in their performance, which is often substantive, is due solely to the different way this information is processed.
• In one dimension the authors construct Filon-type methods similarly to the familiar interpolatory quadrature rules.

4 Critical points

• Worse, in a multivariate setting surprisingly little is known about asymptotic expansions in the presence of high oscillation and critical points (Stein 1993).
• The situation is much clearer and better understood in a single dimension.
• Note that (12) is not a ‘proper’ asymptotic expansion, because of the presence of the function µ0(ω).
• Assuming that µ0 can be computed – and the authors need this anyway for Filon-type methods!.
• (13) We can easily cater for any number of critical points, possibly of different degrees, once the authors include them among the nodes and choose sufficiently large multiplicities.the authors.

5 Conclusions and pointers for further research

• The first and foremost lesson to be drawn from their analysis is that, once the authors can understand the mathematics of high oscillation, they gain access to a wide variety of effective and affordable algorithms.
• Yet, once the authors concern ourselves with bounded domains with boundary and allow for the presence of critical points, a great deal remains to be done.
• Moreover, even if the explicit form of (5) is unavailable, the very existence and known structure of an asymptotic formula allow us to analyse better and more flexible quadrature methods.
• The authors have already touched upon applications of highly oscillatory quadrature to numerical methods for rapidly oscillating differential equations.

Figures

Fig. 7. Scaled error for QL1 (on the left) and Q L 2 (on the right) for Ω = H, f(x, y) = ex−2y and g(x, y) = x3 + x− y.

Fig. 1. Error in Gaussian quadrature with Ω = (−1, 1), f(x) = cosx, g(x) and ν points. Here ν increases by increments of five, from 5 to 30.

Fig. 2. Error, scaled by ωp, in asymptotic quadrature of asymptotic order p with Ω = (−1, 1), f(x) = x sinx, g(x) = x+ 1 4 x2 and s = 1, p = 2 (on the left) and s = 2, p = 3 (on the right).

Fig. 4. Scaled error for QF1 with c = [−1, 1], m = [1, 1] (on the left), Q F 1 with c = [−1,− 3 4 , 3 4 , 1], m = [1, 1, 1, 1] (at the centre) and QA2 with c = [−1, 1], m = [2, 2] (on the right) for Ω = (−1, 1), f(x) = x sinx and g(x) = x+ 1 4 x2.

Fig. 5. Scaled error for QF1 (on the left) and Q F 2 (on the right) for Ω = S2, f(x, y) = ex−2y and g(x, y) = x+ 2y.

Fig. 6. Scaled error for QL1 with c = [−1, 1], m = [1, 1] (on the left), Q L 1 with c = [−1,− 3 4 , 3 4 , 1], m = [1, 1, 1, 1] (at the centre) and QL2 with c = [−1, 1], m = [2, 2] (on the right) for Ω = (−1, 1), f(x) = x sinx and g(x) = x+ 1 4 x2.

Fig. 8. The error for QA1 (on the left) and Q A 2 (on the right), scaled by ω 3 2 and ω 5 2 respectively, for Ω = (−1, 1), f(x) = cosx and g(x) = x2, with a stationary point at the origin.

Fig. 9. Scaled error for QF1 with c = [−1, 0, 1], m = [1, 3, 1] (on the left), Q F 1 with c = [−1,− 1 2 , 0, 1 2 , 1], m = [1, 1, 3, 1, 1] (at the centre) and QA2 with c = [−1, 0, 1], m = [2, 5, 2] (on the right) for Ω = (−1, 1), f(x) = cosx and g(x) = x2.

Fig. 3. Scaled error for QA1 (on the left) and Q A 2 (on the right) for Ω = S2, f(x, y) = ex−2y and g(x, y) = x+ 2y.

Full text

Highly oscillatory quadrature: The story so far
A. Iserles
1
, S.P. Nørsett
2
, and S. Olver
3
1
Department of Applied Mathematics and Theoretical Physics, Centre for
Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge
CB3 0WA, United Kingdom, A.Iserles@damtp.cam.ac.uk
2
Department of Mathematical Sciences, Norwegian University of Science and
Technology, 7491 Trondheim, Norway, S.P.Norsett@math.ntnu.no
3
Department of Applied Mathematics and Theoretical Physics, Centre for
Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge
CB3 0WA, United Kingdom, S.Olver@damtp.cam.ac.uk
Summary. The last few years have witnessed substantive developments in the com-
putation of highly oscillatory integrals in one or more dimensions. The availability
of new asymptotic expansions and a Stokes-type theorem allow for a compre hensive
analysis of a number of old (although enhanced) and new quadrature techniques:
the asymptotic, Filon-type and Levin-type methods. All these methods share the
surprising property that their accuracy increases with growing oscillation. These
developments are described in a unified fashion, taking the multivariate integral
R
Ω
f(x)e
iωg(x)
dV as our point of departure.
1 The challenge of high oscillation
Rapid oscillation is ubiquitous in applications and is, by common consent,
considered a ‘difficult’ problem. Indeed, the standard technique of dealing
with high oscillation is to make it disappear by sampling the signal sufficiently
frequently, and this typically leads to prohibitive cost.
The subject of this article is a review of recent work on the computation
of integrals of the form
I[f, Ω] =
Z
Ω
f(x)e
iωg(x)
dV, (1)
where Ω ⊂ R
n
is a bounded open domain with piecewise-smooth boundary,
while f and the oscillator g are smooth. We assume in (1) that ω ∈ R is large
in modulus, hence I[f, Ω] oscillates rapidly as a function of ω.
A natural technique to compute (1) in a univariate setting is Gaus sian
quadrature. Yet, a moment’s reflection clarifies that it is likely to be absolutely
useless unless |ω| is small. Classical quadrature (with a trivial weight function)
is just an exact integration of a polynomial interpolation of the integrand.

2 A. Iserles, S.P. Nørsett, and S. Olver
However, if the integrand oscillates rapidly, and unless we us e an astronomical
number of function evaluations, polynomial interpolation is useless! This is
vividly demonstrated in Fig. 1. We have computed
Z
1
−1
cos xe
iωx
2
dx = −
π
1
2
2(−iω)
1
2
exp

1
4
1
iω

erf

i
ω +
1
2
(−iω)
1
2

+ erf

i
ω −
1
2
(−iω)
1
2

by Gaussian quadrature with different number of points. The figure displays
the absolute value of the error as a function of ω ∈ [0, 100]. Note that, as
long as ω is small, everything is fine, but as so on as ω is large in comparison
with the number of quadrature points and high oscillation sets in, the error
becomes O(1). As a matter of fact, given that I[f] ∼ O

ω
−
1
2

, the trivial
approximation I[f] ≈ 0 is far superior to Gaussian quadrature with 30 points!
Yet, efficient and cheap quadrature of (1) is perfectly possible. Indeed, once
we understand the mathematical mechanism underlying (1), we can compute
it to high precision with minimal effort and, perhaps paradoxically, the quality
of approximation increases with ω.
6040
ω
80
1.4
200
1
0.8
0.6
0
100
1.2
0.4
0.2
60
0.6
0.4
0
40
ω
80200
1.2
0.8
100
0.2
1
6040
ω
200
1
0.4
0.8
0.2
0
80
0.6
100
600 40
0
20
1
0.8
0.6
10080
0.2
ω
0.4
600 4020
0.7
0.5
0.4
0.3
0
ω
80
0.6
0.2
0.1
100
6040200
0.7
0.6
0.5
0.2
0.1
ω
80
0.4
0.3
0
100
Fig. 1. Error in Gaussian quadrature with Ω = (−1, 1), f(x) = cos x, g(x) and ν
p oints. Here ν increases by increments of five, from 5 to 30.
This article collates a sequence of papers by the authors into unified narra-
tive. In particular, we revisit here the work of (Iserles & Nørsett 2005a, Iserles
& Nørsett 2006, Olver 2005a) and (Olver 2005b), to which the reader is re-

Highly oscillatory quadrature: The story so far 3
ferred for technical details, more comprehensive exposition and a wealth of
further numerical examples.
The conventional organising principle of quadrature is a Taylor expansion.
Once the integrand oscillates rapidly, a Taylor expansion converges very slowly
indeed and is, to all intents and purposes, useless. Instead, we need to exploit
an asymptotic expansion in negative powers of ω. In Section 2 we present an
asymptotic expansion of (1) in the case when the oscillator g has no critical
points: ∇g(x) 6= 0 for all x ∈ cl Ω and subject to the nonresonance condition:
∇g(x) is not allowed to be normal to the boundary ∂Ω for any x ∈ ∂Ω.
The availability of an asymptotic expansion allows us to design and analyse
effective quadrature methods, and this is the subject of Section 3. We single
out for consideration three general techniques: asymptotic methods, consisting
of a truncation of the asymptotic expansion of Section 2, Filon-type methods,
which interpolate just f(x), rather than the entire integral (Filon 1928), and
Levin-type methods, which collocate the integrand (Levin 1982).
In Section 4 we consider the case when critical points are allowed. A com-
prehensive theory exists, as things stand, only in one dimension, hence we
focus on g : [a, b] → R and study the case of g
′
(ξ) = 0 for some ξ ∈ [a, b],
g
′
6= 0 for [a, b] \ {ξ}. (Obviously, we are allowed, without loss of generality,
to assume the existence of just one critical point: otherwise we integrate in a
finite number of subintervals.) An asymptotic expansion in the presence of a
critical point presents us with new challenges. In principle, we could have used
here the standard technique of stationary phase (Olver 1974, Stein 1993), ex-
cept that it is not equal to our task. We present an alternative that leads to an
explicit and workable expansion. It is subsequently used to design asymptotic
and Filon-type methods: unfortunately, Levin-type methods are not available
in this setting.
The purpose of the final section is the sketch gaps in the theory and com-
ment on ongoing challenges and developments. Moreover, we describe there
briefly the recent method of Huybrechs & Vandewalle (2005), as well as the
work in progress in Cambridge and Trondheim.
Quadrature of (1) represents but one problem in the wide range of issues
originating in high oscillation. Quite clearly, a more significant challenge is to
solve highly oscillatory differential equations. It is thus of interest to mention
that the availability of efficient highly oscillatory quadrature is critical to a
number of contemporary methods for ordinary differential equations that ex-
hibit rapid oscillation (Degani & Schiff 2003, Iserles 2002, Iserles 2004, Lorenz,
Jahnke & Lubich 2005).
2 Asymptotic expansion in the absence of critical points
We restrict our analysis to R
2
, directing the reader to (Iserles & Nørsett 2006)
for the general case. Let first Ω = S
2
, the triangle with vertices at (0, 0), (1, 0)
and (0, 1). The nonresonance condition is thus

4 A. Iserles, S.P. Nørsett, and S. Olver
g
y
(x, 0) 6= 0, x ∈ [0, 1], g
x
(0, y) 6= 0, y ∈ [0, 1],
g
x
(x, y) − g
y
(x, y) 6= 0, x, y ≥ 0, x + y ∈ [0, 1].
Integrating by parts in the inner integral,
I[g
2
x
f, S
2
] =
Z
1
0
Z
1−y
0
g
2
x
(x, y)f(x, y)e
iωg(x,y)
dx dy
=
1
iω
Z
1
0
g
x
(1 − y, y)f(1 −y, y)e
iωg(1−y,y)
dy
−
1
iω
Z
1
0
g
x
(0, y)f(0, y)e
iωg(0,y)
dy −
1
iω
I

∂
∂x
(g
x
f), S
2

=
1
iω
Z
1
0
g
x
(x, 1 − x)f(x, 1 − x)e
iωg(x,1−x)
dx
−
1
iω
Z
1
0
g
x
(0, y)f(0, y)e
iωg(0,y)
dy −
1
iω
I

∂
∂x
(g
x
f), S
2

.
By the same token,
I[g
2
y
f, S
2
] =
1
iω
Z
1
0
g
y
(x, 1 − x)f(x, 1 − x)e
iωg(x,1−x)
dx
−
1
iω
Z
1
0
g
y
(x, 0)f(x, 0)e
iωg(x,0)
dy −
1
iω
I

∂
∂y
(g
y
f), S
2

.
Adding up, we have
I[k∇gk
2
f, S
2
] =
1
iω
(M
1
+ M
2
+ M
3
) −
1
iω
I[∇
⊤
(f∇g), S
2
],
where
M
1
=
Z
1
0
f(x, 0)[n
⊤
1
∇g(x, 0)]e
iωg(x,0)
dx,
M
2
=
√
2
Z
1
0
f(x, 1 − x)[n
⊤
2
∇g(x, 1 − x)]e
iωg(x,1−x)
dx,
M
3
=
Z
1
0
f(0, y)[n
⊤
3
g(0, y)]e
iωg(0,y)
dy.
Here
n
1
=

0
−1

, n
2
=
"
√
2
2
√
2
2
#
, n
3
=

−1
0

are outward unit normals at the edges of S
2
.
Since ∇g(x, y) 6= 0 in cl S
2
, we may replace ab ove f by f/k∇gk
2
without
any danger of dividing by zero. The outcome is

Highly oscillatory quadrature: The story so far 5
I[f, S
2
] =
1
iω
Z
∂S
2
n
⊤
(x)∇g(x)
f(x)
k∇g(x)k
2
e
iωg(x)
dS (2)
−
1
iω
I

∇
⊤

f
k∇gk
2
∇g

, S
2

.
Extending this technique to R
n
, it is possible to prove that (2) remains
true once we replace S
2
by S
n
⊂ R
n
, the regular simplex with vertices at 0
and e
1
, . . . , e
n
.
Let
f
0
(x) = f(x), f
m
= ∇
⊤

f
m−1
(x)
k∇g(x)k
2
∇g(x)

, m ∈ N.
We deduce from (2) (extended to R
n
) that
I[f
m
, S
n
] =
1
iω
Z
∂S
n
n
⊤
∇g(x)
f
m
(x)
k∇g(x)k
2
e
iωg(x)
dS −
1
iω
I[f
m+1
, S
n
], m ∈ Z
+
.
Finally, we iterate the above expression to obtain a Stokes-type formula, ex-
pressing I[f, S
n
] as an asymptotic expansion on the boundary of the simplex,
I[f, S
n
] ∼ −
∞
X
m=0
1
(−iω)
m+1
Z
∂S
n
n
⊤
∇g(x)
f
m
(x)
k∇g(x)k
2
e
iωg(x)
dS. (3)
We wish to highlight four important issues. Firstly, a trivial inductive
proof confirms that each f
m
can be expressed as a linear combination of f
and the first m directional derivatives (altogether,

n+m+1
m

quantities), with
coefficients that depend on the oscillator g and its derivatives.
Secondly, the simplest (and most useful) special case is n = 1, whence (3)
reduces to
I[f, (0, 1)] ∼ −
∞
X
m=0
1
(−iω)
m+1

e
iωg(1)
g
′
(1)
f
m
(1) −
e
iωg(0)
g
′
(0)
f
m
(0)

. (4)
Thirdly, using an affine transformation, we can map S
n
to an arbitrary
simplex in R
n
. Applying an identical transformation to (3), we deduce that it
is valid for I[f, S], where S ⊂ R
n
is any simplex.
Fourthly, the boundary of S is itself composed of n + 1 simplices in R
n−1
.
Because of the nonresonance condition, the gradient of the oscillator does not
vanish in any of these simplices and we can apply (3) therein: this expresses
I[f, S] as an asymptotic expansion over (n − 2)-dimensional simplices. We
continue with this procedure until we reach 0-dimensional simplices: the n + 1
vertices of the original simplex. Bearing in mind our first observation, we thus
deduce that
I[f, S] ∼
∞
X
m=0
1
(−iω)
m+n
Θ
m
[f], (5)

Frequently asked questions

Q1. What are the contributions in "Highly oscillatory quadrature: the story so far" ?

In this paper, Olver et al. show how to approximate I [ f, ( a, b ) with narrow wedges to obtain an asymptotic expansion. 

Q2. What are the future works in "Highly oscillatory quadrature: the story so far" ?

Their narrative underlies the importance of further research into quadrature methods for highly oscillatory integrals, in particular in the presence of critical points and when exact moments are unavailable. The underlying idea there, assuming that both f and g can be analytically extended to the complex plane, is to find a path from each endpoint of Ω = ( a, b ) to infinity alongside which g ( z ) − g ( a ) and g ( z ) − g ( b ), respectively, are real and negative. Because of exponential decay of the integrand, each integration can be accomplished by familiar Gauss–Laguerre quadrature and the outcome matches Filon-type and Levin-type methods in its asymptotic behaviour. The authors further note that in the presence of critical points there is a need to integrate also along paths joining them with z = ∞ in a fairly nontrivial manner. 

Q3. What is the natural approach to calculating moments?

The most natural approach is to take a leaf off finite-element theory, tessellate a polytope with simplices (taking care to respect nonresonance) and interpolate in each simplex with suitable polynomials. 

Q4. What is the main benefit of Levin-type methods?

One huge benefit of Levin-type methods is that they work easily on complicated domains and complicated oscillators for which Filon-type methods utterly fail. 

Q5. What is the first lesson to be drawn from the analysis?

The first and foremost lesson to be drawn from their analysis is that, once the authors can understand the mathematics of high oscillation, the authors gain access to a wide variety of effective and affordable algorithms. 

Q6. Why does the gradient of the oscillator not vanish in any of these simplices?

Because of the nonresonance condition, the gradient of the oscillator does not vanish in any of these simplices and the authors can apply (3) therein: this expresses I[f,S] as an asymptotic expansion over (n − 2)-dimensional simplices. 

Q7. What is the simplest way to expand a polytope?

As long as the nonresonance condition is maintained throughout the approximation of Ω by polytopes, their methods can be extended to this setting. 

Q8. Why do the authors have to construct the tessellation via a simplicial complex?

In other words, because of their construction of the tessellation via a simplicial complex, the contributions from neighbouring simplices cancel at internal vertices and each 

Q9. What is the natural way of approximating?

The simplest and most natural means of approximating (1) consists of a truncation of the asymptotic expansion (5) (replacing S by a polytope Ω). 

Q10. how do the authors map fm to an arbitrary simplex?

using an affine transformation, the authors can map Sn to an arbitrary simplex in Rn. Applying an identical transformation to (3), the authors deduce that it is valid for I[f,S], where S ⊂ Rn is any simplex. 

Q11. What is the obvious case when is a polynomial?

In the most obvious case when ϕ is a polynomial, this is equivalent to the explicit computability of relevant moments of the oscillator g,µi(ω) =∫Ωxieiωg(x)dS, xi = xi11 · · ·xinn , i ∈ Zn+