if k < 1

0

Theorem 34 (HERMITE RATE-OF-CONVERGENCE)

(i) If the strip of convergence has a ¬nite width w, then

√

an ∼ O e’w 2n+1

(17.22)

,

that is, subgeometric convergence with exponential convergence index r = 0.5.

(ii) For entire functions, the exponential convergence index r is

= k/(2 [k ’ 1]) (SUPER-GAUSSIAN, that is, k > 2) (17.23)

r

(SUB-GAUSSIAN, that is, k < 2) (17.24)

r = k/2

(iii) Geometric convergence is possible only for entire functions and only for these two cases:

(a) f (y) has k = 2 [Gaussian decay; p need not equal 1/2] (17.25)

(b) f (y) is Super-Gaussian (k > 2) and the scaling factor ± is varied (17.26)

with the truncation N as described in Boyd (1984a)

17.4. HERMITE FUNCTIONS 351

(iv) If the function f (y) decays ALGEBRAICALLY with y so that

f (y) ∼ O 1/|y|δ |y| ’ ∞

as (17.27)

and if the function has p continuous derivatives on [’∞, ∞] (with the (p+1)-st derivative contin-

uous except at a ¬nite number of points), then

1

an ∼ O (17.28)

nq/2

where

q = smaller of [p + 3/2, δ ’ 5/6] (17.29)

Eq. (17.28) only gives a lower bound on the convergence rate; it is known that the Hermite

coef¬cients for f (y) ≡ 1 converge as O(n’1/4 ) even though (17.28) can guarantee only that the

coef¬cients will diverge no faster than n5/12 .

Part (i) was proved by Hille (1939, 1940a,b), (ii) and (iii) by Boyd (1984a), and (iv) by Bain

(1978).

Whew! These theorems are rather complicated, but they do illustrate a couple of im-

portant themes that seem to apply to other types of expansions, too. First, the rate of

convergence depends on how rapidly or how slowly f (y) decays as |y| ’ ∞. When k > 1,

that is to say, when the function is decaying faster than an exponential whose argument is

linear in y, the singularities will determine the width of the strip of convergence, just as for

a Chebyshev series on a ¬nite interval. When k < 1, however, the slowness of the real axis

decay has a more devasting effect on convergence than do the poles and branch points;

when k = 1, either the rate of decay or the location of the nearest pole may determine the

width of the strip of convergence.

Second, for entire functions the convergence is fastest when f (y) is decaying at the same

rate as the exp(’0.5y 2 ) factor built into the basis functions: Supergeometric convergence is

possible for this case. The series converges more and more slowly, however, as the rate of

decay varies more and more from that of exp(’0.5y 2 ), regardless of whether the variation

is to decay more rapidly (as a “Super-Gaussian”) or more slowly (“Sub-Gaussian”).

Third, functions which decay algebraically with y have Hermite series whose coef¬-

cients decrease algebraically with n. (Note that (17.28) gives only a lower bound on the

algebraic index of convergence, so further work is needed). It is striking that whereas the

continuity of one additional derivative would increase the algebraic index of convergence

by 1 for a Chebyshev series, it increases the algebraic index by only (1/2) for a Hermite

√

series ” in other words, we gain only a factor of 1/ N for each additional continuous

derivative for Hermite. This is another penalty of the unboundedness of the interval. The

√

spacing between adjacent collocation points decreases only as 1/ N , so (17.28) is not too

surprising.

Weideman(1992) shows that the spectral radius of the Hermite differentiation matrices

√

is O( N ) for the ¬rst derivative and O(N ) for the second derivative. This is equivalent

to the statement that the maximum timestep for an explicit time-marching algorithm for

a PDE which is second order in space is „max ∼ O(1/N ). The reason is that the Hermite

collocation points, which are almost uniformly spaced as for a Fourier or ¬nite difference

scheme, have an average grid spacing which increases only as N 1/2 instead of N because

the span of the grid points is also increasing as O(N 1/2 ). If we solved the same problem

using Fourier domain truncation with the domain size L increasing as O(N 1/2 ), then the

Fourier grid would be almost indistinguishable from the Hermite grid and the maximum

timestep would also decrease as slowly as O(1/N ).

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

352

Thus, for the same N and for grid points spanning the same interval [’L, L], Hermite

algorithms, the sinc method, Fourier domain truncation, and ¬nite differences with domain

truncation all have similar timestep limits.

Weideman(1992) has proved that Hermite differentiation matrices have additional prop-

erties that are numerically desirable. Both the Galerkin and collocation ¬rst derivative

matrices have purely imaginary eigenvalues while the second derivative matrices have

only real, negative eigenvalues. All these matrices are well-conditioned. In addition, the

Galerkin matrices are banded and are either skew-symmetric (¬rst derivative) or symmet-

ric (second derivative) for the orthonormal basis.

Although we can evaluate the Hermite functions and their derivatives through a three-

term recurrence ” hardly a major programming challenge ” the grid points must be

looked up in a mathematical handbook. Consequently, Hermite functions are not as easy

to program as sinc or rational Chebyshev functions.

The latter have the major advantage that they may be evaluated via the Fast Fourier

Transform. Orszag(1986), Boyd(1992c) and Dutt and Rokhlin(1993) have devised fast trans-

forms applicable to Hermite functions. However, these are much more expensive than a

standard FFT (at least a factor of four for N ∼ O(100)), and none of these algorithms has yet

been applied to any purpose except self-demonstration. The inapplicability of the standard

FFT is another black mark against Hermite functions.

One virtue of Hermite functions is that these are (i) the exact eigenfunctions of the quan-

tum mechanical harmonic oscillator and (ii) are the asymptotic eigenfunctions for Mathieu™s

equation, the prolate spheroidal wave equation, the associated Legendre equation, and

Laplace™s Tidal equation. The reason for this special role is that the Hermite functions

solve the eigenvalue problem

uyy + » ’ y 2 u = 0 (17.30)

u(’∞) = u(∞) = 0

where » is the eigenvalue. If we replace (17.30) by the more general equation

uyy + » ’ y 2 + p(y) u = 0, (17.31)

the exponential decay u(y), forced by the ’y 2 term, makes p(y) irrelevant, at least if this

perturbation is small near y = 0. It follows that (17.30) is “generic” in the sense that its

solutions will approximate those of (17.31) for small but otherwise arbitrary p(y).

This close connection with the physics makes Hermite functions a natural choice of

basis functions for many ¬elds of science and engineering. One modern area of great in-

terest is equatorial oceanography and meteorology. One or two Hermite functions give

very accurate approximations to equatorially-trapped waves when linearized about a rest-

ing state. The same functions have been very valuable in solving equatorial jet and wave

re¬‚ection problems as reviewed in Moore and Philander(1977).

The Vlasov-Maxwell equation of plasma physics has been another important applica-

tion. The velocity u is not an unknown, as in conventional hydrodynamics, but is an inde-

pendent coordinate. This is important because the usual ¬‚uid advection term is multiplica-

tion by u in the Vlasov equation, and the Galerkin representation is a banded matrix, elimi-

nating a slow Hermite grid-to-spectral transform. The u-derivatives also generate a sparse

matrix representation because the derivative of the j-th Hermite function is a weighted

sum of the (j ’ 1)-st and (j + 1)-th Hermite functions. As a result, the Vlasov equation ”

although variable coef¬cient ” may be integrated at a cost which is linear in the Hermite

truncation N , versus a quadratic dependence when Hermite functions are applied to the

Navier-Stokes equations.

The irony is: Although the close connection between simple models and Hermite func-

tions motivated these numerical studies in equatorial ¬‚uid dynamics and plasma physics,

17.5. SEMI-INFINITE INTERVAL: LAGUERRE FUNCTIONS 353

there are serious problems with straightforward application of the Hermite Galerkin method.

In equatorial oceanography, the dif¬culty is that for many important problems such as

Yoshida™s 1959 model of an equatorial jet and Moore™s theory for re¬‚ection from an east-

ern boundary, the ¬‚ow decays algebraically with latitude so that the Hermite coef¬cients

decay algebraically, too. This slow convergence can be ¬xed by applying sum acceleration

methods (Boyd and Moore, 1986, and Holvorcem, 1992).

Similarly, early calculations in plasma physics were bedeviled by ¬lamentation in the

velocity coordinate, which necessitated N >> 100 to compute the solution for large times.

This was solved by a ¬lter introduced by Klimas(1987) and Klimas and Farrell(1994), who

used a tensor product Fourier basis instead of a Fourier-Hermite scheme. Holloway(1996a,b)

and Schumer and Holloway(1998) showed that the convergence rate is greatly improved

by scaling the Hermite functions by a constant. (This is a good strategy for general applica-

tions, too, as suggested by Boyd(1984a) and Tang(1993).) The scaled-and-¬ltered Hermite

algorithm is very effective for the Vlasov equation as shown by Holloway and Schumer.

The irony is that both modi¬cations pull the numerical solution away from the un¬ltered

and unscaled analytical approximation which motivated the choice of Hermite functions

in the ¬rst place.

17.5 Laguerre Functions: Basis for the Semi-In¬nite

Interval

The Laguerre functions are a complete orthogonal basis for the domain x ∈ [0, ∞]. They are

close cousins of the Hermite functions and have a similar form: the product of a decaying

Table 17.3: A Selected Bibliography of Laguerre Function Theory & Applications

References Comments

My Name Is Legion Many applications in quantum chemistry

Schwartz(1963) Effects of singularities on rate-of-convergence

Francis(1972) Vertical structure in weather forecasting model

using Ln , not Laguerre function, as basis.

Showed this scheme needs prohibitively short timestep

Hoskins(1973) Comment on Francis; Legendre basis requires a

constraint which Laguerre does not

Gottlieb&Orszag(1977) Wave equation and heat equation

Maday&Pernaud-Thomas Theory; hyperbolic PDEs

&Vandeven (1985)

Phillips&Davies(1988) Semi-in¬nite spectral elements

Mavriplis(1989) Legendre spectral elements coupled with

one semi-in¬nite element with Laguerre

Coulaud&Funaro Elliptic PDEs in exterior domains

&Kavian(1990)

Funaro (1990b,1992b) Theory

Boyd (1992c) Fast Multipole Methods provide fast Laguerre transform

Iranzo&Falqu´ s (1992)

e Comparisons with rational Chebyshev T Ln

Good examples & attention to practical details

DeVeronico & Funaro Spectral elements for wave propagation with Laguerre basis

& Reali(1994) for the end elements, which extend to in¬nity

Khabibrakhmanov&Summers(1998) Generalized Laguerre for nonlinear equations

Black(1998) spectral element alternative to Laguerre: employs the mapping

of Boyd(1987b) from semi-in¬nite domain into ¬nite interval

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

354

exponential with a polynomial:

φn (x) ≡ exp(’x/2) Ln (x), n = 0, 1, . . . (17.32)

where Ln (x) is the n-th Laguerre polynomial. Three-term recurrence relations for evaluat-

ing both the basis functions and their derivatives are given in Appendix A.

If a function f (x) decays exponentially as |x| ’ ∞, then the Laguerre series for f will

usually converge exponentially fast within the largest parabola, with focus at the origin,

which is free of singularities of f (x). However, it seems likely that very rapidly growing

or decaying functions, such as an exponential-of-an-exponential, may have series that con-

verge only on the positive real x-axis. However, we have not seen any investigation of this

in the literature.

We will not discuss Laguerre functions in detail because (i) their theory is very similar to

Hermite functions and (ii) there have been few applications of Laguerre functions outside

of quantum mechanics as indicated by the brevity of Table 17.3. Within quantum chemistry,

however, there have been lots and lots of applications.

The reason is that the Laguerre functions are the exact radial factors for the eigen-

functions of the hydrogen atom. In the Linear Combinations of Atomic Orbitals (LCAO)

method, basis functions are constructed by using hydrogenic orbitals centered on the ¬rst

atom plus hydrogenic orbitals centered on the second atom and so on. Thus, Laguerre

functions (multiplied by spherical harmonics) are the building blocks for many Galerkin

computations of molecular structures as explained in Chapter 3. The hydrogen s-orbitals,

+

which were employed to approximate the H2 ion in Chapter 3, are actually just the lowest

Laguerre functions, φ0 (r/[a0 /2]) where a0 is the Bohr radius of the hydrogen atom.

Just as for Hermite functions, the Laguerre basis always implicitly contains a scale fac-

tor “ in this case, the Bohr radius.

Iranzo and Falqu´ s (1992) give detailed comparisons between Laguerre functions and

e

the rational Chebyshev functions for the semi-in¬nite interval (Sec. 11 below). Some theory

is given in Funaro™s (1992) book and his two earlier articles listed in Table 17.3.

Mavriplis (1989) is unusual in employing Laguerre polynomials “ without the expo-

nential “ as the basis. This unorthodoxy makes it easier to match the Laguerre series with

Legendre polynomial series in a spectral element (domain decomposition) computation.

However, it is has the disadvantage that the absolute error of her solution is unbounded as

x ’ ∞.

This seems disastrous, but really is not. The true is that no standard basis can resolve

an in¬nite number of oscillations with a ¬nite number of basis functions. The Hermite

or Laguerre function or rational Chebyshev series for a function like exp(’x) cos(x) must

inevitably have large relative error at in¬nity because the function continues to oscillate for

all x whereas a truncated sum of basis functions must have only a ¬nite number of zeros,

and will therefore vary monotonically for some suf¬ciently large x. However, the absolute

error is small everywhere because both the exact solution and the numerical approximation

(except for Mavriplis™ scheme) are exponentially small where the relative error is not small.

The Laguerre polynomial series is guilty of only making this breakdown in relative er-

ror, which is hidden in other methods, obvious through a large absolute error. This break-

down is also present, though implicit, with all the other standard basis sets.

The only way out is to add special basis functions to a standard series as explained in

Sec. 13.

17.6. NEW BASIS SETS VIA CHANGE OF COORDINATE 355

17.6 Mapping Methods: New Basis Sets Through a Change

of Coordinates

By using a transformation that maps an in¬nite interval into a ¬nite domain, it is possible

to generate a great variety of new basis sets for the in¬nite interval that are the images

under the change-of-coordinate of Chebyshev polynomials or Fourier series. This strategy

is merely a generalization of the particular map x = cos(t), which transforms the cosine

functions into the Chebyshev polynomials.

An in¬nite variety of maps is possible, but one can make some useful points by consid-

ering three particular classes of mappings, each illustrated by an example:

1. “Logarithmic maps”:

x ∈ [’1, 1]

y = arctanh(x), (17.33)