Asymptotic Analysis of Numerical Steepest Descent with Path Approximations

TLDR: A loss of asymptotic order is observed, but in the most relevant cases the overall asymaptotic order remains higher than a truncated asymPTotic expansion at similar computational effort.

Abstract: We propose a variant of the numerical method of steepest descent for oscillatory integrals by using a low-cost explicit polynomial approximation of the paths of steepest descent. A loss of asymptotic order is observed, but in the most relevant cases the overall asymptotic order remains higher than a truncated asymptotic expansion at similar computational effort. Theoretical results based on number theory underpinning the mechanisms behind this effect are presented.

Figures

Table 1: Measured order for different numbers n of Gauss-Laguerre points with m terms in the Taylor expansion of the steepest descent path. First column is with the exact path. Numbers in parentheses are orders predicted in Theorem 4.3.

Table 2: Measured order for different numbers n of Gauss points with m terms in the expansion of the steepest descent path. Numbers in parentheses are orders predicted in Theorem 4.3. a)-Case of order one stationary point, b)-Case of order two stationary point.

Figure 1: The contours of the imaginary part of the oscillator g(x) = x2 in the complex plane and the corresponding paths of steepest descent. Two paths emerge from the endpoints x = −1 and x = 1. They are connected by a path passing through the stationary point at x = 0.

Figure 2: Illustration of the exact (continuous) and approximate (dashed) steepest descent paths. The curve Γ connects a truncation of these two paths.

Figure 3: Log-plot of error for different path approximations. Case of regular endpoint.

Full text

Asymptotic analysis of numerical
steepest descent with path
approximations
Andreas Asheim and Daan Huybrechs
Report TW 536, March 2009
n
Katholieke Universiteit Leuven
Department of Computer Science
Celestijnenlaan 200A – B-3001 Heverlee (Belgium)

Asymptotic analysis of numerical
steepest descent with path
approximations
Andreas Asheim and Daan Huybrechs
Report TW 536, March 2009
Department of Computer Science, K.U.Leuven
Abstract
We propose a variant of the numerical method of steepest de-
scent for oscillatory integrals by using a low-cost explicit polynomial
approximation of the paths of steepest descent. A loss of asymp-
totic order is observed, but in the most relevant cases the overall
asymptotic order remains higher than a truncated asymptotic ex-
pansion at similar computational effort. Theoretical results based
on number theory underpinning the mechanisms behind this effect
are presented.
Keywords : oscillatory integral, steepest descent, numerical integration
AMS(MOS) Classification : Primary : 65D30, Secondary : 30E20, 41A60.

Asymptotic analysis of numerical steepest descent
with path approximations
Andreas Asheim
∗
, Daan Huybrechs
†
March 16, 2009
Abstract
We propose a variant of the numerical method of steepest descent
for oscillatory integrals by using a low-cost explicit polynomial approx-
imation of the paths of steepest descent. A loss of asymptotic order
is observed, but in the most relevant cases the overall asymptotic or-
der remains higher than a truncated asymptotic expansion at similar
computational effort. Theoretical results based on number theory un-
derpinning the mechanisms behind this effect are presented.
AMS(MOS) Classification: Primary : 65D30, Secondary : 30E20, 41A60.
Keywords: oscillatory quadrature, steepest descent, numerical integration.
1 Introduction
Consider a highly oscillatory integral of the form
I[f] =
Z
1
−1
f(x)e
iωg(x)
dx, (1.1)
where ω is a large parameter and f and g are smooth functions called the
amplitude function and oscillator of the integral respectively. Such inte-
grals, often referred to as Fourier-type integrals, appear in a wide area of
applications, e.g., highly oscill atory scattering problems in acoustics, elec-
tromagnetics or optics [5, 3, 13, 2]. Numerical evaluation of Fourier-type
integrals with classical techniques becomes expensive as ω becomes large,
which corresponds to a highly oscillatory integral. Typically, a fixed number
of evaluation points per wavelength is required to obtain a fixed accuracy,
which makes the computational effort at least linear in ω [6].
∗
Department of Mathematical Sciences, NTNU, 7491 Trondheim, Norway. E-mail:
Andreas.Asheim@math.ntnu.no
†
Department of Computer Science, KU Leuven, 3001 Leuven, Belgium. E-mail:
Daan.Huybrechs@cs.kuleuven.be. This author is a PostDoctoral Fellow of the Research
Foundation - Flanders (FWO).
1

Asymptotic techniques on the other hand yield approximations that be-
come more accurate as ω increases, making them superior for ω sufficiently
large. One of these techniques, the principle of stationary phase [20, 25],
states that I[f] asymptotically depends only on f and g in a set of special
points as ω → ∞. These points are the endpoints, here x = −1 and x = 1,
and stationary points - points where the derivative of g vanishes. At sta-
tionary points the integral is locally non-oscillatory. The integral has an
asymptotic expansion in inverse powers of ω, with coefficients that depend
on the derivatives of f and g at these critical points [15].
x=−1 x=1
Figure 1: The contours of the imaginary part of the oscillator g(x) = x
2
in
the complex plane and the corresponding paths of steepest descent. Two
paths emerge from the endpoints x = −1 and x = 1. They are connected
by a path passing through the stationary point at x = 0.
A set of particularly effective ways of obtaining the contribution from
a special p oi nt are the saddle point methods[25, 19, 8]. Based on Cauchy’s
integral theorem, the path of integration can be deformed into the complex
plane without changing the value of the integral, provided that f and g
are analytic [9]. The method of steepest descent is obtained by following a
path where g has a constant real part and increasing imaginary part, which
renders the integral (1.1) non-oscillatory and exponentially decreasing. This
procedure yields separate paths originating from each special point that
typically connect at infinity (see Figure 1 for an illustration). The result
is separate contributions corresponding to each special point. Every one of
these contributions is a non-oscillatory integral that can be written as
Z
∞
0
ψ(q)e
−ωq
r
dq, (1.2)
2

where ψ is a smooth function, r = 1 for endpoint contributions, and r > 1
for stationary p oints. These integrals are usually treated with standard
asymptotic techniques like Watson’s Lemma. The larger class of saddle
point methods also contains methods that follow other paths with similar
characteristics as the steepest descent paths, e.g., Perron’s method[25].
The asymptotic expansion of I[f] in general diverges, but it can yield
very accurate approximations if ω is very large. Still, divergence implies that
the error is uncontrollable, which is problematic in the context of numeri-
cal computations. Recent research has however produced several numerical
methods that exhibit convergence. The Filon-type met hods [15, 14, 16] are
based on polynomial interpolation of the amplitude f and can deliver errors
that are O(ω
−p
) for any p, much like truncated asymptotic expansions, but
with controllable error for fixed ω. Filon-type methods require that mo-
ments w
k
= I[x
k
] are available, a serious drawback in some cases. Combin-
ing asymptotic expansions and Filon-type methods[1] can economise on, but
not eliminate the need for moments. Me thods that do not rely on moments
are the Levin-type methods, due to Levin[18] and extended by Olver[22, 23].
Levin-type methods do not work in the presence of stationary points, but a
work-around is provided in [21]. We refer the reader to [11] for a detailed
overview of these and other numerical methods.
One of the al ternatives is the numerical method of steepest descent [12],
which is a numerical adaptation of the above described method of steepest
descent. Relying on classical numerical integration methods applied to an
exact decomposition of the integral, the numerical method of steepest de-
scent has controllable error wherever the exact decomposition is available,
and asymptotic error decay O(ω
−p
) for any p. The paths of steepest descent
can however be difficult to compute, as their computation corresponds to
solving a non-linear problem that can in practice only be solved iteratively.
The metho d of this paper is similar in spirit but based on the practical
observation that the exact choice of path is not essential. This observation
resonates with the theory behind saddle point methods. A Taylor expansion
of the path of steepest descent, which can explicitly be derived from a Taylor
expansion of the oscillator function g, is in many cases sufficient. Iterative
methods to solve a non-linear problem can therefore be entirely avoided. We
obtain a numerical scheme which is relatively simple to implement and cheap
to evaluate. The method exhibits high asymptotic order, and the order is
in fact higher than one would get from a truncated asymptotic expansion
using exactly the same number of derivatives of g.
It is the purpose of this paper to analyse the asymptotic order of the pro-
posed explicit numerical saddle-point method. Unlike the numerical adap-
tation of the steepest descent method and the other methods for highly
oscillatory integrals mentioned above, the asymptotic order does not follow
from standard results in asymptotic analysis. A seemingly irregular relation
between the number of derivatives of g that are used and the number of
3

Frequently asked questions

Q1. What have the authors contributed in "Asymptotic analysis of numerical steepest descent with path approximations" ?

The authors propose a variant of the numerical method of steepest descent for oscillatory integrals by using a low-cost explicit polynomial approximation of the paths of steepest descent. 

Q2. What is the function f( tr )t r?

Since fξ( tr ω )t r−1 is a smooth function, the integral can be efficiently approximated by Gaussian quadrature withweight function e−t r . 

Q3. What is the result of applying an n-point Gaussian quadrature?

The result of applying an n-point Gaussian quadrature leads to an approximation with an error which is O(ω−(2n+1)/r) as ω → ∞ [7]. 

Q4. What is the method of steepest descent?

The method of steepest descent is obtained by following a path where g has a constant real part and increasing imaginary part, which renders the integral (1.1) non-oscillatory and exponentially decreasing. 

Q5. What is the cost of evaluating a Fourier-type integral?

Numerical evaluation of Fourier-type integrals with classical techniques becomes expensive as ω becomes large, which corresponds to a highly oscillatory integral. 

Q6. What is the cost of a fixed number of evaluation points per wavelength?

a fixed number of evaluation points per wavelength is required to obtain a fixed accuracy, which makes the computational effort at least linear in ω [6]. 

Q7. How do the authors test the conclusion of the theorem?

The authors shall test the conclusion of the theorem by using approximate paths with different number of terms and different number of quadrature points, and then measuring the asymptotic order by regression for each combination. 

Q8. How can the authors eliminate the error contribution from the right endpoint?

By using the exact path and a large number of quadrature points, the authors can nearly eliminate the error contribution from the right endpoint. 

Q9. What is the rewrite of the integral I[f ; hx,?

The authors once again rewrite the integral I[f ; h̃x, P ] in the following form:I[f, h̃x, P ] = e iωg(x)∫ Q0 ψ̃(q)eiωRr+m−1(q)e−ωq r dq, (4.4)where r − 1 is the order of the point x, Rβ(q) is a function of the formRβ(q) = q β∞ ∑j=0rjq j , (4.5)andψ̃(q) = r f̃x(q r) qr−1,is a smooth function independent of ω. 

Q10. How do the authors define the steepest descent integral?

The authors define the local path h̃x by truncating the series of hx after m terms,h̃x(p) = x+ m−1 ∑j=1ajp j . (3.4)This means that the left and right hand side of (2.1) match up to order m,g(h̃x(p)) = g(x) + ip+ O(p m), p→ 0. (3.5)¿From this path the authors can define the steepest descent integral with an approximated path, using the notation f̃x(p) = f(h̃x(p))h̃ ′ x(p) and g̃x(p) = g(h̃x(p)),I[f ; h̃x, P ] =∫ P0 f̃x(p)eiωg̃x(p)dp. (3.6)The authors shall later evaluate this integral numerically. 

Q11. What is the numerical method of steepest descent?

The paths of steepest descent can however be difficult to compute, as their computation corresponds to solving a non-linear problem that can in practice only be solved iteratively. 

Q12. What is the rationale for solving the inverse function g1?

This is a non-linear equation and solving it amounts to computing the inverse function g−1, which in practical applications may be difficult to achieve. 

Q13. How do the authors get the first few coefficients?

The first few coefficients are given explicitly by, with evaluation in x implied,a1 = ig′ , a2 =12g′′(g′)3 ,a3 = i121(g′)5( 2g′g(3) − (g′′)2 ) , (3.3)a4 = − 1241(g′)7( g(4)(g′)2 + 10g′g′′g(3) − 15(g′′)3 )