figuresection tablesection

#  Numerical evaluation of oscillatory integrals via automated steepest
descent contour deformation

A. Gibbs11footnotemark: 1, D. P. Hewett11footnotemark: 1, D.
Huybrechs22footnotemark: 2,
11footnotemark: 1 Department of Mathematics, University College London,
London, UK
22footnotemark: 2 Department of Computer Science, KU Leuven, Leuven, Belgium

###### Abstract

Steepest descent methods combining complex contour deformation with numerical
quadrature provide an efficient and accurate approach for the evaluation of
highly oscillatory integrals. However, unless the phase function governing the
oscillation is particularly simple, their application 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 (saddle points) with each other, or with the endpoints of
the original integration contour. In this paper 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. By a series of numerical
experiments we demonstrate that the algorithm is accurate and efficient over a
large range of frequencies, even for examples with a large number of
coalescing stationary points and with endpoints at infinity. As a particular
application, we use our algorithm to evaluate cuspoid canonical integrals from
scattering theory. A Matlab implementation of the algorithm is made available
and is called PathFinder.

##  1 Introduction

In this paper we consider numerical evaluation of the integral

$I=\int_{\Gamma}f(z)\mathrm{e}^{\textrm{i}\omega g(z)}\operatorname{d}\\!{z},$
(1)

where $\Gamma$ is a contour in , possibly starting and/or ending at infinity,
$f$ and $g$ are functions of a complex variable, and $\omega>0$ is a frequency
parameter. Such integrals arise in numerous application areas, particularly in
wave phenomena and quantum mechanics, and are generally challenging to
evaluate numerically, especially when $\omega$ is large, because the presence
of the exponential factor $\mathrm{e}^{\textrm{i}\omega g(z)}$ means that the
integrand may undergo rapid oscillations and/or significant variations in
amplitude along the integration contour.

When $f$ and $g$ are analytic, Cauchy’s theorem provides the possibility of
deforming the integration contour so as to make numerical evaluation easier.
This is the basis of _steepest descent (SD) methods_ , in which one aims to
deform $\Gamma$ onto a contour, or, more typically, a union of contours, which
we term the _steepest descent (SD) deformation_ , on which $\Re[g(z)]$ is
constant, so that the exponential factor $\mathrm{e}^{\textrm{i}\omega g(z)}$
is no longer oscillatory. By the Cauchy-Riemann equations, these contours
coincide with the steepest descent curves of $-\Im[g(z)]$, and they connect
endpoints of the original integration contour, valleys at infinity (sectors in
which the integrand decays rapidly as $|z|\to\infty$), and stationary points
of $g$, which are points $\xi\in\C$ at which $g^{\prime}(\xi)=0$.
111Stationary points are often referred to as “saddle points” because they are
saddle points of the functions $\Im[g(z)]$ and $\Re[g(z)]$, which cannot
possess local maxima or minima by the maximum modulus principle. Along each SD
contour, away from stationary points the integrand typically decays
exponentially, with the rate of decay increasing with increasing $\omega$, and
as $\omega\to\infty$ the value of the integral is dominated by local
contributions close to the endpoints of $\Gamma$ and any stationary points
traversed by the SD deformation. In the _asymptotic_ steepest descent method
(described e.g. in [3, 19]), one exploits this to obtain an asymptotic
expansion for the integral, valid as $\omega\to\infty$, by performing a local
Taylor expansion of the integrand around the endpoints and relevant stationary
points, and reducing the local integrals along the SD contours to simpler
integrals that can be expressed in terms of known special functions.

In the _numerical_ steepest descent (NSD) method (described e.g. in [7, §5])
one evaluates the integrals along the SD contours numerically. This involves
numerically tracing an appropriate segment of each SD contour in the SD
deformation and applying suitable numerical quadrature rules to evaluate the
associated contributions to the integral.

In principle, NSD is a highly accurate and efficient method for evaluating
integrals of the form (1) for moderate or large $\omega$. Indeed, under
appropriate assumptions, the NSD method outputs approximations which, for a
fixed number of quadrature points $N$, are asymptotic to (1) as
$\omega\to\infty$, with the asymptotic accuracy improving with increasing $N$
(see, e.g., [7, Thm 5.7]). Furthermore, if $f$ and $g$ are sufficiently well
behaved it can also be the case that the NSD approximation converges to (1) as
$N\to\infty$, for fixed $\omega>0$, with a cost that remains bounded as
$\omega\to\infty$.

In practice, however, applying the NSD method to an integral of the form (1)
often requires significant expert user input. This is because:

  1. (P1)

Determining the SD contour deformation corresponding to a given $g$ and
$\Gamma$ requires careful a priori analysis.

  2. (P2)

Parametrizing SD contours from or near stationary points, and evaluating
integrals along them, is fraught with numerical difficulties, especially when
stationary points are close to other stationary points or endpoints of
$\Gamma$.

The issues described in (P1) and (P2) are particularly troublesome when one
wishes to evaluate (1) for multiple instances of a phase function
$g(z)=g(z,\mathbf{c})$ depending on a set of parameters $\mathbf{c}\in{}^{q}$.
This is because, firstly, the number and location of the stationary points,
and the nature of the SD deformation, have to be determined for each different
value of $\mathbf{c}$, and, secondly, stationary points may coalesce as
$\mathbf{c}$ approaches certain regions in parameter space, leading to a non-
uniformity in the accuracy of the resulting NSD approximations.

The problem of stationary point coalescence in the context of NSD was studied
in detail in [10] in the special case of the cubic phase function
$g(z,c)=\frac{z^{3}}{3}-cz$, for $c\in\C$, which for $c\neq 0$ has a pair of
order one stationary points which coalesce as $c\to 0$ (at $z=0$) into a
single stationary point of order two for $c=0$.222The order of a stationary
point $\xi$ is the multiplicity of $\xi$ as a root of $g^{\prime}$. In this
case, the SD deformation and contour parametrization were carried out manually
by analytically inverting the phase (illustrating (P1)), but the resulting
integrals were found to be nearly singular for small $c$, leading to poor
accuracy of standard NSD approximations (illustrating (P2)). It was shown in
[10] how to construct a family of non-standard quadrature rules for this
integral which perform uniformly well for $c\approx 0$ using complex-valued
Gaussian quadrature, producing quadrature nodes that in general lie off the SD
deformation. In principle, similar rules could be developed for more
complicated coalescences involving higher order stationary points and/or
endpoints of $\Gamma$. However, for each type of coalescence a bespoke
quadrature rule would have to be developed, and a general catalogue of such
rules is not yet available in the literature.

In contrast to [10], our aim is not to develop an optimized method for a
specific instance of (1), but rather to present a relatively simple algorithm
that can evaluate (1) accurately, for a general class of $f$ and $g$, without
the need for expert user input or a priori analysis, even in the case of
coalescing stationary points, thus addressing problems (P1) and (P2). Our
specific focus in this paper is on the case where $f$ is entire and $g$ is a
polynomial. The extension of our approach to more general cases where $f$
and/or $g$ have pole or branch point singularities is the subject of ongoing
research. Necessarily, in aiming for generality and robustness we will
sacrifice some efficiency. However, our method is designed to be rapidly
convergent as $N\to\infty$ with approximately $\omega$-independent cost, and
the fact that this is realised in practice is demonstrated by extensive
numerical experiments in §5.

Our algorithm follows the basic principles of NSD, combining complex contour
deformation with numerical quadrature. However, in contrast to standard NSD
our algorithm does not trace SD contours directly from stationary points.
Instead, stationary points are enclosed in a bounded “non-oscillatory region”
within which the integrand is guaranteed to undergo at most a fixed number of
oscillations. The original contour $\Gamma$ is replaced by a “quasi-SD
deformation” comprising a union of straight-line contours in the non-
oscillatory region, for which numerical quadrature is straightforward, and SD
contours outside the non-oscillatory region, on which standard NSD quadrature
techniques can be applied. By excluding a neighbourhood of the stationary
points from the region in which SD contours are traced, our algorithm avoids
the problems mentioned in (P2) associated with stationary-point/stationary-
point and/or stationary-point/endpoint coalescence. This not only
“uniformizes” the accuracy of our algorithm compared to standard NSD, but it
also enables us to tackle the problem (P1) by automating the contour
deformation step. For the latter, we first perform low accuracy SD contour
tracing outside the non-oscillatory region to build a graph describing the
global connections (via SD contours) between the endpoints of $\Gamma$, the
different components of the non-oscillatory region, and the valleys at
infinity, and then determine the quasi-SD deformation using a shortest path
algorithm, before refining the accuracy of the SD contour tracing at the
quadrature stage.

One other problem with standard NSD is that it typically degenerates as
$\omega\to 0$, because the quadrature points diverge to infinity [7, §5.2.4].
This issue has been addressed in the special case $g(z)=z$ for bounded
$\Gamma$ in [2, 6]; however, it remains an open problem for general $g(z)$.
Our algorithm is well-behaved in the limit as $\omega\to 0$ for general
polynomial $g(z)$, since it reduces to standard non-oscillatory quadrature for
sufficiently small $\omega$ for any bounded $\Gamma$.

Our algorithm is implemented in the open-source Matlab code “PathFinder”,
available at github.com/AndrewGibbs/PathFinder [8]. The basic user input to
the code is a function handle for the amplitude f, the coefficients of the
polynomial phase g, endpoints a and b (complex numbers, or angles in the case
of infinite endpoints), the frequency parameter omega, and a parameter N
controlling the number of quadrature points to be used. Approximating the
integral (1) using PathFinder can be done with the following Matlab command:

PathFinder(a,b,f,g,omega,N,'infcontour',[A B]) (2)

Here 'infcontour' is an optional input for which the user should supply a
Boolean array [A B] (whose default value is [false false]) such that A
(respectively B) is true if the endpoint a (resp. b) is infinite and false if
it is finite. Examples of PathFinder code will be given in §5. Advanced users
can also adjust a small number of other tuning parameters, whose role will be
discussed during the presentation of our algorithm.

An outline of the paper is as follows. In §2 we provide a detailed description
of our algorithm, first presenting an overview of the main steps, and then
providing details of how each step is realised in PathFinder. In §3 we present
some theoretical results underpinning our approach. In §4 we discuss some
further implementation details, and in §5 we exhibit numerical results
demonstrating the performance of our algorithm on a range of challenging
integrals with large $\omega$ and complicated stationary point configurations.

We end this introduction by remarking that integrals with coalescing
stationary points are of fundamental importance in numerous applications,
including the study of optics and high frequency (short wavelength) acoustics,
where they describe the wave field in the vicinity of geometrical
singularities (or “catastrophes”) in the classical ray-theoretic framework,
Kelvin’s celebrated ship-wave problem, and the theory of molecular collisions
in quantum mechanics and theoretical chemistry. A catalogue of such integrals,
along with links to relevant literature, can be found in [1, §36]. In §5.5 we
show how PathFinder can be applied to accurately calculate these types of
integrals.

##  2 Algorithm description

In this section we present our algorithm for the numerical approximation of
(1) when $f$ is entire333When $\Gamma$ is infinite we additionally implicitly
assume that $f$ is sufficiently well-behaved at infinity for the integral (1)
to converge. Note that in many cases of interest, numerical evaluation of the
integral to high accuracy after path deformation only requires $f$ to be
analytic in a small (and shrinking with increasing $\omega$) neighbourhood of
the stationary points. and $g$ is a polynomial.

We start with some definitions and basic facts. Let

$g(z)=\sum_{j=0}^{J}\alpha_{j}z^{j},$ (3)

for some $J\in\N$, $J\geq 1$, and $\alpha_{j}\in\C$, $j=0,\ldots,J$, with
$\alpha_{J}\neq 0$. Then $g$ has at most $J-1$ stationary points, which are
the solutions of

$\displaystyle g^{\prime}(z)=\sum_{j=1}^{J}j\alpha_{j}z^{j-1}=0.$ (4)

We denote the set of all stationary points by $\mathcal{P}_{\rm stat}$. We
define the valleys at infinity to be the sectors of angular width $\pi/J$
centred on the angles

$\mathcal{V}:={\left\\{\frac{(2(m-1)+1/2)\pi-\arg{(\alpha_{J})}}{J}:\quad
m=1,\ldots,J\right\\}}.$ (5)

These have the property that if $z=r\mathrm{e}^{\textrm{i}\theta}$ with
$\theta\in(v-\pi/(2J),v+\pi/(2J))$ for some $v\in\mathcal{V}$ then
$\mathrm{e}^{\textrm{i}\omega g(z)}\to 0$ as $r\to\infty$. For each
$\eta\in\C\setminus\mathcal{P}_{\rm stat}$ there exists a unique SD contour
$\gamma_{\eta}$ beginning at $\eta$ and ending either at a stationary point
$\xi\in\mathcal{P}_{\rm stat}$ or at a valley $v\in\mathcal{V}$, on which $\Re
g(z)=\Re g(\eta)$ for $z\in\gamma_{\eta}$ (see, e.g., [4]).

We let $\mathcal{P}_{\rm endp}$ denote the set of finite endpoints of the
integration contour $\Gamma$, which could have zero, one or two elements. We
assume for now that any infinite endpoint of $\Gamma$ is at one of the valleys
$v\in\mathcal{V}$; see §4.4 for extensions.

We now provide a high-level overview of our algorithm. The following steps
will be explained in more detail in sections 2.1-2.6.

  1. 1.

Compute the set of stationary points $\mathcal{P}_{\rm stat}$ (the solutions
of (4)).

  2. 2.

For each $\xi\in\mathcal{P}_{\rm stat}$, select a radius $r_{\xi}>0$ for which
the function $\textrm{e}^{\textrm{i}\omega g(z)}$ is considered “non-
oscillatory” on the closed ball $\Omega_{\xi}$ of radius $r_{\xi}$ centred at
$\xi$. These balls may overlap. However, if two balls overlap significantly,
indicating near coalescence, one of the stationary points (along with its
associated ball) is removed from the set $\mathcal{P}_{\rm stat}$. This
removal process continues recursively until no pair of balls is judged to
overlap too much.

We call $\\{\Omega_{\xi}\\}_{\xi\in\mathcal{P}_{\rm stat}}$ the _non-
oscillatory balls_ , and their union

$\displaystyle\Omega:=\bigcup_{\xi\in\mathcal{P}_{\rm stat}}\Omega_{\xi}$ (6)

the _non-oscillatory region_.

  3. 3.

Find the local minima of $|\textrm{e}^{\textrm{i}\omega g(z)}|$ on the
boundary of the non-oscillatory region $\Omega$. We call these points _exits_
, and denote by $\mathcal{P}_{\rm exit}$ the set of all exits.

  4. 4.

For each $\eta\in\mathcal{P}_{\rm exit}\cup(\mathcal{P}_{\rm
endp}\setminus\Omega)$, trace the SD contour $\gamma_{\eta}$ from $\eta$, and
determine whether

     * •

(i) $\gamma_{\eta}$ enters $\Omega$ at some point
$z\in\partial\Omega\setminus\\{\eta\\}$, or

     * •

(ii) $\gamma_{\eta}$ converges towards a valley $v\in\mathcal{V}$ without
entering $\Omega$.

We call points $z\in\partial\Omega$ determined in case (i) entrances, and
denote by $\mathcal{P}_{\rm entr}$ the set of all entrances.

  5. 5.

Construct a graph $G$ with a vertex for each of the elements of
$\mathcal{P}_{\rm stat}$, $\mathcal{P}_{\rm endp}$, $\mathcal{P}_{\rm exit}$,
$\mathcal{P}_{\rm entr}$ and $\mathcal{V}$. Add edges between the vertices of
$G$ as follows:

     * •

For each $\xi\in\mathcal{P}_{\rm stat}$, add an edge between each pair of
elements of $(\mathcal{P}_{\rm stat}\cup\mathcal{P}_{\rm
endp}\cup\mathcal{P}_{\rm exit}\cup\mathcal{P}_{\rm entr})\cap\Omega_{\xi}$.

     * •

For each pair $\xi,\xi^{\prime}\in\mathcal{P}_{\rm stat}$,
$\xi\neq\xi^{\prime}$, for which
$\Omega_{\xi}\cap\Omega_{\xi^{\prime}}\neq\emptyset$, add an edge between
$\xi$ and $\xi^{\prime}$, if not already added in the previous step.

     * •

For each $\eta\in\mathcal{P}_{\rm exit}\cup(\mathcal{P}_{\rm
endp}\setminus\Omega)$, add an edge between $\eta$ and the entrance
$z\in\mathcal{P}_{\rm entr}$ or the valley $v\in\mathcal{V}$ to which the SD
contour $\gamma_{\eta}$ leads.

Find the shortest path (in the graph-theoretic sense) between the vertices
corresponding to the endpoints of $\Gamma$.

  6. 6.

Generate quadrature nodes and weights for the evaluation of each of the
contour integrals corresponding to the edges in the shortest path. For an edge
between two points in the non-oscillatory region, use a straight-line contour.
For an edge between an exit or an endpoint of $\Gamma$ to an entrance or a
valley, use a refined version of the SD contour traced in step 4. The union of
all the contours corresponding to the edges of the shortest path defines the
“quasi-SD deformation” of the original integration contour.

Finally, use the quadrature nodes and weights to approximate the integrals
over the contours in the quasi-SD deformation and sum them to obtain an
approximation of the original integral (1).

In Figures 1 and 2 we illustrate the outcome of the above steps for the
particular choice of phase function

$\begin{split}g(z)=&\frac{z^{7}}{7}+z^{6}\,\left(\frac{7}{20}+\frac{13}{30}{}\mathrm{i}\right)+z^{5}\,\left(-\frac{1047}{2000}+\frac{543}{1000}{}\mathrm{i}\right)+z^{4}\,\left(-\frac{4409}{8000}-\frac{5077}{8000}{}\mathrm{i}\right)\\\
&+z^{3}\,\left(\frac{711}{2000}-\frac{4441}{6000}{}\mathrm{i}\right)+z^{2}\,\left(\frac{237}{800}-\frac{207}{800}{}\mathrm{i}\right)+z\,\left(\frac{63}{1000}-\frac{77}{2000}{}\mathrm{i}\right)\end{split}$
(7)

and the parameters $\omega=40$, $a=-1.5$, $b=2$, $N=10$, using the default
parameter set for PathFinder (see Table 1). For this choice of $g$ there is
one order 2 stationary point and 4 order one stationary points. In Figure 1 we
plot these stationary points, along with their non-oscillatory balls, and the
SD contours traced from the exits. Such plots can be generated in PathFinder
by adding the optional 'plot' flag. The ball centred at the stationary point
$\xi=-\textrm{i}$ contains two entrances, reached by SD contours from the
balls above. In Figure 2 we plot the graph $G$, using the optional PathFinder
input flag 'plot graph'. This graph, in addition to edges corresponding to the
SD contours shown in Figure 1, contains edges corresponding to contours
between points in the non-oscillatory region, including connections within the
two overlapping balls. The shortest path between $a$ and $b$, which is
highlighted with thick lines in Figure 2, corresponds to the quasi-SD
deformation, the integral over which is equal to (1) by Cauchy’s Theorem. The
integral is discretised using $N$ quadrature points on each contour in the
quasi-SD deformation that makes a non-negligible contribution to the integral
(see §2.6) - these points are plotted in Figure 1 in red.

The process of computing all the SD contours and the selection of a subset
thereof via the shortest path algorithm addresses problem (P1). Surrounding
stationary points by balls, and only tracing SD contours outside the balls,
means that we avoid having to determine the local structure of the SD contours
and compute integrals along them near stationary points, addressing problem
(P2).

(a) SD contours and quasi-SD deformation

(b) Zoom of boxed region near the origin

Figure 1: Output of algorithm when applied with phase (7), $\omega=40$,
$a=-1.5$, $b=2$, $N=10$, and the default parameter set (see Table 1). Here we
observe stationary points (black stars) surrounded by balls (grey), SD
contours traced from exits and finite endpoints (black lines), and quadrature
points allocated along the appropriate contours in the quasi-SD deformation
(red points). The “region of no return” (see §3.2) around the valleys at
infinity is also shaded grey.

(a) Graph $G$

(b) Zoom of the graph $G$

Figure 2: The graph $G$ corresponding to the problem considered in Figure 1.
The thick line represents the shortest path between the endpoints, which in
this case are both finite. The balls (shaded grey) are included for ease of
comparison with Figure 1. The lower figure zooms in on the centre of the upper
figure, showing the multiple edges that are constructed inside the balls.

In the following subsections we provide further details of how we carry out
the steps outlined above in PathFinder.

###  2.1 Step 1 - Computing stationary points

Computing the stationary points of $g$ (the roots of $g^{\prime}(z)$) requires
us to find the complex roots of the polynomial (4). In our implementation we
compute stationary points using the Matlab `roots` command, which applies a
companion matrix approach. We note that obtaining highly accurate values for
the positions of stationary points is not critical to our algorithm, since the
stationary points are enclosed in the non-oscillatory region and we never
trace SD contours from them. Indeed, the difficulty in distinguishing
numerically between multiple roots and roots of higher order contributes to
the motivation for considering such non-oscillatory regions.

###  2.2 Step 2 - Defining the non-oscillatory region

The non-oscillatory region $\Omega$ was defined in (6) to be a union of balls
centred at the elements of $\mathcal{P}_{\rm stat}$. We choose the radii of
the balls as follows. First fix some user-defined constant $C_{\rm ball}>0$.
Then, given $\xi\in\mathcal{P}_{\rm stat}$, define

$\displaystyle r_{\xi}:=\max\\{r>0:|z-\xi|\leq
r\Rightarrow\omega|g(z)-g(\xi)|\leq C_{\rm ball}\\}.$ (8)

This definition enforces an upper bound on the number of oscillations within
each ball. Accordingly, the region $\Omega$ shrinks to the stationary points
as $\omega\to\infty$ and expands to fill the whole complex plane as $\omega\to
0$.

In our implementation we approximate $r_{\xi}$ numerically as follows. Let
$N_{\rm ball}\in\N$ be a user-defined parameter. For each
$n\in\\{1,\ldots,N_{\rm ball}\\}$ we consider the ray
$\\{z=\xi+r\textrm{e}^{\textrm{i}2\pi n/N_{\rm ball}},\,r>0\\}$, and compute
the smallest positive root $r_{n}>0$ of the function
$u_{n}(r):=\omega^{2}|g(\xi+r\textrm{e}^{\textrm{i}2\pi n/N_{\rm
ball}})-g(\xi)|^{2}-C_{\rm ball}^{2}$, which is a polynomial in $r$ of degree
$2J$. For this root-finding problem we use the Matlab `roots` command; in case
this command produces no positive real roots (because of stability issues) we
resort to a bisection approach instead. We then take as our approximation to
$r_{\xi}$ the positive number $\max_{n\in\\{1,\ldots,N_{\rm ball}\\}}r_{n}$.

When elements of $\mathcal{P}_{\rm stat}$ are close it is natural to
amalgamate their respective non-oscillatory balls. To do this systematically
we adopt an iterative approach. Let $\delta_{\rm ball}>0$ be a user-defined
parameter.

  * •

For each pair $\xi_{1},\xi_{2}\in\mathcal{P}_{\rm stat}$ compute

$d_{\xi_{1},\xi_{2}}:=|\xi_{1}-\xi_{2}|/\max(r_{\xi_{1}},r_{\xi_{2}}).$

  * •

If $\min_{\xi_{1},\xi_{2}}d_{\xi_{1},\xi_{2}}<\delta_{\rm ball}$ let
$\xi_{1},\xi_{2}$ be a pair realising the minimum. Remove from
$\mathcal{P}_{\rm stat}$ whichever of $\xi_{1},\xi_{2}$ has the smaller
associated ball radius (or choose arbitrarily between them if
$r_{\xi_{1}}=r_{\xi_{2}}$).

  * •

Repeat the previous step until either
$\min_{\xi_{1},\xi_{2}}d_{\xi_{1},\xi_{2}}\geq\delta_{\rm ball}$, or there is
only one element of $\mathcal{P}_{\rm stat}$ remaining.

###  2.3 Step 3 - Determining the exits

The exits associated with each $\xi\in\mathcal{P}_{\rm stat}$ are defined to
be the local minima on
$\partial\Omega_{\xi}\setminus\bigcup_{\xi^{\prime}\in\mathcal{P}_{\rm
stat},\xi^{\prime}\neq\xi}\Omega_{\xi^{\prime}}^{\circ}$ of the function
$|\textrm{e}^{\textrm{i}\omega g(z)}|$, equivalently of the function $-\Im
g(z)$.

For each $\xi\in\mathcal{P}_{\rm stat}$ the function $-\Im g(z)$ restricted to
the $\partial\Omega_{\xi}$ is a trigonometric polynomial. Using this fact, in
our implementation we determine the local minima of $-\Im g(z)$ on
$\partial\Omega_{\xi}$ by finding the roots of the derivative of $-\Im g(z)$
in the angular direction (which is also a trigonometric polynomial) by the
companion matrix approach of [5, §2.2], and keep only the real roots
corresponding to local minima. We discard all those minima corresponding to
points inside $\bigcup_{\xi^{\prime}\in\mathcal{P}_{\rm
stat},\xi^{\prime}\neq\xi}\Omega_{\xi^{\prime}}^{\circ}$, and add the
remaining minima to the set $\mathcal{P}_{\rm exit}$.

###  2.4 Step 4 - Tracing the SD contours

Given $\eta\in\mathcal{P}_{\rm exit}\cup(\mathcal{P}_{\rm
endp}\setminus\Omega)$, the SD contour $\gamma_{\eta}$ beginning at $\eta$ is
the unique curve on which $\Re g(z)$ is constant, with $-\Im g(z)$ decreasing
along $\gamma_{\eta}$. It can be parametrized in terms of a parameter $p\geq
0$ as $z=h_{\eta}(p)$, where $h_{\eta}(p)$ is defined implicitly by

$\displaystyle g(h_{\eta}(p))=g(\eta)+\textrm{i}p,\qquad h_{\eta}(0)=\eta.$
(9)

Differentiating (9) with respect to $p$ gives

$\displaystyle
h_{\eta}^{\prime}(p)=\frac{\textrm{i}}{g^{\prime}(h_{\eta}(p))}=:F(h_{\eta}(p)),\qquad
h_{\eta}(0)=\eta,$ (10)

which is a first order ODE initial value problem for $h_{\eta}(p)$. By solving
(10) numerically one can trace the contour $\gamma_{\eta}$ until it either (i)
enters the non-oscillatory region $\Omega$, or (ii) one can be sure that it
will tend to a valley $v\in\mathcal{V}$, without entering $\Omega$. For (ii)
we appeal to the theoretical result in Theorem 3.3, which provides a “region
of no return” $R_{v}$ associated with each valley $v\in\mathcal{V}$ for which
it is guaranteed that if an SD contour enters $R_{v}$ it will never leave
$R_{v}$, and will converge to $v$.

Staying away from stationary points ensures that the factor $1/g^{\prime}$ in
the right-hand side of (10) does not get too large.

In our implementation we trace the SD contour using a predictor-corrector
approach, combining a forward Euler step for (10) and a Newton iteration for
(9), to generate approximations $h_{\eta}^{(n)}\approx h_{\eta}(p_{n})$ on a
mesh $0=p_{0}<p_{1}<p_{2}<\ldots<p_{n_{\rm max}}$, where the total number of
steps $n_{\rm max}$ is determined as part of the algorithm, as discussed
below.

As the initial value we take $h_{\eta}^{(0)}=\eta$. Then, given
$h_{\eta}^{(n)}$, to compute $h_{\eta}^{(n+1)}$ we first apply a forward Euler
step for the ODE (10), with adaptive step length

$p_{n+1}-p_{n}=\delta_{\rm
ODE}\min\left(2\frac{|g^{\prime}(h_{\eta}^{(n)})|^{2}}{|g^{\prime\prime}(h_{\eta}^{(n)})|},|g^{\prime}(h_{\eta}^{(n)})|\mathop{\rm
dist}(h_{\eta}^{(n)},\mathcal{P}_{\rm stat})\right),$

where $\delta_{\rm ODE}\in(0,1)$ is a user-specified parameter. The first
argument of the minimum is included to ensure stability of the solver - note
that
$F^{\prime}(h)=-\frac{\textrm{i}g^{\prime\prime}(h)}{(g^{\prime}(h))^{2}}$ and
we might expect instability if the local step length were as large as
$2/|F^{\prime}(h)|=2\frac{|g^{\prime}(h)|^{2}}{|g^{\prime\prime}(h)|}$. The
second argument is included to ensure that the solver “slows down” as it
approaches the non-oscillatory region $\Omega$, so that we can detect whether
$\gamma_{\eta}$ enters $\Omega$ or not. To ensure that
$|h_{\eta}^{(n+1)}-h_{\eta}^{(n)}|\leq\delta_{\rm ODE}d$, where
$d:=\mathop{\rm dist}(h_{\eta}^{(n)},\mathcal{P}_{\rm
stat})=\min_{\xi\in\mathcal{P}_{\rm stat}}|h_{\eta}^{(n)}-\xi|$, we require
that $p_{n+1}-p_{n}\leq\frac{\delta_{\rm
ODE}d}{|F(h_{\eta}^{(n)})|}=\delta_{\rm ODE}d|g^{\prime}(h_{\eta}^{(n)})|$.
This also ensures that $h_{\eta}^{(n+1)}$ remains far enough from
$\mathcal{P}_{\rm stat}$, so that (10) doesn’t get too large.

After each forward Euler step, we correct $h_{\eta}^{(n+1)}$ by running a
Newton iteration to enforce (9) (with $p=p_{n+1}$ fixed), until the Newton
step size
$|\frac{g(h_{\eta}^{(n+1)})-g(\eta)-\textrm{i}p_{n+1}}{g^{\prime}(h_{\eta}^{(n+1)})}|$
is smaller than $\delta_{\rm coarse}\mathop{\rm
dist}(h_{\eta}^{(n+1)},\mathcal{P}_{\rm stat})$, for some user-specified
tolerance $\delta_{\rm coarse}>0$.

We repeat this process for $n=0,1,2,\ldots$ until either

  * (i)

$h_{\eta}^{(n)}\in\Omega_{\xi}$ for some $\xi\in\mathcal{P}_{\rm stat}$, in
which case we add $z=h_{\eta}^{(n)}$ to the set $\mathcal{P}_{\rm entr}$ of
entrances. Note that in general the point $z=h_{\eta}^{(n)}$ will lie inside
$\Omega_{\xi}^{\circ}$ rather than on $\partial\Omega_{\xi}$, but will be
closer to $\partial\Omega_{\xi}$ the smaller $\delta_{\rm ODE}$ is;

  * or

  * (ii)

$h_{\eta}^{(n)}\in R_{v}$ for some $v\in\mathcal{V}$, in which case, by
Theorem 3.3, $\gamma_{\eta}$ converges to the valley $v$. Here the “region of
no return” $R_{v}$ is defined by

$\displaystyle R_{v}:=\\{z\in\C:|\arg z-v|_{2\pi}<\pi/(2J)\text{ and
}G(|z|,|\arg z-v|_{2\pi})>0\\},$ (11)

where

$\displaystyle|\theta|_{2\pi}:=\min_{m\in\Z}|\theta-2\pi m|,$ (12)

and, for $r>0$ and $\theta\in(0,\pi/(2J))$,

$\displaystyle
G(r,\theta):=J|\alpha_{J}|r^{J-1}\min\left(1/\sqrt{2},\cos{J\theta}\right)-\sum_{j=1}^{J-1}j|\alpha_{j}|r^{j-1}.$
(13)

For further explanation of the meaning of $R_{v}$ see §3.2 below.

A necessary condition for a point $z$ to lie in $R_{v}$ is that $|z|\geq
r_{*}$, where $r_{*}>0$ is the unique positive solution of the polynomial
equation $G(r_{*},\pi/(4J))=0$, i.e. the solution of

$\frac{J|\alpha_{J}|r_{*}^{J-1}}{\sqrt{2}}=\sum_{j=1}^{J-1}j|\alpha_{j}|r_{*}^{j-1}.$

Having found $r_{*}$ once and for all (using the Matlab roots commnad), to
check that a point $z$ lies in $R_{v}$ we first check that $|z|\geq r_{*}$. If
so, we then check that $|\arg z-v|_{2\pi}<\pi/(2J)$. If so, we then check that
$G(|z|,|\arg z-v|_{2\pi})>0$. The point of introducing $r_{*}$ is so that we
don’t compute $G(|z|,|\arg z-v|_{2\pi})$ unless absolutely necessary.

In either case, tracing of the SD contour stops and we set $n_{\rm max}=n$ for
this contour.

###  2.5 Step 5 - Finding the shortest path

The construction of the graph $G$ requires no further explanation. To find the
shortest path in $G$ between the endpoints of the original contour $\Gamma$ we
apply the standard Dijkstra shortest path algorithm [13, §10.6.2].

###  2.6 Step 6 - Evaluating the contour integrals

The quasi-SD contour deformation corresponding to the graph-theoretic shortest
path between the endpoints of $G$ calculated during step 5 involves integrals
over three types of contour:

  * Type 1:

Straight line contours between points in the non-oscillatory region;

  * Type 2:

Infinite SD contours from exits/endpoints to valleys;

  * Type 3:

Finite SD contours from exits/endpoints to entrances.

Some of these contours will make a larger contribution to the value of the
original integral (1) than others. It is natural to neglect contours that make
a negligibly small contribution. In our implementation, we only compute the
contribution from a contour $\gamma$ in the quasi-SD deformation if at least
one of the finite endpoints $\eta$ of $\gamma$ satisfies
$|\textrm{e}^{\textrm{i}\omega g(\eta)}|/M>\delta_{\rm quad}$, where
$\delta_{\rm quad}\geq 0$ is a small, user-specified parameter and

$M:=\max|\textrm{e}^{\textrm{i}\omega g(\xi)}|,$

where the maximum is taken over all $\xi\in\mathcal{P}_{\rm
stat}\cup\mathcal{P}_{\rm endp}\cup\mathcal{P}_{\rm exit}$ appearing in the
shortest path corresponding to the quasi-SD deformation.

In our implementation, for Type 1 contours we use Gauss-Legendre quadrature,
for Type 2 contours we use either Gauss-Laguerre quadrature (which is the
default choice in PathFinder) or truncated Gauss-Legendre quadrature, and for
Type 3 contours we use (possibly truncated) Gauss-Legendre quadrature, as
detailed below. By default our implementation uses the same number $N$ of
quadrature points on each contour in the quasi-SD deformation whose
contribution we compute, regardless of the type of integral (we comment on
this in §3.3.4). Accordingly, if $N_{\rm cont}$ is the number of these
contours then the total number of quadrature points used in the algorithm,
$N_{\rm tot}$, is given by

$\displaystyle N_{\rm tot}=NN_{\rm cont}.$ (14)

####  2.6.1 Evaluation of integrals over Type 1 contours

Let $z_{0},z_{1}\in\Omega$, and let $\gamma$ be the straight-line contour in
starting at $z_{0}$ and ending at $z_{1}$, parametrized by

$\displaystyle
z_{[z_{0},z_{1}]}(t)=\frac{1}{2}\Big{(}(z_{1}-z_{0})t+(z_{0}+z_{1})\Big{)},\qquad
t\in[-1,1].$ (15)

Given $N\in\N$, let $t_{m}^{\rm Leg}$ and $w_{m}^{\rm Leg}$, for
$m=1,\ldots,N$, denote the nodes and weights for standard $N$-point Gauss-
Legendre quadrature on $[-1,1]$. Our quadrature approximation to the integral
over $\gamma$ is then:

$\displaystyle\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z\approx\frac{z_{1}-z_{0}}{2}\sum_{m=1}^{N}w_{m}^{\rm
Leg}f(z_{[z_{0},z_{1}]}(t_{m}^{\rm Leg}))\textrm{e}^{\textrm{i}\omega
g(z_{[z_{0},z_{1}]}(t_{m}^{\rm Leg}))}.$ (16)

####  2.6.2 Evaluation of integrals over Type 2 contours

Let $\eta\in\mathcal{P}_{\rm exit}\cup(\mathcal{P}_{\rm endp}\setminus\Omega)$
be such that the SD contour $\gamma$ from $\eta$ leads to a valley.
Parametrizing $\gamma$ by (9), for $p\in[0,\infty)$, noting (10), and
rescaling $p=\tilde{p}/\omega$, we have

$\displaystyle\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z=\frac{\textrm{e}^{\textrm{i}\omega
g(\eta)}}{\omega}\int_{0}^{\infty}\tilde{f}(\tilde{p})\textrm{e}^{-\tilde{p}}\,\operatorname{d}\\!\tilde{p},$
(17)

where

$\tilde{f}(\tilde{p}):=f(h_{\eta}(\tilde{p}/\omega))h_{\eta}^{\prime}(\tilde{p}/\omega)=\textrm{i}\frac{f(h_{\eta}(\tilde{p}/\omega))}{g^{\prime}(h_{\eta}(\tilde{p}/\omega))}.$

By tracing contours outside of $\Omega$, the contours remain a positive
distance from $\mathcal{P}_{\rm stat}$. This ensures that $\tilde{f}$ is
analytic in a complex neighbourhood of $[0,\infty)$.

By default, PathFinder evaluates the integral on the right-hand side of (17)
by Gauss-Laguerre quadrature. Let $t_{m}^{\rm Lag}$ and $w_{m}^{\rm Lag}$, for
$n=1,\ldots,N$, denote the standard Gauss-Laguerre nodes and weights on
$[0,\infty)$. Our quadrature approximation to the integral over $\gamma$ is
then:

$\displaystyle\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z\approx\frac{\textrm{e}^{\textrm{i}\omega
g(\eta)}}{\omega}\sum_{m=1}^{N}w_{m}^{\rm Lag}\tilde{f}(t_{m}^{\rm Lag}).$
(18)

To evaluate $\tilde{f}(t_{m}^{\rm Lag})$ we need accurate computations of
$h_{\eta}(t_{m}^{\rm Lag}/\omega)$ for $m=1,\ldots,N$. For this, for each $m$
we run a Newton iteration on (9) with $p=t_{m}^{\rm Lag}/\omega$ fixed, until
the magnitude of the increment is smaller than a user-specified tolerance
$\delta_{\rm fine}>0$. Typically we take $\delta_{\rm fine}$ to be
considerably smaller than the tolerance $\delta_{\rm coarse}$ used in the
Newton iteration in step 4, since when carrying out quadrature we require
higher accuracy in our approximation of the SD contour than is required for
determining the global structure of the quasi-SD deformation in step 4. As the
initial guess for the Newton method we use a piecewise linear interpolant of
the points $\\{(p_{0},h_{\eta}^{(0)}),(p_{1},h_{\eta}^{(1)}),\ldots,(p_{n_{\rm
max}},h_{\eta}^{(n_{\rm max})})\\}$ computed in step 4, where $n_{\rm max}$
denotes the total number of steps taken in the ODE solve in step 4 before the
contour tracing algorithm terminated. If $p_{n_{\rm max}}<t_{N}^{\rm
Lag}/\omega$ then before running the Newton iteration we first need to restart
the contour tracing algorithm of step 4 to extend the SD contour until
$p_{n_{\rm max}}\geq t_{N}^{\rm Lag}/\omega$, so that there are points to
interpolate between.

As an alternative, one can evaluate the integral over a Type 2 contour using
truncated Gauss-Legendre quadrature, as suggested in [16]. To activate this
alternative in PathFinder one should add the optional input 'inf quad rule',
'legendre'. In this case we truncate the integral to

$\displaystyle\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z\approx\frac{\textrm{e}^{\textrm{i}\omega
g(\eta)}}{\omega}\int_{0}^{P}\tilde{f}(\tilde{p})\textrm{e}^{-\tilde{p}}\,\operatorname{d}\\!\tilde{p},$
(19)

for some $P>0$, then apply Gauss-Legendre quadrature on $[0,P]$, to obtain the
approximation

$\displaystyle\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z\approx\frac{P\textrm{e}^{\textrm{i}\omega
g(\eta)}}{2\omega}\sum_{m=1}^{N}w_{m}^{\rm Leg}\tilde{f}(z_{[0,P]}(t_{m}^{\rm
Leg}))\textrm{e}^{-z_{[0,P]}(t_{m}^{\rm Leg})},$ (20)

where we compute $h_{\eta}(z_{[0,P]}(t_{m}^{\rm Leg}/\omega))$ (which is
required for the evaluation of $\tilde{f}(z_{[0,P]}(t_{m}^{\rm Leg}))$) by the
same Newton iteration discussed above for $h_{\eta}(t_{m}^{\rm Lag}/\omega)$.
For the truncation point $P$ we take

$\displaystyle P=L,$ (21)

where

$\displaystyle L:=-\log\left(\delta_{\rm quad}M/|\textrm{e}^{\textrm{i}\omega
g(\eta)}|\right),$ (22)

which describes the point at which the magnitude of the exponential part of
the integrand drops below $\delta_{\rm quad}$ times its maximum value $M$ on
the quasi-SD deformation.

####  2.6.3 Evaluation of integrals over Type 3 contours

Let $\eta\in\mathcal{P}_{\rm exit}\cup(\mathcal{P}_{\rm endp}\setminus\Omega)$
be such that the SD contour $\gamma$ from $\eta$ leads to an entrance
$z\in\mathcal{P}_{\rm entr}$. In this case we apply (possibly truncated)
Gauss-Legendre quadrature as in formulas (19) and (20), but now with

$\displaystyle P=\min(p_{n_{\rm max}}/\omega,L),$ (23)

where $p_{n_{\rm max}}$ is defined as in §2.4 and $L$ is defined as in §2.6.2.

In the case where the minimum is attained by $p_{n_{\rm max}}/\omega$, so that
the whole contour is considered, a potential inconsistency arises, because the
application of the higher accuracy Newton iteration described in §2.6.2 for
the calculation of $h_{\eta}(z_{[0,P]}(t_{m}^{\rm Leg}/\omega))$ corresponds
implicitly to a slight shifting of the endpoint of the contour $\gamma$ away
from the entrance $z=h_{\eta}^{(n_{\rm max})}$ added to the graph $G$ in step
5. To avoid this inconsistency, in our implementation, in step 4, whenever the
contour tracing terminates in case (i), we run a Newton iteration on the final
point $h_{\eta}^{(n_{\rm max})}$ with the high accuracy tolerance $\delta_{\rm
fine}$, before adding it to the list of entrances $\mathcal{P}_{\rm entr}$.
Note that this may mean that $h_{\eta}^{(n_{\rm max})}$ lies very slightly
outside $\Omega$.

##  3 Theoretical results

In this section we collect some theoretical results that motivate the design
of our algorithm.

###  3.1 Removal of stationary points

In §2.2 we described our algorithm for removing stationary points from the set
$\mathcal{P}_{\rm stat}$ when they are close. When removing stationary points
and their associated non-oscillatory balls, we need to ensure that the removed
stationary points still lie inside one of the remaining non-oscillatory balls,
so that we don’t encounter any stationary points along the trajectory in our
ODE solve for the SD contour tracing (see the discussion in §3.3.2 below). In
this section we provide a sufficient condition on the parameter $\delta_{\rm
ball}$ for this to be guaranteed.

######  Proposition 3.1.

Suppose that in the removal algorithm of §2.2, $n$ stationary points have been
removed from $\mathcal{P}_{\rm stat}$. Then for any stationary point $\xi$
that was removed, there exists $\xi^{\prime}\in\mathcal{P}_{\rm stat}$ such
that $|\xi-\xi^{\prime}|\leq n\delta_{\rm ball}r_{\xi^{\prime}}$.

###### Proof.

We proceed by induction on $n$. The result is trivially true for $n=0$. Assume
that it is true after the removal of $n$ points, and suppose that the
$(n+1)$st point is now to be removed. Let $\xi_{1},\xi_{2}$ denote the pair of
points selected as realising $\min_{\xi_{1},\xi_{2}}d_{\xi_{1},\xi_{2}}$, and
without loss of generality suppose that $\xi_{2}$ is the point to be removed
(so that $r_{\xi_{1}}\geq r_{\xi_{2}}$). Then
$|\xi_{2}-\xi_{1}|\leq\delta_{\rm ball}r_{\xi_{1}}\leq(n+1)\delta_{\rm
ball}r_{\xi_{1}}$, so the claimed property holds for $\xi_{2}$. Furthermore,
by the inductive hypothesis, for each point $\xi$ previously removed, there
exists $\xi^{\prime}\in\mathcal{P}_{\rm stat}$ such that
$|\xi-\xi^{\prime}|\leq n\delta_{\rm ball}r_{\xi^{\prime}}$. If
$\xi^{\prime}\neq\xi_{2}$ then $\xi^{\prime}$ will still be present in
$\mathcal{P}_{\rm stat}$ after the removal of $\xi_{2}$, and
$|\xi-\xi^{\prime}|\leq n\delta_{\rm ball}r_{\xi^{\prime}}\leq(n+1)\delta_{\rm
ball}r_{\xi^{\prime}}$. On the other hand, if $\xi^{\prime}=\xi_{2}$ then
$\xi^{\prime}$ will not be present in $\mathcal{P}_{\rm stat}$ after the
removal of $\xi_{2}$, but $\xi_{1}$ will be, and by the triangle inequality

$|\xi-\xi_{1}|\leq|\xi-\xi_{2}|+|\xi_{2}-\xi_{1}|\leq n\delta_{\rm
ball}r_{\xi_{2}}+\delta_{\rm ball}r_{\xi_{1}}\leq(n+1)\delta_{\rm
ball}r_{\xi_{1}},$

completing the inductive step. ∎

As a consequence, we obtain the following.

######  Corollary 3.2.

If $J>2$ and $0<\delta_{\rm ball}\leq 1/(2(J-2))$ then, after the removal
algorithm has run, for every stationary point $\xi$ there exists
$\xi^{\prime}\in\mathcal{P}_{\rm stat}$ such that
$\xi\in\Omega_{\xi^{\prime}}$ and $\mathop{\rm
dist}(\xi,\partial\Omega_{\xi^{\prime}})\geq r_{\xi^{\prime}}/2$.

###  3.2 Region of no return for SD contours

The following result establishes a _region of no return_ : once an SD contour
enters this region, we can say with certainty which valley it will converge
to. The idea behind this result is that in the region of no return the highest
degree term $\alpha_{J}z^{J}$ of the polynomial $g$ is sufficiently dominant
over the lower degree terms that the SD contours inside the region converge to
the same valley as those corresponding to the monomial phase
$\alpha_{J}z^{J}$.

######  Theorem 3.3 (Region of no return).

Let $g$, $\mathcal{V}$ and $R_{v}$, for $v\in\mathcal{V}$, be as in (3), (5)
and (11). The regions $R_{v}$, $v\in\mathcal{V}$, contain no stationary points
of $g$. Furthermore, if an SD contour enters $R_{v}$ for some
$v\in\mathcal{V}$, it never leaves $R_{v}$.

###### Proof.

That $R_{v}$ contains no stationary points follows because if $G(r,\theta)>0$
then

$J|\alpha_{j}||z|^{J-1}>\sum_{j=1}^{J-1}j|\alpha_{j}||z|^{j-1}\geq\left|\sum_{j=1}^{J-1}j\alpha_{j}z^{j-1}\right|,$

so that $g^{\prime}(z)\neq 0$.

Now fix $v\in\mathcal{V}$. Given $\theta^{\prime}\in(0,\pi/(2J))$ and $R>0$ we
define the sector

$S_{v}(R,\theta^{\prime}):=\\{z\in\C:|\arg{z}-v|_{2\pi}<\theta^{\prime}\text{
and }|z|>R\\},$

with $|\cdot|_{2\pi}$ defined as in (12). We also define the function

$\displaystyle\tilde{G}(R,\theta^{\prime}):=|J||\alpha_{J}|R^{J-1}\min\left(\sin{J\theta^{\prime}},\cos{J\theta^{\prime}}\right)-\sum_{j=1}^{J-1}j|\alpha_{j}|R^{j-1},$
(24)

which for each fixed $\theta^{\prime}$ is a polynomial in $R$ of degree $J-1$.

We claim that if $\theta^{\prime}\in(0,\pi/(2J))$ and
$\tilde{G}(R,\theta^{\prime})>0$, then if an SD contour enters
$S_{v}(R,\theta^{\prime})$ it never leaves $S_{v}(R,\theta^{\prime})$. To
prove this, we show that if an SD contour intersects $\partial
S_{v}(R,\theta^{\prime})$ then the direction of descent always points into
$S_{v}(R,\theta^{\prime})$. Since $\partial S_{v}(R,\theta^{\prime})$ is the
union of the sets

$\\{z\in\C:|\arg z-v|_{2\pi}\leq\theta^{\prime}\text{ and }|z|=R\\}$

and

$\\{z\in\C:|\arg z-v|_{2\pi}=\theta^{\prime}\text{ and }|z|\in[R,\infty)\\},$

it suffices to show that, in polar coordinates $(r,\theta)$,

$\displaystyle\Im\frac{\partial g}{\partial r}>0,\quad\text{for
}|\theta-v|_{2\pi}\leq\theta^{\prime}\text{ and }r=R,$ (25)
$\displaystyle\mp\Im\frac{1}{r}\frac{\partial
g}{\partial\theta}>0,\quad\text{for }\theta=v\pm\theta^{\prime}\text{ (mod
$2\pi$) and }r\geq R.$ (26)

For (25), let $|\theta-v|_{2\pi}\leq\theta^{\prime}$. Since

$\displaystyle\frac{\partial g(r\mathrm{e}^{\textrm{i}\theta})}{\partial r}$
$\displaystyle=\sum_{j=1}^{J}j\alpha_{j}\mathrm{e}^{\textrm{i}j\theta}r^{j-1}$

and
$\Im[\alpha_{J}\mathrm{e}^{\textrm{i}J\theta}]=|\alpha_{J}|\cos{(J|\theta-v|_{2\pi}})$
(using the definition of $v$) we have that

$\displaystyle\Im\frac{\partial g(r\mathrm{e}^{\textrm{i}\theta})}{\partial
r}$ $\displaystyle\geq
J|\alpha_{J}|r^{J-1}\cos(J|\theta-v|_{2\pi})-\sum_{j=1}^{J-1}j|\alpha_{j}|r^{j-1},$

so a sufficient condition for (25) to hold is that

$\displaystyle
J|\alpha_{J}|R^{J-1}\cos(J\theta^{\prime})-\sum_{j=1}^{J-1}j|\alpha_{j}|R^{j-1}>0.$
(27)

For (26), let $\theta=v\pm\theta^{\prime}$ (mod $2\pi$). Since

$\displaystyle\frac{1}{r}\frac{\partial
g(r\mathrm{e}^{\textrm{i}\theta})}{\partial\theta}$
$\displaystyle=\sum_{j=1}^{J}\textrm{i}j\alpha_{j}\mathrm{e}^{\textrm{i}j\theta}r^{j-1},$

and
$\Im[\textrm{i}\alpha_{J}\mathrm{e}^{\textrm{i}J\theta}]=\Im[\textrm{i}\alpha_{J}\mathrm{e}^{\textrm{i}J(v\pm\theta^{\prime})}]=\mp|\alpha_{J}|\sin(J\theta^{\prime})$
we have that

$\displaystyle\left.\mp\Im\frac{1}{r}\frac{\partial
g(r\mathrm{e}^{\textrm{i}\theta})}{\partial\theta}\right|_{\theta=v\pm\theta^{\prime}}$
$\displaystyle\geq
J|\alpha_{J}|r^{J-1}\sin(J\theta^{\prime})-\sum_{j=1}^{J-1}j|\alpha_{j}|r^{j-1}=:\phi(r).$

The function $\phi(r)$ has the property that if $R>0$ and $\phi(R)>0$ then
$\phi(r)>0$ for all $r\geq R$. To see this, note that

$\phi(r)=r^{J-1}\Big{(}J|\alpha_{J}|\sin(J\theta^{\prime})-\sum_{j=1}^{J-1}j|\alpha_{j}|r^{j-J}\Big{)},$

and that the term in brackets is a strictly decreasing function of $r$, which
tends to $-\infty$ as $r\to 0$ and to $J|\alpha_{J}|\sin(J\theta^{\prime})>0$
as $r\to\infty$. Hence a sufficient condition for (26) is that $\phi(R)>0$,
i.e.

$\displaystyle
J|\alpha_{J}|R^{J-1}\cos(J\theta^{\prime})-\sum_{j=1}^{J-1}j|\alpha_{j}|R^{j-1}>0.$
(28)

Since the assumption $\tilde{G}(R,\theta^{\prime})>0$ implies both (27) and
(28), our claim is proved.

The statement of the theorem then follows by noting that the region $R_{v}$ is
the union of all the sectors $S_{v}(R,\theta^{\prime})$ such that
$0<\theta^{\prime}\leq\pi/(2J)$ and $\tilde{G}(R,\theta^{\prime})>0$. We note
that if $0<\theta^{\prime}<\pi/(4J)$ then
$\sin{J\theta^{\prime}}<\sin{\pi/4}$, so that if
$\tilde{G}(R,\theta^{\prime})>0$ then $\tilde{G}(R,\pi/(4J))>0$. This implies
that the union can actually be taken over
$\pi/(4J)\leq\theta^{\prime}<\pi/(2J)$ only, justifying the definition of the
function $G$ in (13). ∎

###  3.3 Quadrature error

In §2.2 we defined the non-oscillatory region as a union of balls on which the
exponential $\textrm{e}^{\textrm{i}\omega g(z)}$ undergoes a bounded number of
oscillations. Here we show that the definition (8) strikes a balance between
the accuracy of our quadrature approximations to the integrals outside and
inside this region.

####  3.3.1 Quadrature in the non-oscillatory region

The Type 1 straight line contour integrals between points in the non-
oscillatory region are evaluated using Gauss-Legendre quadrature, as detailed
in §2.6.1. To assess the accuracy of this we note the following theorem, which
is a simple consequence of the standard error analysis presented in [15, Chap.
19].

######  Theorem 3.4.

Let $z_{0},z_{1}\in\C$. Suppose that $\gamma$ is a straight-line contour in
starting at $z_{0}$ and ending at $z_{1}$ and that there exists $\rho>0$,
$C>0$ and $\xi_{*}\in\C$ such that $f$ is analytic and bounded in
$z_{[z_{0},z_{1}]}(B_{\rho})$, where $B_{\rho}$ is a standard Bernstein
ellipse (relative to $[-1,1]$) and $z_{[z_{0},z_{1}]}$ is defined as in (15),
and

$\displaystyle\omega|g(\xi_{*})-g(z)|\leq C,\qquad z\in
z_{[z_{0},z_{1}]}(B_{\rho}).$ (29)

Let $I$ and $Q$ denote the left- and right-hand sides of (16), respectively.
Then, for some $\tilde{C}>0$, depending only on $\rho$,

$\displaystyle\left|I-Q\right|\leq\tilde{C}|z_{1}-z_{0}|\|f\|_{L^{\infty}(z_{[z_{0},z_{1}]}(B_{\rho}))}\textrm{e}^{-\omega\Im[g(\xi_{*})]}\textrm{e}^{C}\rho^{-2N}.$
(30)

###### Proof.

Noting that

$I=\textrm{e}^{\textrm{i}\omega
g(\xi_{*})}\int_{\gamma}f(z)\textrm{e}^{\textrm{i}\omega(g(z)-g(\xi_{*}))}\,\operatorname{d}\\!z$

and that

$|f(z)\textrm{e}^{\textrm{i}\omega(g(z)-g(\xi_{*}))}|\leq\|f\|_{L^{\infty}(z_{[z_{0},z_{1}]}(B_{\rho}))}\textrm{e}^{C},\qquad
z\in z_{[z_{0},z_{1}]}(B_{\rho}),$

the result follows from [15, Thm 19.3]. ∎

Theorem 3.4 motivates the definition of the non-oscillatory region in (8).
Indeed, if the assumptions of Theorem 3.4 hold with $\rho$ and $C$ independent
of $\omega$ then the bound (30) guarantees $\omega$-independent exponential
convergence for $\omega$ bounded away from zero. However, even when (8) is
satisfied, the relationship between $\xi_{*}$, $\rho$, $C$ and $\omega$ is
beyond our control in general because the ellipse may extend beyond the non-
oscillatory region, so that $C>C_{\rm ball}$. Thus we cannot control the
factor $\textrm{e}^{C}$ entirely based on condition (8).

Still, the bound (30) shows that the quadrature error decreases with
increasing $N$. The precise rate of decrease depends on a balance between the
decay of $\rho^{-2N}$ and the growth of $\textrm{e}^{C}$ and
$\|f\|_{L^{\infty}(z_{[z_{0},z_{1}]}(B_{\rho}))}$ for increasing $\rho$. We
quantify this in the special case of monomial phase in §3.3.3.

####  3.3.2 Quadrature for the SD contours

For Type 2 or Type 3 integrals along SD contours we use either Gauss-Laguerre
or (possibly truncated) Gauss-Legendre quadrature, as detailed in §2.6.2 and
§2.6.3. We expect these rules to converge rapidly to the true value of the
integral as the number of quadrature points $N$ tends to infinity, provided
that the integrand is analytic and bounded in a suitable region of the complex
$\tilde{p}$ plane.

For Gauss-Laguerre the following result appeared recently in [17, Thm 6.3].

######  Theorem 3.5.

Suppose that $\tilde{f}$ is analytic inside and on the parabola
$P_{\rho}:=\\{z\in\C:\sqrt{-z}=\rho\\}$ for some $\rho>0$, where the branch
cut is along the positive real axis and $\sqrt{-z}$ is real and positive on
the negative real axis, that $\tilde{f}$ grows at most algebraically as
$z\to\infty$ inside the parabola, and that the integral

$\mathcal{K}_{\rho}:=\int_{P_{\rho}}|\textrm{e}^{-z}\sqrt{-z}\tilde{f}(z)|\,\operatorname{d}\\!z$

is finite. Let $I$ and $Q$ denote the left- and right-hand sides of (18),
respectively. Then

$\displaystyle\left|I-Q\right|\leq\mathcal{K}_{\rho}\frac{\textrm{e}^{-\omega\Im[g(\eta)]}}{\omega}\textrm{e}^{-4\rho\sqrt{N}}.$
(31)

This result implies that our Gauss-Laguerre quadrature approximation should
converge root-exponentially as $N\to\infty$, provided that $f$ is sufficiently
well-behaved at infinity. The presence of singularities in the complex
$\tilde{p}$-plane limits the size of $\rho$, and hence the convergence rate.
We know from (17) that our integrand is singular at points $\tilde{p}\in\C$
where $g^{\prime}(h_{\eta}(\tilde{p}/\omega))=0$, i.e. where
$h_{\eta}(\tilde{p}/\omega)=\xi$ for some stationary point $\xi$. Since we
only trace SD contours outside the non-oscillatory region (which contains the
stationary points), we know that there cannot be singularities on the SD
contour itself. If the start point $\eta$ lies on an SD contour emanating from
a stationary point $\xi$ then we expect there to be a singularity in the
$\tilde{p}$-plane at $\tilde{p}=\omega\Im[g(\xi)-g(\eta)]<0$. We show in
§3.3.3 that in the special case of monomial phase this singularity lies at
$\tilde{p}=-C_{\rm ball}$, which implies root-exponential convergence
independent of $\omega$ for $\omega$ bounded away from zero. Determining the
locations of the other possible singularities in the complex $\tilde{p}$-plane
is more challenging, since it involves study of the (multivalued) inverse of
$g$. We leave further theoretical investigation of this to future work.

####  3.3.3 Results for monomial phase

It is instructive to consider the special case of a monomial phase
$g(z)=z^{J}$ for some $J\in\N$. In this case there is a single stationary
point of order $J-1$ at $\xi=0$, and $g(0)=0$. Following the prescription (8),
we obtain a ball radius

$r_{0}=(C_{\rm ball}/\omega)^{1/J}.$

We first consider a Type 1 integral in the non-oscillatory region. For
simplicity we choose $f(z)\equiv 1$. Specifically, we consider the evaluation
of the integral

$\int_{0}^{r_{0}\textrm{e}^{\textrm{i}\theta}}\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z,$

for some $\theta\in[0,2\pi]$. Taking $\xi_{*}=0$, we can apply Theorem 3.4
with any $\rho>1$, and the resulting scaled and translated Bernstein ellipse
surrounding $[0,r_{0}\textrm{e}^{\textrm{i}\theta}]$ is tightly contained in
the disc $|z|\leq sr_{0}$, where $\rho$ and $s$ are related by

$\rho=2s-1+\sqrt{(2s-1)^{2}-1}=2s-1+2\sqrt{s^{2}-s}.$

Hence condition (29) is satisfied, independently of $\theta$, with

$C=C_{\rm ball}s^{J},$

which is independent of $\omega$ but dependent on $J$. When $s$ is large, we
have $\rho\approx 4s$, and in this regime the error bound provided by (30) for
Gauss-Legendre quadrature is approximately proportional to

$(C_{\rm ball}/\omega)^{1/J}\textrm{e}^{C_{\rm ball}s^{J}}(4s)^{-2N}.$

As a function of $s$, with $J$ and $N$ fixed, this quantity is minimised where
its $s$-derivative vanishes, which occurs where

$C_{\rm ball}Js^{J}-2N=0,$

i.e. where

$s=\left(\frac{2N}{C_{\rm ball}J}\right)^{1/J}.$

Accordingly, the error bound is approximately proportional to

$(C_{\rm ball}/\omega)^{1/J}16^{J}\left(\frac{8eN}{C_{\rm
ball}J}\right)^{-2N/J}.$

Thus we expect super-exponential convergence as $N\to\infty$ for fixed $J$.
However, we expect the convergence to be slower the larger $J$ is.

Next we consider a Type 2 integral over an SD contour, again with $f(z)\equiv
1$. Specifically, we consider the evaluation of the integral

$\int_{r_{0}\textrm{e}^{\textrm{i}v}}^{\infty\textrm{e}^{\textrm{i}v}}\textrm{e}^{\textrm{i}\omega
g(z)}\,\operatorname{d}\\!z,$

where $v=((2j+1/2)\pi)/J$ for some $j\in\\{1,\ldots,J\\}$. Following our
method, the contour is parametrized by

$h_{\eta}(p)=(r_{0}^{J}+p)^{1/J}\textrm{e}^{\textrm{i}v},\qquad
p\in[0,\infty),$

and, recalling (10) and (17), after rescaling $p=\tilde{p}/\omega$ the
integral becomes

$\frac{\textrm{e}^{-C_{\rm
ball}}\textrm{e}^{\textrm{i}v}}{\omega^{1/J}J}\int_{0}^{\infty}(C_{\rm
ball}+\tilde{p})^{1/J-1}\textrm{e}^{-\tilde{p}}\,\operatorname{d}\\!\tilde{p}.$

The integrand has a branch point at

$\tilde{p}=-C_{\rm ball},$

but we note that the distance between the branch point and the positive real
$\tilde{p}$-axis equals $C_{\rm ball}$, which is independent of both $\omega$
and $J$.

For truncated Gauss-Legendre the relevant theory can be found in [15, Chap.
19] (and see also [16]). Due to the branch point at $\tilde{p}=-C_{\rm ball}$,
as $N\to\infty$ we obtain exponential convergence to the integral over the
interval $[0,P]$, where $P$ is given by either (21) or (23). In the case where
$P=L$, by the definition of $L$ in (22), we expect the truncation error to
have relative order $\delta_{\rm quad}$.

####  3.3.4 Number and distribution of quadrature points

PathFinder uses a fixed number $N$ of quadrature points on each contributing
contour, and that number is the same both for integrals within and outside the
non-oscillatory region, i.e., for Gauss–Legendre and Gauss–Laguerre
quadrature. Thus, increasing the single parameter $N$ provides a way of
uniformly improving accuracy.

The theoretical results in this section (specifically, Theorems 3.4 and 3.5)
imply that the precise rate of improvement with respect to $N$ depends on the
type of integral being approximated. They suggest even that a different
strategy for the distribution of quadrature points may be superior. Indeed,
exponential convergence of Gauss–Legendre for Type 1 integrals in the non-
oscillatory region is not balanced with root-exponential convergence of
Gauss–Laguerre for Type 2 integrals outside. Similarly, convergence rates of
Gauss-Laguerre and truncated Gauss-Legendre outside the non-oscillatory region
are different. Our choice of a fixed parameter $N$ is inspired on the one hand
by simplicity, and on the other hand by the lack of robust methods to optimize
parameters in alternative schemes. For example, we have shown in §3.3.3 that
the convergence rate of Gauss-Legendre for Type 1 integrals may depend on the
order of nearby stationary points. While this can be quantified precisely for
the case of monomial phase, it is not at all clear how to generalise this
analysis when a cluster of multiple stationary points is present. Hence,
stationary point order is a quantity that we deliberately do not explicitly
compute, estimate or rely on in any way. Implicitly, of course, it plays a big
role, and it does so mainly via the definition of the ball of the radius in
(8).

The main practical benefit of the theoretical analysis of quadrature error in
this section is the guarantee that $N$ is a robust parameter for improving
accuracy. Concerning possible future improvements, rather than attempting to
optimize the quadrature point distribution a priori, we believe a more
promising development would be the ability to invoke standard adaptive
quadrature schemes along the contours for a given function $f$. However, it
should be borne in mind that quadrature forms just one step in our algorithm,
and that the other steps (particularly the SD path tracing) incur a non-
negligible cost overhead, that should also be considered when trying to
further optimize performance.

##  4 Further implementation aspects

In this section we discuss some additional aspects of the implementation of
our algorithm in PathFinder.

###  4.1 Default parameter values

In Table 1 we list the user-specified parameters in our algorithm, along with
the default values used in all our numerical results in §5. These were
determined as the result of extensive numerical experiments on a range of
examples, not detailed here. Instructions on how to adjust these parameters
away from their default values can be found at
github.com/AndrewGibbs/PathFinder.

Parameter | Domain | Meaning | Default
---|---|---|---
$C_{\rm ball}$ | $(0,\infty)$ | Governs maximum number of | $2\pi$
|  | oscillations across each non-oscillatory |
|  | ball (and hence the ball radius) |
$N_{\rm ball}$ |  | Number of rays used | 16
|  | when determining the ball radius |
$\delta_{\rm ball}$ | $(0,1)$ | Governs when overlapping balls | $10^{-3}/(2\max(J-2,1))$
|  | should be amalgamated |
$\delta_{\rm ODE}$ | $(0,1)$ | Governs the local step size in | $0.1$
|  | the ODE solver for SD path tracing |
$\delta_{\rm coarse}$ | $(0,1)$ | Tolerance for the increment in the | $10^{-2}$
|  | Newton iteration in the SD path tracing |
$\delta_{\rm fine}$ | $(0,1)$ | Tolerance for the increment in the | $10^{-13}$
($<\delta_{\rm coarse}$) |  | Newton iteration in the quadrature |
$\delta_{\rm quad}$ | $(0,1)$ | Governs when the contribution | $10^{-16}$
|  | from an integral on the quasi-SD |
|  | deformation is computed |
$N$ |  | Number of quadrature points to use | no default
|  | in each integral evaluated in step 76 |
Table 1: User-specified parameters and their default values in PathFinder.

###  4.2 Small $\omega$

While our algorithm is geared towards the case where $\omega$ is moderate or
large, we make a brief comment on the case where $\omega$ is small. If
$\Gamma$ is infinite then the integral (1) typically diverges for $\omega=0$.
However, if $\Gamma$ is finite then the integral converges for $\omega=0$ and
for small enough $\omega$ it is non-oscillatory. In PathFinder we detect and
deal with this case in the following way. If both endpoints are finite, then
before starting step 1 of the algorithm we construct non-oscillatory balls
around the endpoints (using the process in §2.2) and check whether the balls
intersect non-trivially. If so, we apply standard Gauss-Legendre quadrature to
evaluate (1); if not, the balls are discarded and we proceed with the rest of
the algorithm.

###  4.3 The case $J=1$

In the case $J=1$ (linear phase) there are no stationary points, and our
algorithm simplifies dramatically. Furthermore, the SD contours are simply
parallel straight lines in the direction of the single valley at angle
$\pi/2-\arg(\alpha_{1})$, and there is no need to trace them numerically.
Hence when $J=1$ PathFinder skips the ODE contour tracing step and exploits
the exact characterization of the SD contours mentioned above.

###  4.4 Specifying infinite endpoints

In the description of our algorithm in §2 we made the assumption that any
infinite endpoint of the contour $\Gamma$ should be at a valley
$v\in\mathcal{V}$. PathFinder is actually more flexible than this. The user is
permitted to specify an infinite endpoint at any
$\theta\in[v-\pi/(2J),v+\pi/(2J)]$ and the code will automatically adjust this
to equal $v$. The case $\theta=v\pm\pi/(2J)$ is delicate because the highest
order term in the phase does not provide exponential decay along the contour.
Nonetheless, we include it, because in applications one often encounters this
case, with the integral converging conditionally (under appropriate
assumptions on $f$) and the contour deformation to $v$ being justified by
Jordan’s Lemma.

##  5 Numerical results

In this section we present numerical results illustrating the performance of
our algorithm and its implementation in PathFinder. All results in this
section were produced using PathFinder Version 1.0 [8].

###  5.1 A “generic” example

(a) $\omega=0.01$

(b) $\omega=1$

(c) $\omega=5$

(d) $\omega=50$

Figure 3: PathFinder output (cf. Figure 1) with $N=10$ for the approximation
of (32).

We begin by illustrating the performance of PathFinder on the integral

$\displaystyle
I=\int_{-1}^{1}(2z^{4}+7z^{3}+z^{2}+8z+2)\textrm{e}^{\textrm{i}\omega(3z^{9}+z^{8}+4z^{7}+z^{6}+5z^{5}+9z^{4}+2z^{3}+6z^{2}+5z+3)}\,\operatorname{d}\\!z,$
(32)

where, to convey the message that our approach is applicable to truly
“generic” amplitudes and polynomial phase functions, the coefficients of $f$
and $g$ are chosen to be the first 5 digits of e and the first 10 digits of
$\pi$, respectively. This can be approximated by PathFinder via the Matlab
code (cf. (2))

PathFinder(-1,1,@(z) 2*z.^4+7*z.^3+z.^2+8*z+2,... [3 1 4 1 5 9 2 6 5
3],omega,N)

In Figure 3 we plot the quasi-SD deformations and quadrature point
distributions (using the PathFinder 'plot' option) for (32) for
$\omega\in\\{0.01,1,5,50\\}$ and $N=10$. As explained in §2.2, for smaller
$\omega$ the non-oscillatory balls are larger, and can overlap, while for
larger $\omega$ they shrink around the stationary points. In more detail, in
Figure 3(a) ($\omega=0.01$), $\omega$ is small enough that both endpoints are
inside the same non-oscillatory ball. Hence the integral is treated as non-
oscillatory and is approximated by Gauss-Legendre quadrature along a single
straight-line contour. In Figure 3(b) ($\omega=1$), $\omega$ is still small
enough that many of the balls overlap, and the quasi-SD deformation comprises
two SD contours (one from an exit and one from an endpoint) plus four
straight-line contours in the non-oscillatory region. In Figure 3(c)
($\omega=5$), $\omega$ is large enough that only two balls overlap, and the
quasi-SD deformation comprises five SD contours (two from endpoints, two from
exits to valleys, and one from an exit to an entrance), plus four straight-
line contours in the non-oscillatory region. Finally, in Figure 3(d)
($\omega=50$), $\omega$ is so large that none of the balls overlap, and the
quasi-SD deformation comprises eight contributing SD contours (two from
endpoints and six from exits to valleys), plus three straight-line contours in
the non-oscillatory region. However, in this case the two SD contours and one
straight-line contour associated with the stationary point near
$0.2+0.5\textrm{i}$ are judged to make a negligible contribution to the
integral, so are not assigned any quadrature points. We emphasize that this
intricate behaviour is fully automated, with no expert input required from the
user.

(a) Accuracy

(b) Timings

Figure 4: Accuracy (a) and timings (b) of the PathFinder approximation of
(32).

In Figure 4(a) we plot the error in the PathFinder approximation of (32),
compared to reference values computed using the Julia QuadGK package when
$\omega<500$, and using PathFinder with $N=500$ when $\omega\geq 500$. For
fixed $\omega$ we observe rapid convergence as $N\to\infty$, at a rate that
appears independent of $\omega$. In Figure 4(b) we show the associated
computation times, which remain bounded as $\omega$ increases.

###  5.2 Coalescence and the Airy function

The canonical example of an integral with two coalescing stationary points is
provided by the integral representation for the Airy function, viz. (see [1,
9.5.4])

$\displaystyle\mathrm{Ai}(x)=\int_{\infty\textrm{e}^{-\textrm{i}\pi/3}}^{\infty\textrm{e}^{\textrm{i}\pi/3}}\textrm{e}^{z^{3}/3-xz}\,\operatorname{d}\\!z=\int_{\infty\textrm{e}^{-\textrm{i}\pi/3}}^{\infty\textrm{e}^{\textrm{i}\pi/3}}\textrm{e}^{\textrm{i}(-\textrm{i}(z^{3}/3-xz))}\,\operatorname{d}\\!z,\qquad
x\in\C,$ (33)

which is of the form (1) with $f\equiv 1$, $\omega=1$ and
$g(z;x)=-\textrm{i}(z^{3}/3-xz)$. Up to a change of variable this is the same
example for which, as mentioned in §1, a bespoke, complex Gaussian quadrature
rule was developed in [10]. $\mathrm{Ai}$ can be approximated by PathFinder
via the Matlab code

Ai = @(x) PathFinder(-pi/3,pi/3,[],... -1i*[1/3 0 -x 0],1,N,'infcontour',[true
true])

where the input [] for f indicates that $f\equiv 1$.

(a) $x=-5$

(b) $x=-1$

(c) $x=-0.5$

(d) $x=0$

(e) $x=5$

Figure 5: PathFinder output with $N=20$ for the approximation of
$\mathrm{Ai}(x)$ via (33) at various $x$, showing the stationary point
coalescence at $x=0$.

In Figure 5 we plot the quasi-SD deformations, along with the distribution of
quadrature points for $N=20$, for the evaluation of $\mathrm{Ai}(x)$ at
$x\in\\{-5,-1,-0.5,0,5\\}$. Here one observes in detail how our algorithm
deals with stationary point coalescence, as the non-oscillatory balls overlap,
merge, then split. In Figure 5(a) the quasi-SD deformation comprises four SD
contours from exits, plus two straight-line contours inside balls (which do
not go via stationary points). In Figure 5(b) the balls overlap and this
changes to two SD contours from exits plus three straight-line contours inside
balls (which go via both stationary points). In Figure 5(c) the balls overlap
enough that both stationary points are contained in both balls, so we get two
SD contours from exits plus just two straight-line contours inside balls
(which go via only one of the stationary points). In Figure 5(d) the balls
have merged completely and in addition to the two SD contours from exits there
is just one straight-line contour inside a ball (which does not go via the
stationary point). In Figure 5(e) the balls have split again, but we see the
same deformation structure as in Figure 5(d). Again, we emphasize that these
calculations are fully automated.

Figure 6: Accuracy of PathFinder approximation of $\mathrm{Ai}(x)$ for
different $N$.

In Figure 6 we show the accuracy of the PathFinder approximation for this
example as a function of $x\in[-10,4]$, for different $N$. Our reference is
the built-in Matlab command airy. We note that between $x=-3$ and $x=0$ the
error for the smaller values of $N$ undergoes some jumps. These are due to the
fact that near stationary point coalescence the topology of the quasi-SD
deformation, the number of contours constituting it, and hence the total
number of quadrature points along it (recall (14)), all change discontinuously
as a function of $x$ (as illustrated in Figure 5). However, as $N$ increases
we see a clear, approximately exponential decrease in the error, and, although
the rate of decrease depends slightly on $x$ (because of the factors mentioned
above), for $N=30$ we achieve approximately $10^{-13}$ error uniformly across
the interval.

###  5.3 A high order stationary point - comparison with Mathematica’s
implementation of Levin quadrature

Figure 7: PathFinder output for (34) with $\omega=100,000$ and $N=50$.

We now consider the integral

$I=\int_{-1}^{1}\sin(z)\textrm{e}^{\textrm{i}\omega
z^{9}}\operatorname{d}\\!z,$ (34)

which has a stationary point of order 8 at the origin. The integral (34) can
be approximated by PathFinder via the command

PathFinder(-1,1,@(z) sin(z),[1 0 0 0 0 0 0 0 0 0],omega,N)

Figure 7 shows the quasi-SD deformation and quadrature point distribution
obtained by PathFinder for $\omega=100,000$ and $N=50$. There are small
contributions from the endpoints, but the main contribution comes from the
ball containing the stationary point.

In the Mathematica documentation [18, pp75-86], it is stated that oscillatory
integrals with monomial phase functions such as (34) can be evaluated
efficiently using the built-in Mathematica function NIntegrate, via its
implementation of Levin quadrature (which is described, e.g., in [7, §3.3]).
To do this one can use the Mathematica command:

NIntegrate[Sin[x]Exp[omega*I*x^9],{x,-1,1},
Method->{"LevinRule","Kernel"->Exp[omega*I*x^9]}]

(a) Accuracy (compared to Mathematica)

(b) Timings

Figure 8: Accuracy (a) and timings (b) of the PathFinder approximation of
(34), compared to Mathematica’s NIntegrate command.

In Figure 8(a) we show a plot of the relative accuracy of our PathFinder
approximation, compared to the Mathematica approximation (using the default
settings), as a function of $\omega$, for different $N$ values. For all three
$N$ values the accuracy of our approximation is approximately uniform in
$\omega$, and for $N=50$ our approximation agrees with Mathematica’s to approx
$13$ digits. In Figure 8(b) we report the corresponding computation times
(averaged over 100 identical runs) for the Mathematica routine and for the
PathFinder approximation with $N=50$. These results were obtained on a laptop
(i7-1185G7, 32GB RAM) running Mathematica v13.0 and Matlab v2021b. The results
suggest that PathFinder is highly competitive with Mathematica for this
problem, especially for large $\omega$.

###  5.4 Coalescence to a high order stationary point

(a) $r=0.001$

(b) $r=0.35$

(c) $r=0.5$

(d) $r=1$

Figure 9: PathFinder output for the approximation of (35) with $\omega=1000$
and $N=15$.

We now investigate the robustness of our algorithm in the presence of a large
number of coalescing stationary points. Specifically, we consider the integral

$\displaystyle\int_{-1}^{1}\textrm{e}^{\textrm{i}\omega(z^{7}/7-r^{6}z)}\,\operatorname{d}\\!z,$
(35)

where $r\geq 0$ is a parameter controlling the coalescence. For $r>0$ there
are 6 stationary points with $|\xi|=r$, namely the solutions of
$\xi^{6}=r^{6}$, and for $r=0$ there is a single stationary point of order
$6$. To evaluate this integral in PathFinder for a given $r$, one can use the
command

PathFinder(-1,1,[],[1/7 0 0 0 0 0 -r^6 0],omega,N)

In Figure 9 we plot the quasi-SD deformations and quadrature point
distributions for a some different $r$ values, showing how the balls first
intersect and then merge as $r\to 0$.

(a) Convergence

(b) Timings

Figure 10: Accuracy (a) and timings (b) for (35) for $p=6$ and $r=0.01$.

(a) $N=10$

(b) $N=50$

Figure 11: Accuracy for (35) as $r\to 0$ for $p=6$, for $N=10$ (a) and $N=50$
(b).

In Figure 10 we show convergence (with respect to a PathFinder reference with
$N=500$) and CPU times (averaged over 100 runs) for fixed $r=0.01$. We see
that both the error and the CPU time are essentially independent of $\omega$
in this case. In Figure 11 we plot errors for two fixed $N$ values $N=10,50$,
as a function of $r$. We observe that as $r\to 0$, the error stays bounded.
For $N=10$ the error jumps up between $r=10^{-3}$ and $r=10^{-2}$, at a point
depending on $\omega$. This represents the point at which the balls around the
stationary points merge, resulting in a reduction of $N_{\rm tot}$, and hence
a reduction in accuracy. But after this point we observe no further reduction
in accuracy as $r\to 0$. We remark that for sufficiently small $r>0$ the six
stationary points are numerically indistinguishable, but this isn’t a problem
for our algorithm because in that case the problem will be treated identically
to that of a monomial phase.

###  5.5 Canonical cuspoid integrals and their generalisations

In this section we show how our algorithm can be applied to the computation of
some of the canonical integrals catalogued by Berry and Howls in [1, §36],
which, as mentioned already in §1, are of fundamental importance in numerous
application areas including optics, acoustics and quantum mechanics.

In this context, our algorithm is related to that of [11], where an adaptive
contour deformation approach was applied to evaluate the cuspoid integrals
considered in §5.5.1. The algorithm in [11] is similar in spirit to our
approach, in that it deforms the integration contour so that it terminates in
valleys at infinity, and splits the contour into portions close to stationary
points plus portions away from stationary points. However, in contrast to our
approach, the algorithm in [11] does not attempt to trace SD contours, and
hence is susceptible to rounding errors associated with the “violent”
behaviour of the exponential factor $\textrm{e}^{\textrm{i}\omega g(z)}$ when
one is not on a true SD contour - see [11, §2]. Furthermore, while the
algorithm in [11] was specialised to the case of integration over the real
line, our algorithm can handle much more general contours, as we illustrate in
§5.5.2.

####  5.5.1 Cuspoid integrals

The so-called “cuspoid integrals” listed in [1, §36.2.4] are all of the form
(1) with polynomial phase $g$ and unit amplitude $f\equiv 1$, unit frequency
$\omega=1$, and integration along the real line. Our algorithm is ideally
suited to the evaluation of these integrals, and to demonstrate this we
compute two of them. In the notation of [1, §36], we consider the cusp
catastrophe integral

$\displaystyle\Psi_{2}(x,y)=P(y,x)=\int_{-\infty}^{\infty}\textrm{e}^{\textrm{i}(t^{4}+yt^{2}+xt)}\,\operatorname{d}\\!t,$
(36)

where $P$ is the Pearcey function, and the swallowtail catastrophe integral

$\displaystyle\Psi_{3}(x,y,z)=\int_{-\infty}^{\infty}\textrm{e}^{\textrm{i}(t^{5}+zt^{3}+yt^{2}+xt)}\,\operatorname{d}\\!t.$
(37)

Both exhibit coalescence of stationary points on certain algebraic varieties
(see [1, §36.5(ii)]) on which both the first and second derivatives of the
phase function vanish. In the case of (36) this occurs when

$\displaystyle y=-\frac{3}{2}|x|^{2/3},$ (38)

and for (37) this occurs when

$\displaystyle
400x^{3}-360x^{2}z^{2}-135y^{4}-27y^{2}z^{3}+540xy^{2}z+81xz^{4}=0.$ (39)

(a) $|\Psi_{2}(x,y)|$

(b) $|\Psi_{3}(x,y,-7.5)|$

Figure 12: Magnitude plots of (36) and (37), with coalescence curves (38) and
(39) (the latter with $z=-7.5$) superimposed in black. The computational grid
was of size $100\times 100$ for (a) and $500\times 500$ for (b).

The integrals (36) and (37) can be computed in Pathfinder via the commands

Psi2 = @(x,y) PathFinder(pi,0,[],[1 0 y x 0],... 1,N,'infcontour',[true true])
Psi3 = @(x,y,z) PathFinder(pi,0,[],[1 0 z y x 0],... 1,N,'infcontour',[true
true])

Figure 12 shows plots of the magnitude of (36) and (37) (the latter over the
plane $z=-7.5$), computed using PathFinder with the default settings and
$N=50$. The plots agree qualitatively with those presented in [1, Figs 36.3.1
& 36.6.5], and, for (36), agree quantitatively (to all five decimal places
presented) with the values presented in [11, Table 1]. Computation times on a
small desktop computer (Intel i7-4790, 32GB RAM) were less than a minute for
the cusp (which required the computation of 10000 instances of (36), averaging
$0.005$s per instance) and less than an hour for the swallowtail (which
required 250000 instances of (37), averaging $0.01$s per instance).

####  5.5.2 Generalisations

In [9] the authors considered a family of generalisations of certain canonical
cuspoid integrals, with integration no longer over the real line, but rather
over a complex contour starting and ending at valleys at infinity, and
possibly with a non-unit amplitude function.

A specific aim of [9] was to investigate the relevance of such integrals to
the study of the so-called “inflection point problem”, a canonical problem in
wave scattering originally introduced over 50 years ago by Popov in [12]. This
problem, which remains unsolved in closed form, concerns two-dimensional time-
harmonic wave propagation near a boundary with an inflection point, and seeks
a solution for the wave field near the inflection point that describes the
transition from an incoming “whispering gallery wave” supported on the concave
portion of the boundary, to outgoing “creeping waves” along the convex portion
of the boundary, along with a scattered “searchlight” beam (for details and
further references see [14]).

In this context, in [9, §3.3] the authors studied the family of integrals

$\displaystyle
A_{ij}(x,y)=\int_{\Gamma_{ij}}f(t)\textrm{e}^{\textrm{i}(2t^{5}/5-xt^{4}/2-yt^{2})}\,\operatorname{d}\\!t,$
(40)

where $f(t)$ is some amplitude to be specified, and $\Gamma_{ij}$ denotes any
contour from valley $v_{i}$ to valley $v_{j}$, where
$v_{j}:=(2(j-1)+1/2)\pi/5$, $j=1,\ldots,5$. These integrals have stationary
point coalescence on the cubic curve $y+4x^{3}/27=0$, which suggests that, by
appropriately choosing $f$ and $\Gamma_{ij}$, they might exhibit certain
features of the solution of the inflection point problem. Indeed, in [9, §4]
it was shown that as $x\to-\infty$ near the cubic curve, the integral $A_{32}$
has the character of an incoming whispering gallery type wave, and that, as
$x\to+\infty$ near the cubic curve, the integral $A_{52}$ has the character of
an outgoing creeping wave. However, plots of the resulting fields could not be
presented in [9] due to the lack of a suitable numerical evaluation method and
implementation.

(a) $|A_{32}(x,y)|$

(b) $|A_{52}(x,y)|$

(c) $\Re[A_{32}(k^{1/5}x_{0},k^{3/5}y_{0})\textrm{e}^{\textrm{i}kx_{0}}]$

(d) $\Re[A_{52}(k^{1/5}x_{0},k^{3/5}y_{0})\textrm{e}^{\textrm{i}kx_{0}}]$

Figure 13: Plots of (40) with $f\equiv 1$, along with the associated
approximate Helmholtz equation solutions for $k=40$.

Using PathFinder we are able to remedy this. In Figures 13(a) and 13(b) we
provide plots of the magnitude of $A_{32}$ and $A_{52}$ with $f\equiv 1$. To
evaluate the integrals we used the PathFinder code

A32 = @(x,y) PathFinder(9*pi/10,pi/2,[],[2/5 -x/2 0 -y 0 0],...
1,N,'infcontour',[true true]) A52 = @(x,y) PathFinder(17*pi/10,pi/2,[],[2/5
-x/2 0 -y 0 0],... 1,N,'infcontour',[true true])

We only plot $A_{52}$ above the cubic curve $y+4x^{3}/27=0$, because below
this curve $A_{52}$ becomes exponentially large (cf. [9, Fig. 12(i)]). In
Figures 13(c) and 13(d) we present corresponding plots of the modulated plane
wave

$u(x_{0},y_{0})=A_{ij}(x,y)\textrm{e}^{\textrm{i}kx_{0}},$

where $(x_{0},y_{0})$ are outer variables, related to the inner variables
$(x,y)$ by $x=k^{1/5}x_{0}$, $y=k^{3/5}y_{0}$, which is an asymptotic solution
of the Helmholtz equation $\Delta u+k^{2}u=0$ as $k\to\infty$ in the region
$x_{0}=O(k^{-1/5})$, $y_{0}=O(k^{-3/5})$ [9, §1]. Here one observes the
predicted incoming whispering gallery type behaviour of $A_{32}$ near the top
of Figure 13(c) between $x_{0}=-2$ and $x_{0}=-1$, with oscillations giving
way to an exponentially small field in the caustic shadow, and the predicted
creeping wave type behaviour of $A_{52}$ near the bottom of Figure 13(d)
between $x_{0}=1$ and $x_{0}=2$, with waves propagating along the cubic curve,
shedding rays tangentially.

In ongoing and future studies we plan to use PathFinder to further investigate
the properties of integrals of the form (40), and generalisations involving
different choices of $f$ and higher degree phase functions (see [9]), which we
hope may shed new light on the inflection point problem and related problems
in high frequency wave propagation.

##  6 Acknowledgements

DPH and AG were supported by EPSRC grants EP/S01375X/1 and EP/V053868/1. DPH
and AG would like to thank the Isaac Newton Institute for Mathematical
Sciences for support and hospitality during the programme Mathematical Theory
and Applications of Multiple Wave Scattering when work on this paper was
undertaken. This work was supported by EPSRC grant number EP/R014604/1.
Additionally, AG acknowledges support from KU Leuven project C14/15/05, and DH
acknowledges support FWO-Flanders project G.088.622N.

## References

  * [1] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, release 1.1.3 of 2021-09-15.
  * [2] A. Asheim, A. Deaño, D. Huybrechs, and H. Wang, A Gaussian quadrature rule for oscillatory integrals on a bounded interval, Discr. Cont. Dyn. Sys., 34 (2013), pp. 883–901.
  * [3] N. Bleistein and R. A. Handelsman, Asymptotic Expansions of Integrals, Dover, 1986.  2nd Edition.
  * [4] W. M. Boothby, The topology of the level curves of harmonic functions with critical points, Amer. J. Math., 73 (1951), pp. 512–538.
  * [5] J. P. Boyd, Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding, J. Eng. Math., 56 (2006), pp. 203–219.
  * [6] A. F. Celsus, A. Deaño, D. Huybrechs, and A. Iserles, The kissing polynomials and their Hankel determinants, T. Math. Appl., 6 (2021), p. tnab005.
  * [7] A. Deaño, D. Huybrechs, and A. Iserles, Computing Highly Oscillatory Integrals, vol. 155, SIAM, 2018.
  * [8] A. Gibbs, PathFinder v1.0, https://github.com/AndrewGibbs/PathFinder
/releases/tag/v1.0, (2023).

  * [9] D. P. Hewett, J. R. Ockendon, and V. P. Smyshlyaev, Contour integral solutions of the parabolic wave equation, Wave Motion, 84 (2019), pp. 90–109.
  * [10] D. Huybrechs, A. B. J. Kuijlaars, and N. Lejon, A numerical method for oscillatory integrals with coalescing saddle points, SIAM J. Numer. Anal., 57 (2019), pp. 2707–2729.
  * [11] N. Kirk, J. Connor, and C. Hobbs, An adaptive contour code for the numerical evaluation of the oscillatory cuspoid canonical integrals and their derivatives, Comp. Phys. Comm., 132 (2000), pp. 142 – 165.
  * [12] M. M. Popov, The problem of whispering gallery waves in a neighbourhood of a simple zero of the effective curvature of the boundary, J. Sov. Math. (now J. Math. Sci.), 11 (1979), pp. 791–797.
  * [13] K. H. Rosen, Discrete Mathematics and its Applications, McGraw-Hill, 2019.
  * [14] V. P. Smyshlyaev and I. V. Kamotski, Searchlight asymptotics for high-frequency scattering by boundary inflection, St. Petersburg Math. J., 33 (2022), pp. 387–403.
  * [15] L. N. Trefethen, Approximation theory and approximation practice, SIAM, 2013.
  * [16] , Exactness of quadrature formulas, SIAM Review, 64 (2022), pp. 132–150.
  * [17] H. Wang, Optimal convergence analysis of Laguerre spectral approximations for analytic functions , arXiv preprint 2304.05744, (2023).
  * [18] Wolfram Research, Inc., Advanced Numerical Integration In Mathematica, 2008.  Downloaded 6 April 2023 from library.wolfram.com.
  * [19] R. Wong, Asymptotic Approximations of Integrals, SIAM, 2001.