but-nonananalytic on the expansion are “subgeometric with exponen-

tial index of convergence β”:

an ∼ ν(n) exp ’µnβ (124)

where ν is a prefactor that varies algebraically rather than exponen-

tially with n. Elliott [124, 125] pioneered this method, but Miller [227]

¬rst applied it to estimate Chebyshev coe¬cients for functions with

divergent power series about one endpoint of the expansion interval.

For example, denoting the coe¬cients of the corresponding divergent

power series by bn and the Chebyshev expansion interval by ∈ [0, γ],

Miller found for S( ) and (( )), respectively,

n2/3

16π 1

’ an ∼ (’1) exp ’3 1/3

n n

bn = (’1) n! exp

3γ 3γ γ

(125)

16π 1

’ an ∼ i (’1)n exp ’ —

bn = n! (126)

3γ 3γ

33/2 n2/3 3 n2/3

exp ’i exp ’ 1/3

2 γ 1/3 2γ

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

69

Exponential Asymptotics

0

10

-2

10

-4

10

-6

10

-8

10

-10

10

0 10 20 30 40 50

n

Figure 15. Chebyshev coe¬cients for the expansion on ∈ [0, 1] of the Stieltjes

function, S( ) [lower curve, circles] and [S(’ )] (upper curve, solid).

(One can show that the error in truncating an exponentially convergent

series after the N -th term is proportional to aN as explained in [55, 73].)

Note that even when the asymptotic series is monotonic, corresponding

to a principal value integral with a simple pole at t = ’1/ on the path

of integration, the Chebyshev series still happily converges. However,

roughly twice as many terms are needed to achieve the same accura-

cy for (S(’ )) as for S( ) (Fig. 15). Asymptotic approximations to

the subgeometrically decreasing Chebyshev coe¬cients for many other

special functions are given by Miller [227] and in the books by Luke

[187] and N´meth [236].5

e

It easy to derive asymptotic approximations for Chebyshev or other

spectral coe¬cients for speci¬c functions, but few general results are

known. One such theorem applies to the class of Stieltjes functions,

which is de¬ned to be the set of all functions that can be written in

the form

∞ ρ(t)

f( ) ≡ dt (127)

1+ t

0

5

It appears from his English-language monograph that N´m´th independently

ee

derived many asymptotic approximations in Hungarian-language articles in the mid-

60™s.

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

70 John P. Boyd

for some positive semi-de¬nite weight function ρ(t) such that the moment

integrals

∞

bn ≡ tn ρ(t) dt (128)

0

exist for all nonnegative integers n. (These moments are also the coef-

¬cients of the power series expansion of f( ) about = 0.) W. G. C.

Boyd has developed a general theory for hyperasymptotics for Stieltjes

functions.

J. P. Boyd [49, 51] showed that if the power series coe¬cients diverge

as (n!)r or (rn)!, i. e.,

log | bn |

lim sup =r (129)

n log n

n’∞

then the Chebyshev coe¬cients satisfy the inequalities

log | (log | an |) |

2 r

≥ lim sup ≥1’ (130)

r+2 log n 2

n’∞

Less precisely, the theorem implies that if the Chebyshev coe¬cients

are decreasing like O(exp(’µnβ )) for some constants µ and β, then β

must be smaller 2/(r + 2), implying subgeometric convergence for all

r > 0, i. e. all factorial divergence. However, the exponential index of

convergence β cannot be smaller than 1 ’ r/2.

The integration-by-parts theorem shows that even for r > 2, the

convergence of Chebyshev series is beyond all orders in 1/N . However,

it would be highly desirable to extend Boyd™s theorem to more general

classes of asymptotic functions, and perhaps sharpen it, too.

Berry [29] has shown that his error-function smoothing of Stokes

phenomenon applies even to several broad classes of functions whose

coe¬cients diverge faster than any factorial, so-called “superfactorial

asymptotics”. His numerical illustration is the function

∞

n2

BeS( ; A) ∼ n

exp

A

n=0

∞ exp ’A (t ’ (1/2) log( ))2

A

≡ PV dt(131)

1 ’ exp(2 t)

π ’∞

where we have replaced his z by ’(1/2) log( ) so that the asymptotic

series is a conventional power series. Fig. 16 shows that Chebyshev poly-

nomials have no problem in providing a highly accurate approximation

even though the power series coe¬cients are blowing up like Gaussians

of n. Unfortunately, there are a hundred papers on the asymptotics and

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

71

Exponential Asymptotics

0

10

-2

10

-4

10

-6

10

-8

10

-10

10

0 10 20 30 40

n

Figure 16. Solid with circles: The absolute values of the Chebyshev coe¬cients an

for the expansion of Berry™s superfactorial function, BeS( ; A = 4) on ∈ [0, 1].

The thin dashed curve, closely tracking the solid curve, illustrates the coe¬cients of

exp(’2/ ), which are known to decay as exp(’3 21/3 n2/3 ).

hyperasymptotics of the con¬‚uent hypergeometric function for every

one on the asymptotics of spectral series.

When exponentially small e¬ects are present, there are often algo-

rithmic challenges present, too. Numerical checking of the prefactors

in front of exp(’q/ ) as obtained by the PKKS matched asymptotics

(or whatever) may be impossible in single precision because the expo-

nentially small quantity may fall below the single precision threshold

for moderate , making it impossible to determine whether numerical

di¬erences are due to errors in the asymptotics, or the neglected e¬ects

of higher order terms at not-so-small . When ± was smaller by a factor

of 10’48 than the core of the solitary wave, Boyd [71] numerically com-

puted the radiation coe¬cient ± to a relative precision of six decimal

places by using 70 decimal place arithmetic in Maple.

These calculations would seem to be very expensive since (i) spectral

methods generate full matrices instead of the sparse matrices produced

by ¬nite di¬erences and (ii) multiple precision arithmetic is excruci-

atingly slow in comparison to single precision operations, which are

executed directly in silicon. However, most spectral solutions to bound-

ary value and eigenvalue problems are performed with preconditioning.

That is to say, the bulk of the code for solving a di¬erential equation in

multiple precision with spectral accuracy is to write a low order ¬nite

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

72 John P. Boyd

di¬erence or ¬nite element code to solve the inhomogeneous version of

the problem in single precision. By repeatedly calling this ¬nite dif-

ference solver, evaluating the residual in multiple precision with spec-

tral methods at the end of each iteration, one can obtain a multiple

precision, spectrally accurate solution without ever factoring (or even

computing) the full, dense spectral matrix. By use of the Fast Fourier

Transform, the spectral evaluation of the residual of an ordinary di¬er-

ential equation can be performed at a cost that grows as O(N log2 N )

operations where N is the number of degrees of freedom.

Another special di¬culty that arises mostly in exponentially small

phenomena is that of solutions on the in¬nite interval which do not

decay to zero for | x |,but rather to sinusoidal oscillations. Two good

strategies have been developed.

The ¬rst is to approximate the in¬nite interval by a large but spa-

tially periodic interval, and then expand the solution as a Fourier series.

The drawback is that the radiation coe¬cient ± is sensitive to the spa-

tial period P (modulo the wavelength of the far ¬eld oscillations, W ).

However, the periodic solutions themselves are often interesting. (In

the atmosphere, for example, the solutions are always periodic in lat-

itude and longitude.) In addition, the parameter P/W is actually a

manifestation of a genuine degree of freedom, the “phase factor” ¦,

of the in¬nite interval. Consequently, it is possible to trace the entire

parameter space for the unbounded domain by using the device of a

large but periodic computational interval.

The second strategy is add one or more additional basis functions

which are chosen to mimic the required asymptotic behavior of the

“wings” of the weakly nonlocal solitary wave (or whatever). When

the width of the core structure is inversely proportional to and the

wavenumber of the wing oscillations is kf , an e¬ective radiation basis

function is

φrad (x) ≡ H(x+¦; ) sin(kf x+¦)+H(’x+¦; ) sin(kf ’x+¦) (132)

where a smoothed approximation to the step function is de¬ned by

H(x; ) ≡ (1/2){1 + tanh( x)} (133)

Boyd [60] has successfully applied a mixed basis of the rational Cheby-

shev polynomials[53] plus a single “radiation function” to compute

quantum scattering in one dimension. Boyd [62] shows that one can

construct a basis function that depends nonlinearly on its coe¬cient

” which is an approximation to ± ” so as to apply this method even

when the oscillations for large | x | are allowed to weakly self-interact,

as is inevitable for nonlinear di¬erential equations.

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