What's the state-of-the-art in highly oscillatory integral computation?

Question

What's the state-of-the-art in the approximation of highly oscillatory integrals in both one dimension and higher dimensions to arbitrary precision?

Answer 1

I'm not entirely familiar with what's now done for cubatures (multidimensional integration), so I'll restrict myself to quadrature formulae.

There are a number of effective methods for the quadrature of oscillatory integrals. There are methods suited for finite oscillatory integrals, and there are methods for infinite oscillatory integrals.

For infinite oscillatory integrals, two of the more effective methods used are Longman's method and the modified double exponential quadrature due to Ooura and Mori. (But see also these two papers by Arieh Iserles.)

Longman's method relies on converting the oscillatory integral into an alternating series by splitting the integration interval, and then summing the alternating series with a sequence transformation method. For instance, when integrating an oscillatory integral of the form

$$\int_0^\infty f(t)\sin\,t\mathrm dt$$

one converts this into the alternating sum

$$\sum_{k=0}^\infty \int_{k\pi}^{(k+1)\pi} f(t)\sin\,t\mathrm dt$$

The terms of this alternating sum are computed with some quadrature method like Romberg's scheme or Gaussian quadrature. Longman's original method used the Euler transformation, but modern implementations replace Euler with more powerful convergence acceleration methods like the Shanks transformation or the Levin transformation.

The double exponential quadrature method, on the other hand, makes a clever change of variables, and then uses the trapezoidal rule to numerically evaluate the transformed integral.

For finite oscillatory integrals, Piessens (one of the contributors of QUADPACK) and Branders, in two papers, detail a modification of Clenshaw-Curtis quadrature (that is, constructing an Chebyshev polynomial expansion of the nonoscillatory part of the integrand). Levin's method, on the other hand, uses a collocation method for the quadrature. (I am told there is now a more practical version of the old standby, Filon's method, but I've no experience with it.)

---

These are the methods I remember offhand; I'm sure I've forgotten other good methods for oscillatory integrals. I will edit this answer later if I remember them.

Answer 2

Besides "multidimensional vs. single-dimensional" and "finite range vs. infinite range", an important categorization for methods is "one specific type of oscillator (usually Fourier-type: $\sin(t)$, $\exp(it)$, etc, or Bessel-type: $J_0(t)$, etc.) vs. more general oscillator ($\exp(i g(t))$ or even more general oscillators $w(t)$)".

At first, oscillatory integration methods focused on specific oscillators. As J. M. said, prominent ones include Filon's method and the Clenshaw-Curtis method (these two are closely related) for finite range integrals, and series extrapolation based methods and the double-exponential method of Ooura and Mori for infinite range integrals.

More recently, some general methods have been found. Two examples:

1. Levin's collocation-based method for any $\exp(i g(t))$ (Levin 1982), or later for any oscillator $w(t)$ defined by a linear ODE (Levin 1996 as linked by J. M.). Mathematica uses Levin's method for integrals not covered by the more specialized rules.

2. Huybrechs and Vandewalle's method based on analytic continuation along a complex path where the integrand is non-oscillatory (Huybrechs and Vandewalle 2006).

No distinction is necessary between methods for finite and infinite range integrals for the more general methods, since a compactifying transformation can be applied to an infinite range integral, leading to a finite range oscillatory integral that can still be addressed with the general method, albeit with a different oscillator.

Levin's method can be extended to multiple dimensions by iterating over the dimensions and other ways, but as far as I know all the methods described in literature so far have sample points that are an outer product of the one-dimensional sample points or some other thing that grows exponentially with dimension, so it rapidly gets out of hand. I'm not aware of more efficient methods for high dimensions; if any could be found that sample on a sparse grid in high dimensions it would be useful in applications.

Creating automatic routines for the more general methods may be difficult in most programming languages (C, Python, Fortran, etc) in which you would normally expect to program your integrand as a function/routine and pass it to the integrator routine, because the more general methods need to know the structure of the integrand (which parts look oscillatory, what type of oscillator, etc) and can't treat it as a "black box".

Answer 3

You could also check the work of Marnix Van Daele and co-authors. See for example this and this.