Fig. 17.1 is a schematic of the optimum relationship between h and N . If we keep h

¬xed, and merely use a larger and larger number of terms in (17.12), our approximation

will converge not to f (y) but to fW (y), and we are stuck with the ¬xed error |f (y) ’ fW (y)|

even with million grid points. Conversely, if we only decrease h and truncate the sum in

(17.12) by omitting all grid points such that |yj | > L for ¬xed L, we would be left with

an O(f [L]) error even if the grid spacing is ridiculously small. Because of the need to

√

simultaneously increase L and decrease h, the error is O(exp[’p N ]) for a typical function

like f (y) = sech(y): Subgeometric convergence with exponential index of convergence =

1/2.

There are several alternatives for representing functions on in¬nite or semi-in¬nite in-

tervals, but they do no better. Sub-geometric convergence is a way of life when the interval

is unbounded. The grid of Gaussian quadrature points for the Hermite functions, for ex-

ample, automatically becomes both wider (in total span) and more closely spaced as N

increases as shown in Fig. 17.2. There is always this need to serve two masters ” increas-

ing the total width of the grid, and decreasing the spacing between neighboring points ”

for any method on an unbounded interval. The cardinal function series differs from or-

thogonal polynomial methods only in making this need explicit whereas it is hidden in the

automatic behavior of the interpolation points for the Hermite pseudospectral algorithm.

For more information about sinc numerical algorithms, consult the review article by

Stenger (1981) and the books by Lund and Bowers (1992) and Stenger (1993) or the articles

in Table 17.1.

17.3. WHITTAKER CARDINAL OR “SINC” FUNCTIONS 345

Table 17.1: A Selected Bibliography of Sinc (Whittaker Cardinal) Applications

References Comments

Whittaker(1915) A few words from our founder: theory only

Stenger (1979) Boundary value problems

Lundin (1980) Nonlinear PDE boundary value problem

Stenger (1981) Review article; theory plus numerical applications

Saffman&Szeto(1981) Vortex stability

Lund&Riley(1984) Sinc with mapping to solve radial Schroedinger equation

Elliott&Stenger(1984) Numerical solution of singular integral equations

Robinson&Saffman (1984) Nonlinear eigenvalue problem in 2 space dimensions

Lund (1986) Symmetrization of sinc-Galerkin matrix

Boyd (1986b) Shows that exponential mapping of Stenger to in¬nite

interval is key to solving problems with

endpoint singularities on [’1, 1];

sinc basis is effective, but other

in¬nite interval bases work just as well

Schaffer&Stenger (1986) Multigrid iteration with sinc

Bowers & Lund(1987) Singular Poisson equations

Eggert&Jarratt&Lund(1987) Finite and semi-in¬nite intervals, too, through mapping;

singularities at end of ¬nite interval

Lewis&Lund&Bowers(1987) Sinc in both x and t for parabolic PDE

McArthur&Bowers&Lund(1987a) Second order hyperbolic PDEs

McArthur&Bowers&Lund(1987a) Numerical methods for PDEs

Bialecki(1989, 1991) Boundary value problems

Lund&Bowers&McArthur (1989) Elliptic PDEs

Smith&Bowers&Lund(1989) Fourth order problems in control theory

Jarratt&Lund&Bowers(1990) Sinc basis; endpoint singularities , ¬nite interval

Lund&Vogel (1990) Inverse problem for parabolic PDE

Boyd (1991f) Sinc solution of boundary value problems with a banded

matrix obtained by Euler sequence acceleration

Smith&Bogar&Bowers&Lund(1991) Fourth order differential equations

Lund&Bowers&Carlson(1991) Boundary feedback stabilization (control theory)

Smith&Bowers(1991) Inverse problem for Euler-Bernoulli beams

Smith&Bowers&Lund(1992) Inverse problem for Euler-Bernoulli beams

Boyd (1992c) Fast Multipole Method provides a fast sinc transform

Boyd (1992d) Fast off-grid sinc interpolation

Lund & Bowers (1992) Monograph

Stenger (1993) Monograph

Yin (1994) Poisson equation through tensor-product sinc basis

Koures (1996) Coulomb Schroedinger equation

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

346

Appendix F gives formulas for the derivatives of the Whittaker cardinal functions at the

grid points. The ¬rst derivative matrix is skew-symmetric while the second derivative is

symmetric (Stenger, 1979). These are desirable properties. Skew-symmetry implies that all

eigenvalues of the ¬rst derivative matrix are pure imaginary; thus the sinc discretization

of the one-dimensional advection equation, ut + cux = 0, is nondissipative. Similarly, the

symmetry of the second derivative matrix implies that all eigenvalues are real; the sinc

discretization of the diffusion equation gives a system of ODEs in time whose eigenmodes

all decay monotonically in time.

The sinc expansion is the easiest of all pseudospectral methods to program and to un-

derstand. One weakness, shared with the Hermite functions, is that the Whittaker cardinal

series is useful only when f (y) decays exponentially as |y| ’ ∞. In contrast, the T Bn (y)

will give rapid, exponential convergence for functions which decay algebraically with |y| or

asymptote to a constant, provided certain conditions are met. Thus, the mapped Cheby-

shev polynomials have a wider range of applicability than the sinc expansion.

A more serious weakness of sinc series is that ” unlike the mapped Chebyshev func-

tions ” they cannot be summed via the Fast Fourier Transform. This is especially unfortu-

nate for time-dependent problems and for iterative solution of multi-dimensional bound-

ary value problems.

A third limitation is that sinc functions do not give the exact solution to any classical

eigenvalue problem (to my knowledge). In contrast, the Hermite functions are the ide-

alized or asymptotic solutions to many physical problems. Thus, the Whittaker cardinal

functions, despite their simplicity, have not eliminated the competing basis sets discussed

in the next two sections.

The “Big Sky School” of Stenger, Lund, and their students have applied the sinc expan-

sion to problems on a ¬nite interval with endpoint singularities ” always in combination

with a map that exponentially stretches the bounded domain to an in¬nite interval. As ex-

plained in Chap. 16, Sec. 5, this is a good strategy for coping with endpoint branch points.

Eggert, Jarratt, and Lund (1987) also compute some solutions without endpoint singulari-

ties, but this is not recommended. Chebyshev series are much more ef¬cient when u(x) is

analytic at the boundaries.

In summary, it is best to regard the Whittaker cardinal functions as a good basis for

[’∞, ∞] only. Although it has some important defects in comparison to Hermite and

rational Chebyshev functions, the sinc series is the simplest and it converges just as fast as

the alternatives.

17.4 Hermite functions

The members of this basis set are de¬ned by

ψn (y) ≡ e’0.5 y Hn (y)

2

(17.14)

where Hn (y) is the Hermite polynomial of degree n. Like the familiar Chebyshev poly-

nomials, the n-th Hermite polynomial is of degree n. The members of the set satisfy a

three-term recurrence relation. The derivative of Hn (y) is simply (2 n Hn’1 ), so it is easy to

use Hermite functions to solve differential equations.

Unlike most other sets of orthogonal functions, the normalization constants for the Her-

mite functions increase exponentially with n, so it is dif¬cult to work with unnormalized

functions without encountering over¬‚ow (unless N is very small). In contrast, normalized

Hermite functions satisfy the bound

¯

|ψn (y)| ¤ 0.816 all n, all real y (17.15)

17.4. HERMITE FUNCTIONS 347

Table 17.2: Hermite Function Theory & Applications

References Comments

Hille(1939, 1940a,b, 1961) Convergence theorems

Englemann&Feix 1D Vlasov eq. (plasma physics); Fourier in x

& Minardi&Oxenius(1963) Hermite in velocity coordinate u

Grant&Feix(1967) Fourier-Hermite for 1D Vlasov equation

Fokker-Planck collision term to improve Hermite convergence

Armstrong(1967) 1D Vlasov equation; Fourier-Hermite Galerkin

Hermite truncation reduced with increasing time to mitigate

slow convergence of Hermite part of tensor basis

Birkhoff&Fix(1970) Eigenproblems

Joyce&Knorr 1D Vlasov eq.; Fourier-Hermite basis compared with

&Meier(1971) double Fourier & method-of-moments (!)

Anderson (1973) Time-dependent model of the tropical ocean

Moore & Philander (1977) Equatorial waves and jets in the ocean (review)

Shoucri&Gagn´ (1977)

e 2D Vlasov; tensor Hermite basis in velocity

¬nite differences in space; plasma physics

Bain (1978) Convergence theorems

Banerjee (1978) Eigenvalue problem: quantum anharmonic operator

Banerjee et al. (1978) ψn (±y) with variable ±

Boyd (1980b) Rate of convergence theorems

Boyd (1984a) Asymptotic theory for Hermite coef¬cients through

steepest descent and the calculus of residues

Tribbia (1984) Hough-Hermite Galerkin method

Boyd & Moore (1986) Euler sequence acceleration of algebraically-converging

Hermite series with applications to oceanography

Marshall & Boyd (1987) Equatorial ocean waves

Smith(1988) Hermite-Galerkin model of tropical ocean

Funaro & Kavian (1991) Time-dependent diffusion problems

Weideman (1992) Hermite differentiation matrices are well-conditioned

Holvorcem(1992) Summation methods for Hermite series including

Holvorcem&Vianna(1992) Green™s functions, which depend on two variables

Vianna&Holvorcem(1992)

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

Boyd (1993) Hermite methods in symbolic manipulation languages

Tang (1993) Gaussian functions; role of the scale factor

Holloway(1996a,b) Vlasov-Maxwell equations (plasma physics)

Introduces “asymmetric” Hermite basis, which is very

ef¬cient for this application

Tse&Chasnov(1997) Hermite for vertical coordinate

Convection in unstable layer surrounded by

unbounded stable layers; Fourier in x, y

Schumer&Holloway(1998) 1D Vlasov equation; Hermite-Fourier basis

asymmetric Hermite is unstable; standard Hermite is good

with use of anti-¬lamentation ¬ltering

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

348

where the overbar denotes that the ψn (y) have been normalized. As for Chebyshev and

Fourier series, this bound (17.15) implies that the error in truncating a Hermite series is

bounded by the sum of the absolute values of all the neglected terms.

The optimum pseudospectral points are the roots of the (N +1)-st Hermite function. The

points are not exactly evenly spaced ” the outermost points are a little farther apart than

the separation between points near y = 0 ” but the asymptotic approximations to the Her-

mite functions show that the Hermite grid does not differ very much from the sinc grid

as shown in Fig. 17.2. In the previous section, we found that it was necessary to simul-

taneously (i) decrease the spacing between adjacent grid points and (ii) increase the span

between the outermost points on the grid if the sinc series was to give increasing accuracy

as N increased. Exactly the same √ thing happens with the Hermite grid: the separation be-

√

tween grid points decreases as 1/ N while the span of the grid increases as N . It is not

necessary to choose h and the span of the grid explicitly; using the roots of the Hermite

function as the grid points automatically insures the correct behavior.

Fig. 17.3 illustrates the ¬rst six (normalized) Hermite functions.

However, it is still necessary to pick a scaling factor ± since we can can always use

ψn (±y) as the basis set for any ¬nite ±. This freedom does not exist for a ¬nite interval, but

after a change of scale in y, an in¬nite domain is still an in¬nite domain. Banerjee (1978)

and Banerjee et al.(1978), Boyd (1984a), and Tang(1993) have discussed the choice of scale

factor. The theory (Boyd, 1984a) involves heavy use of steepest descent to asymptotically

approximate the Hermite coef¬cient integrals as n ’ ∞. One general conclusion is impor-

tant since it is also true of mapped Chebyshev methods: The optimum ± is a function of

N . A good choice of scaling parameter for N = 5 is a mediocre choice for N = 50; at least

for the cases studied, ± should increase with N .

The rates-of-convergence theory for Hermite functions is now well-developed for a ¬xed

(N -independent) scale factor ±. It is worth summarizing since it indicates the subtleties in-

troduced by the unboundedness of the expansion interval. We must begin with a de¬nition

because the rate of convergence depends not only on the location of singularities, as true of

all expansions, but also upon the rate at which the function f (y) decays as |y| tends to in¬n-

ity.

Figure 17.2: The pseudospectral grid for Hermite functions. For all N , the grid points are

the roots of the (N +1)-st Hermite polynomial. The grid spacing varies slightly, but is almost

uniform over the interval. The grid closely resembles the sinc grid: when N is quadrupled,

the span of the grid points is roughly doubled and the grid spacing is halved.

17.4. HERMITE FUNCTIONS 349

Figure 17.3: Graphs of the ¬rst 6 Hermite functions, ψn (y).

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

350

De¬nition 39 (ORDER OF REAL AXIS DECAY)

The ORDER of REAL AXIS DECAY k is the least upper bound of j for which

f (y) = O e’p|y|

j

(17.16)

for some constant p as |y| ’ ∞ along the real axis. The function is said to be “SUB-GAUSSIAN”

or “SUPER-GAUSSIAN” respectively if

[“SUB-GAUSSIAN”] (17.17)

k<2

[“SUPER-GAUSSIAN”] (17.18)

k>2

Theorem 33 (HILLE™S HERMITE WIDTH-OF-CONVERGENCE)

(i) The domain of convergence of a Hermite series is the in¬nite STRIP about the real y-axis bounded

by

|Im(y)| = w (17.19)

(ii) For an ENTIRE function, that is, f (y) which has no singularities except at in¬nity, the width is

∞ if k > 1

(17.20)

w=

if k < 1

0

(iii) For SINGULAR functions, let „ denote the absolute value of the imaginary part of the location

of that singularity which is closest to the real y-axis, and let p denote the constant in the de¬nition

of real axis decay. Then

±

if k > 1

„

smaller of [p, „ ] if k = 1 (17.21)