User Guide
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 5
Abstract
This document is intended to serve as a user guide for the Numerical Steepest Descent (NSD)
Matlab package created by Daan Huybrechs and Andrew Gibbs at KU Leuven.
Contents
1 Setup 1
2 Overview 1
3 Additional options and features 2
3.1 Polynomialphase.................................... 2
3.2 Infinite and semi-infinite contour integrals . . . . . . . . . . . . . . . . . . . . . . 2
3.3 Weakly singular f................................... 3
3.4 Improving performance by providing stationary points and/or a search radius . . 3
3.5 Visualisation ...................................... 3
4 Examples 3
4.1 Simplephase ...................................... 3
4.2 Singular amplitude function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
4.3 Contourintegral .................................... 4
1 Setup
Our NSD code is available to download from https://github.com/AndrewGibbs/NSDpackage.
Once downloaded, open Matlab, navigate to NSDpackage-master and run StandardPaths.m
to add all of the necessary subfolders to the Matlab path.
2 Overview
This package is designed to solve integrals of the form
Iω[f, g] = Zb
a
f(x)eiωg(x)dx, (2.1)
where
•aand bcan be real numbers. Under certain restrictions, aand bmay denote infinite points
in the complex plane, see §3.2. For a(or b) finite, we represent this by the variable a(or
b). For a(or b) infinite, we represent this by the angle of the infinite point by the variable
a(or b).
•Frequency parameter ω > 0, represented by the variable freq.
•gcomplex analytic in some neighbourhood of [a, b].
•fis non-oscillatory, has a finite number of weak singularities inside of [a, b], otherwise
analytic inside a neighbourhood of [a, b]
Requiring only the derivatives of g, our package will return a quadrature rule
Iω[f, g]≈Qω[f, g] =
N
X
i
wif(zi) (2.2)
where the number of weights and nodes Ndoes not depend on ω, however Nwill typically
increase with the number of stationary points close to the region of integration. The only
1
requirement for our quadrature rule to work is a cell array containing vectorised anonymous
functions corresponding to increasing derivatives of g. Hence,
G{j+1}= g j ;,where g j represents g(j)(z).
For example, if g(x) = x2, then we would have
G = {@(x) x.^2, @(x) 2*x, @(x) 2, @(x) 0};
Each derivative of gshould return an upright vector of values, if the input is also an upright
vector of values. Without this, the code will most likely error. The minimum number of
derivatives required depends on the order of the stationary points, which in general is not
known. Typically four or five derivatives should be more than enough, the user is notified if
more are required. To obtain the weights and nodes of the quadrature rule (2.2), one can type
into Matlab:
[z,w] = PathFinder( a, b, freq, N, G);
This can then be evaluated by typing:
Q=w.’*f(z);
3 Additional options and features
A range of optional arguments can be input to PathFinder. Some of these can improve perfor-
mance, whilst others broaden the type of problems solvable by our routine.
3.1 Polynomial phase
If gis an order mpolynomial,
g(x) = c1xm+c2xm−1+. . . +cmx+cm+1,(3.1)
then the derivatives in G, the stationary points and their orders can all be determined from the
coefficients of g. We can define a vector Cin Cm+1 such that the ith entry is equal to ciof (3.1),
and the routine will compute derivates as required.
[z,w] = PathFinder( a, b, freq, N, [],’gpolycoeffs’,C,...);
3.2 Infinite and semi-infinite contour integrals
If aand/or bare/is an infinite point in the complex plane, then the optional argument of ’ainf’,
and/or ’binf’ flags this. The inputs aand bshould be the angle of the infinite point at which
the contour ends.
In the case of infinite endpoint(s), it is also necessary to provide a radius R>0, such that
outside of this region the steepest descent paths can be assumed to be straight lines. With a
polynomial, this is simply where the leading order term will become dominant. This argument
should be preceeded by ’settledRad’.
In summary, for a fully infinite contour we would type:
[z,w] = PathFinder( a, b, freq, N, G,...,’ainf’,’binf’,’settledRad’,R,...);
where the dots ... denote any other optional arguments discussed in this section.
2
3.3 Weakly singular f
If fhas a weak singularity inside of [a, b], this is handled via Generalised Gaussian Quadra-
ture [2]. The user must create an instance of the singularity class, containing all of the necessary
information about the singularity, which has the inputs singInfo=Singularity1D(position,
blowUpType). Here position denotes the point in Rwhere the singularity occurs, and blowUpType
is a string describing the type of singularity, e.g. blowUpType=’log’. To use this function, the
extra argument should be added like so:
[z,w] = PathFinder( a, b, freq, N, G,...,’fSingularities’,S,...);
where the dots ... denote any other optional arguments discussed in this section.
3.4 Improving performance by providing stationary points and/or a
search radius
By default, the routine will search for them, inside of a rectangle R, for some search radius R
such that
R:= {z∈C|dist(z, [a, b]) < R}.
The default option is R∼1/ω, justified as stationary points which are far from [a, b] will have
little contribution to the integral. Alternatively the user can specify a search radius via the
option ’rectRad’,
[z,w] = PathFinder( a, b, freq, N, G,...,’rectRad’,R,...);
where the dots ... denote any other optional arguments discussed in this section.
The user can provide stationary points if they are known. This will significantly decrease
the time taken to construct the quadrature rule, as searching for stationary points is often the
most computationally expensive part of the method. This can be achieved simply by
[z,w] = PathFinder( a, b, freq, N, G,...,’stationary points’, C, ’order’, D,,...);
where Cis a vector of stationary points, and Dis a vector of their orders (see e.g. [1] for definitions
of these terms).
3.5 Visualisation
Add the optional argument ’visuals on’ to produce a plot describing the process of the NSD
package. The concentric rectangles used to find the roots are visible, alongside the nodes in the
complex plane. Chosen paths are plotted in pink, and discarded paths are plotted in blue. For
example, the following code
G = @(x) x.^3, @(x) 3*x.^2, @(x) 6*x, @(x) 6, @(x) 0;
PathFinder( -1, 1, 10, 15, G, ’visuals on’);
produced Figure 1.
4 Examples
Here we provide a series of example Matlab code for some typical oscillatory integrals.
4.1 Simple phase
The integral
Z1
−1
sin(x)e100ix2dx
3
Figure 1
can be approximated using only a few lines with:
G={@(x) x.^2, @(x) 2*x, @(x) 2, @(x) 0};
[z,w] = PathFinder( -1,1,100,10,G);
Q = w.’*cos(z)
Alternatively, as the phase is a polynomial and we know the only coefficient, we may replace
the first two lines with [z,w] = PathFinder( -1,1,100,10,[],’gpolycoeffs’,[1 0 0]);.
As explained in §3.1, this bypasses the need to search for stationary points. The fourth input to
PathFinder ensures that there will be (roughly speaking) 10 quadrature points at each endpoint,
and 2 ×10 on each path from the stationary point at zero.
4.2 Singular amplitude function
The integral
Z1
0
log(x)e12345ixdx
can be approximated in a few lines:
S=Singularity1D(0, ’log’);
G = {@(x) x, @(x) 1, @(x) 0 };
[ z, w ] = PathFinder( 0,1,12345,10,G, ’fSingularities’,S);
Q = w.’*log(z)
4.3 Contour integral
The following non-oscillatory integral over an infinite contour
Zeiπ/10
ei9π/10
ei(2x2/5−x4/2−x2)dx
can be approximated using five lines:
polyCoeffs=[2/5 -1/2 0 -1 0 0];
4
R = monomialSettleRadius(polyCoeffs);
a=exp( 1i*pi/10);
b=exp( 1i*pi*9/10);
[z,w] = PathFinder( a,b,1,10,[], ’ainf’, ’binf’, ’settleRad’, R, ’gpolycoeffs’, polyCoeffs);
Q = sum(w)
References
[1] Alfredo Dea˜no Gamallo, Daan Huybrechs, and Arieh Iserles. Computing Highly Oscillatory
Integrals, volume 155. SIAM, 2018.
[2] D. Huybrechs and R. Cools. On generalized Gaussian Quadrautre rules for singular and
nearly singular integrals. SIAM J. Numer. Anal., 47:719–739, 2009.
5
