it that 1/h and the order of Richardson extrapolation simultaneously

tend to in¬nity, the quadrature scheme converges.

Lyness and Ninham [190] noted this exponential dependence on 1/h

in quadrature errors nearly thirty years ago. Lyness has emphasized

that the Euler-MacLaurin series for the trapezoidal rule is closely relat-

ed to a general formula for the coe¬cients of a Fourier series. His

“FCAE” [Fourier Coe¬cient Asymptotic Expansion] takes the form

[188, 189]

M

’ g 2k’1 (0)

2k’1 (1)

n+k’1 g

an = (’1) (120)

(2πn)2k

k=1

(’1)M +1 1

g (2M +1) (x) sin(2π n[x ’ 1/2] )dx

+

(2πn)2M +1 0

plus a similar series for the sine coe¬cients. It is derived by integration-

by-parts. As noted by Lyness, it is usually divergent.

As explained in the books by Boyd [55] and Canuto et al. [91], spec-

tral methods usually employ a basis set so that the error is exponential-

ly small in 1/h or equivalently, in the number of degrees of freedom M .

Fourier series, for example, are restricted to periodic functions. Cheby-

shev polynomials give exponential convergence even for non-periodic

functions, provided only that f(x) is free of singularities on x ∈ [’1, 1].

These polynomials are de¬ned by Tn (cos[θ]) ≡ cos(n θ) so that the

ActaApplFINAL_OP92.tex; 21/08/2000; 16:16; no v.; p.63

64 John P. Boyd

expansion of f (x) as a Chebyshev series is identical, under this change

of variable, with the Fourier expansion of f (cos(θ)). The transformed

function is always periodic in θ, so the error in truncating a Chebyshev

series after M terms decreases exponentially fast with M . For Cheby-

shev and Fourier spectral methods, the power series in h is always the

trivial one with zero coe¬cients.

It follows that all asymptotic approximations to Fourier, Chebyshev,

and other spectral coe¬cients for large M are implicitly hyperasymp-

totic (Table V). One might object that this catalogue of asymptotics for

orthogonal series is out of place here because it is “beyond all orders”

[in h] only because the coe¬cients of all powers of h are zero. The

main tools for the work of Elliott, Miller, Luke, Weideman and Boyd

were steepest descents and the calculus of residues ” no explicit use

of hyperasymptotic thinking at all.

Nevertheless, Table V makes several important points. First, much

of hyperasympotics is steepest descent and the calculus of residues. In

discussing the Stieltjes function earlier, for example, we noted that one

could go beyond the superasympotic approximation by applying steep-

est descent to the error integral of the optimally-truncated power

series. Similarly, the heart of the PKKS technique of matched asymp-

totic expansions in the complex plane is the notion that the singular-

ities or other critical points closest to the real axis control the hyper-

asymptotic behavior. The asymptotic behavior of the coe¬cients of

orthogonal expansions is likewise controlled by complex singularities.

In both cases, the constant q inside the exponential factor, exp(’q/ )

or exp(’q/h), is simply the distance from the dominant singularity to

the real x-axis.3

Second, an important hyperasymptotic strategy is to isolate the

exponentially small contributions. For the large degree behavior of

Chebyshev and Fourier coe¬cients, this isolation is a free gift, the result

of choosing the sensible spectral basis ” Chebyshev or Legendre poly-

nomials for non-periodic problems, spherical harmonics for problems

on the surface of a sphere and so on. For less trivial problems, the key

to isolation is to subtract an optimally truncated asymptotic expansion.

This is the justi¬cation for applying steepest descent to the error inte-

gral for Stieltjes function, for applying Euler™s method or other sum

acceleration scheme, for Dingle™s universal terminants and for Berry™s

smoothing of Stokes phenomenon. The di¬erence between a quantity

and its superasymptotic approximation is as trivially isolated as the

3

For Chebyshev series, “distance from the real axis” means distance in the trans-

formed coordinate θ, the argument of the equivalent Fourier series, rather than the

argument of the Chebyshev polynomials, x = cos(θ).

ActaApplFINAL_OP92.tex; 21/08/2000; 16:16; no v.; p.64

65

Exponential Asymptotics

spectral coe¬cients of a Fourier series, for which the superasymptotic

approximation is zero.

Third, the asymptotics of spectral series and other numerical pro-

cesses is a relatively under-cultivated area. Can the ideas reviewed here

lead to optimal order Richardson extrapolation for numerical quadra-

ture, various classes of di¬erential equations, and so on?

Fourth, there have been some limited but important excursions

beyond conventional asymptotics in the analysis of the convergence

of spectral series. For example, Boyd™s 1982 article on the optimization

of Chebyshev methods on an unbounded domain notes that there are

two di¬erent species of contributions to the asymptotic spectral coef-

¬cients: (i) saddle point contributions that depend on how rapidly the

function f (x) being expanded decays as | x |’ ∞ and (ii) contribu-

tions from the poles and other singularities of f (x). In the limit that

the degree n of the rational Chebyshev coe¬cient tends to in¬nity for

¬xed value of the “map parameter” L, one type of contribution will be

exponentially large compared to the other, and it is inconsistent (in the

Poincar´ sense of asymptotics) to retain the other. Boyd points out that

e

to optimize numerical e¬ciency, one should allow L to vary with the

truncation N of the Chebyshev series. Convergence is optimized when

N and L simultaneously tend to in¬nity in a certain way so that both

the pole and saddle point contributions are of equal order. This sort

of analysis does not explicitly use exponentially-improved asymptotics

of the Dingle-Berry-Olver sort. Nevertheless, hyperasymptotic thinking

” considering the role of terms that at ¬rst glance are exponentially

small compared to the dominant terms ” is absolutely essential to this

kind of numerical optimization.

A few other interesting studies of the role of exponential smallness

in numerical analysis have already been made. For example, the usual

second order di¬erential equation for the nonlinear pendulum, qtt +

sin(q) = 0, can be written as the equivalent system

pt = sin(q), qt = p (121)

where q is the angle of de¬‚ection of the pendulum with q = 0 when the

pendulum is standing (unstably) on its head and p is the momentum.

Hakim and Mallick [136] show that when the ¬rst equation is discretized

by a forward di¬erence and the second equation by a backward di¬er-

ence, the system (121) becomes what in dynamical systems theory, with

a slight change in notation, is called the “standard mapping”:

pn+1 = pn + „ sin(qn ), qn+1 = qn + „ pn+1 (122)

The usual numerical analysis description begins and ends with the

statement that this algorithm is second order accurate, that is, has

ActaApplFINAL_OP92.tex; 21/08/2000; 16:16; no v.; p.65

66 John P. Boyd

3

3

p

p

0

0

-3

-3

π 2π

0 π 2π 0 q

q

Figure 14. Phase plane for the nonlinear pendulum. Left panel: trajectories for the

exact solution. Right panel: trajectories when the di¬erential equation is integrated

by a ¬nite di¬erence scheme, i. e, the “standard mapping”, with „ = 1. The cross-

hatched area is the region spanned by a single chaotic trajectory. To avoid printer

overload, the individual iterates were erased and replaced by a uniform texture.

an error which is proportional to „ 2 . 4 Hakim and Mallick point out

that there are also changes which are exponentially small in 1/„ ” and

qualitatively di¬erent from the e¬ects of ¬nite di¬erencing which are

proportional to powers of the timestep.

These changes are easiest to explain by examining the trajectories

of the di¬erential equation and the di¬erence system in the phase plane

(Fig. 14). The nonlinear pendulum is an exactly integrable system, and

all trajectories are periodic. The closed curves in the phase plane rep-

resent side-to-side, small amplitude oscillations of the pendulum. The

open curves at top and bottom show trajectories in which the pendu-

lum swings through complete loops like a propeller. These two species of

trajectories are divided by the “separatrix”, which is a trajectory that

passes through q = 0, the unstable equilibrium, with zero momentum

p. The separatrix and trajectories near it are super-sensitive to per-

turbations because only a tiny additional amount of momentum will

su¬ce to push a large amplitude oscillation over the top and thereby

convert it into a propeller-like motion.

The di¬erence system, alias “standard mapping”, is not integrable

and has chaotic solutions. When the time step „ is very small, however,

one would expect that in some sense the di¬erence and di¬erential sys-

4

The one-sided di¬erences are only ¬rst order, but by eliminating p, one can

show that the system is equivalent to applying centered, 2d order di¬erences to the

second order di¬erential equation.

ActaApplFINAL_OP92.tex; 21/08/2000; 16:16; no v.; p.66

67

Exponential Asymptotics

tems would be close to one another. Indeed, trajectories away from the

separatrix are not drastically altered by the discretization; ¬nite di¬er-

ences give a good approximation to these trajectories of the nonlinear

pendulum. The neighborhood of the separatrix, however, dissolves into

chaotic motion. The incoming and outgoing separatrices from the equi-

librium split with a splitting angle at the point where the separatrices

cross of approximately [136], [180]

π2

3514.9

φ∼ exp ’ (123)

„3 „

The width of the region of chaos around the separatrices has a similar

exponential dependence.

The bland statement “the method is second order” or even a formal

expansion of the errors in powers of the timestep „ completely misses

the spawning of this region of chaos. Such exponentially small quali-

tative changes are perhaps less important in the real world, where the

original di¬erential equation probably has regions of chaotic motion

anyway, than in the idealized world of exactly integrable systems such

as the nonlinear pendulum, the Korteweg-deVries equation and so on.

Still, it reiterates the theme that hyperasymptotics is important to

numerical analysis.

Hakim and Mallick observe that the discretized system can be inter-

preted in two ways: (i) a second order accurate approximation to the

pendulum system or (ii) a fourth order accurate approximation to

a nonlinear di¬erential equation which is obtained by modifying the

pendulum equation by the addition of a higher derivative with a „ -

dependent coe¬cient. Although the second interpretation seems rather

arti¬cial, it is also illuminating. Weakly nonlocal solitary waves and

hydrodynamic boundary layers arise in this same way through addi-

tion of a higher derivative with a coe¬cient proportional to the small

parameter. The result of such a singular perturbation is that the power

series in the small parameter is divergent, and there are e¬ects which

depend on the exponential of the reciprocal of the small parameter.

17. Numerical Methods for Exponential Smallness or:

Poltergeist-Hunting by the Numbers, I: Chebyshev &

Fourier Spectral Methods

Because of the messiness of hyperasymptotic methods even for clas-

sical special functions and ordinary di¬erential equations, numerical

algorithms are important both as checks and as alternatives to hyper-

asymptotics. The exponential dependence on 1/ for << 1 cries out

ActaApplFINAL_OP92.tex; 21/08/2000; 16:16; no v.; p.67

68 John P. Boyd

for numerical schemes whose error also falls exponentially with the

number of degrees of freedom M . Fortunately, Chebyshev and Fourier

spectral methods [55, 91] and also Pad´ approximants [9, 19] have this

e

property. In this section, we shall discuss spectral methods while Pad´ e

algorithms are described in the following section.

However, when the unknown function f ( ) has only a divergent pow-

er series and also has contributions that lie beyond all orders in , both

the rate of convergence and (sometimes) the methodology are altered.

For example, when a function f (x) which is analytic on x ∈ [’1, 1] is

expanded in a Chebyshev series, the error decreases geometrically with

M, the truncation of the Chebyshev series. In other words, the error

EM = O(exp(’µM )) as M ’ ∞.

When f ( ) has only a divergent power series about = 0, the Cheby-

shev or other spectral series on an interval that includes this point

will lack geometric convergence. However, as long as the function is

in¬nitely di¬erentiable (with bounded derivatives), it is easy to prove

by integration-by-parts that the error must decrease faster than any

¬nite power of M . “C ∞ ” singularities do not defeat spectral methods,

but merely slow them down.

By using the method of steepest descents, one can often show that