. 73
( 136 .)



if k < 1

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

an ∼ O e’w 2n+1

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)

(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)

(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

an ∼ O (17.28)

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

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
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 ).

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,

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
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
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

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
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 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)


. 73
( 136 .)