Numerical evaluation of oscillatory integrals via automated steepest descent contour deformation

Dave Hewett (University College London)

Abstract. Oscillatory integrals arise in many areas of mathematics, science and technology, including Fourier analysis, special function theory, signal processing, quantum mechanics, and acoustic and electromagnetic wave propagation. By an oscillatory integral we mean an integral over some subset of R^n in which the integrand undergoes a large number of oscillations. Typically, the integrand takes the form of a non-oscillatory prefactor multiplied by an oscillating exponential of the form exp(ikg(x)), for some “phase function” g(x) governing the oscillations. Such integrals can be prohibitively expensive to evaluate numerically using a conventional interpolatory quadrature method (e.g. Gauss quadrature, trapezoidal rule etc), because one needs a very large number of quadrature points to properly capture the oscillations in the integrand. Thankfully, a number of alternative quadrature approaches have been developed to bypass this problem. Our attention in this talk will be on one-dimensional oscillatory integrals, and the numerical steepest descent (NSD) method. The basic idea of NSD is to deform the integration contour into the complex plane onto a contour on which the integrand is non-oscillatory and exponentially decaying, so that conventional quadrature rules can be applied. However, unless the phase function governing the oscillation is particularly simple, the application of NSD requires a significant amount of a priori analysis and expert user input, to determine the appropriate contour deformation, and to deal with the non-uniformity in the accuracy of standard quadrature techniques associated with the coalescence of stationary points (zeros of the derivative of the phase function) with each other, or with the endpoints of the original integration contour. In this talk we present a novel algorithm for the numerical evaluation of oscillatory integrals with general polynomial phase functions, which automates the contour deformation process and avoids the difficulties typically encountered with coalescing stationary points and endpoints. The inputs to the algorithm are simply the phase and amplitude functions, the endpoints and orientation of the original integration contour, and a small number of numerical parameters. As an application, we use our algorithm to evaluate cuspoid canonical integrals from scattering theory. A Matlab implementation of the algorithm is available online and is called PathFinder.

This is joint work with Andrew Gibbs (UCL) and Daan Huybrechs (KU Leuven)