. 31
( 136 .)


the continuous spectrum or “continuum”. Since the sech2 potential decays exponentially
fast to zero as |y| ’ ∞, the continuum eigenmodes are asymptotically proportional to
exp(iE 1/2 ).
The discrete modes are spatially localized around the origin and have negative E. In-
troducing the auxiliary parameter ν(») via

» ≡ ν (ν + 1) (7.14)

the discrete eigenvalues are

Ej = ’(ν ’ j)2 , j = 0, 1, ..., jmax (7.15)

where jmax is the largest integer smaller than ν, i. e., the number of allowed modes is one
plus the integer part of ν.
The good news is that both the continuum and discrete eigenmodes can be calculated
by spectral methods. The continuous spectrum requires special tricks; see Sec. 4 of Chapter
19 for an example. The discrete modes can be computed by exactly the same tricks as used
for the quantum harmonic oscillator, whose countable in¬nity of modes are the Hermite
functions. The only complication is that the number of “good” modes is now ¬xed, and
will not increase above jmax + 1 even if we used 10,000 spectral basis functions.
Laplace™s Tidal Equations are even trickier, an eigenproblem of the “Third Kind” in
Boyd™s terminology:

»u ’ xv ’ sζ =0
xu ’ »v + Dζ =0
su ’ Dv ’ » (1 ’ x2 )ζ (7.16)

where x is the coordinate (cosine of colatitude) and D ≡ (1 ’ x2 )d/dx as for Legendre™s
equation, u and v are the horizontal ¬‚uid velocities, and ζ is the sea surface height (in
oceanographic applications) or the pressure (in meteorological use). Depending on the
application, any of the set of three parameters (», s, ) may be the eigenvalue where » is the
nondimensional frequency, s is the zonal wavenumber and is “Lamb™s parameter”, which
is proportional to the depth of the water or (in meteorology) to the vertical wavelength.
The eigenproblem is of the Third Kind because for modes of low frequency, the differ-
ential system is singular at those latitudes (“inertial latitudes” or “critical latitudes”) where
» = x. However, the eigenfunctions themselves are well-behaved and analytic even at the
inertial latitudes.

When the goal is to calculate the free tidal oscillations of a global ocean, s and are
known and the frequency is the eigenvalue. All modes are discrete and this case is rela-
tively uncomplicated.
For forced atmospheric tides, s and the frequency » are known. However, the earliest
calculations agreed badly with observations, even the limited data of the 1930™s. All sorts
of imaginative (but wrong!) theories were proposed, but the correct explanation was given
independently after thirty years of confusion by Lindzen and by Kato. It had been assumed
by analogy with Sturm-Liouville eigenproblems of the First Kind that all the eigenvalues
were positive. Actually, because of the inertial latitudes, the diurnal tide also has an in¬nite
number of modes with negative . One class of modes is oscillatory at low latitudes (be-
tween the inertial latitudes) and decays exponentially beyond the inertial latitudes as the
poles are approached. The other class of modes is con¬ned mostly poleward of the inertial
latitudes and has little amplitude in the tropics. The singularities of the differential equa-
tion are crucial even though the eigenfunctions themselves are paragons of mathematical
A third application (O™Connor, 1995, 1996, Boyd, 1996c) is ocean oscillations in a basin
bounded by meridians where now the eigenparameter is the zonal wavenumber s. Com-
plex eigenvalues are possible, and merely describe oscillations whose amplitude decays
exponentially with increasing distance from the coast. The surprise, as yet unsupported by
rigorous proof, is that the eigenvalues are a mixture of a discrete and a continuous spec-
trum. In addition, unless s is an integer, the eigenfunctions have weak singularities at the
poles. This implies that for the spectral series of an eigenmode (using any standard basis),
the asymptotic rate of convergence is algebraic rather than exponential unless special tricks
such as an exponential change-of-coordinate are used (Chapter 16, Sec. 5).
Lastly, Boyd (1981a, 1982a, 1985a) has studied Sturm-Liouville eigenproblems of the
Fourth Kind such as

uxx + (1/x ’ »)u = 0, (7.17)
u(a) = u(b) = 0

where the interval x ∈ [a, b] spans the origin. Super¬cially, this looks like a self-adjoint SL
problem of standard form, which would be expected to have only real discrete eigenvalues.
In reality, the presence of the pole in the coef¬cient of the undifferentiated term changes
everything. The pole must be interpreted as the limit

1/(x ’ iδ) (7.18)

as δ tends to zero where δ represents viscosity. For non-zero δ, there is a small viscous
layer of thickness proportional to δ. Outside this layer, the eigenfunctions are effectively
inviscid. In the limit δ ’ 0, the thickness of the viscous layer shrinks to zero and the
eigenfunctions are singular as x log(x) at x = 0. The eigenfunctions and eigenvalues are
both complex-valued; one must circle the branchpoint below the real axis to obtain the
branch which is the correct limit of the viscous solution.
The good news is that Boyd (1981a, 1982a) and Gill and Sneddon(1995, 1996) developed
special spectral methods to compute the singular eigenfunctions. The bad news is that
Boyd (1985a) showed that the author™s earlier papers had missed a root.
Clearly, Sturm-Liouville eigenproblems can be full of surprises. Hydrodynamic stabil-
ity problems, which are not self-adjoint and usually have complex-valued eigenvalues and
nearly-singular, complex-valued eigenfunctions, merely reiterate this theme: Eigenprob-
lems can be nasty, and it is terribly easy to be speared by internal layers or singularities, or
miss modes entirely.

7.5 Winnowing the Chaff: Criteria for Rejecting Numeri-
cally Inaccurate or Physically Non-Existent Eigenvalues
It is always a good idea to repeat each calculation twice with different N to verify that
the solution is well-resolved. With eigenproblems, this is doubly important because, as
expressed by the Eigenvalue Rule-of-Thumb, many of the eigenvalues of the discretization
matrix are numerical nonsense unrelated to the eigenvalues of the underlying differential
or integral eigenproblem. To know how many modes are “good”, that is, have been com-
puted accurately, one must compare the list of eigenvalues for two different N and accept
only those which are the same (within some user-set tolerance) on both lists.
One dif¬culty is that if one is performing a large number of eigencalculations, comparing-
by-eye can lead to eyestrain, headaches, and a ¬rm resolve to ¬nd a new, more exciting pro-
fession like accounting or tax law. This is one motive for inventing some simple graphical
and numerical methods for reliably separating “good” eigenvalues from trash.
Another is that we often do not know a priori how many discrete eigenmodes even exist,
as illustrated by the Legendre equation.
Our recommended strategy is to make a plot on a logarithmic scale of the reciprocal of
the difference between corresponding eigenvalues as calculated at different resolutions N ,
scaled by some measure of the size of the eigenvalues. The reciprocal of the difference is
plotted so that the “good” eigenvalues are at the top of the graph; a semilogarithmic plot is
recommended because the accuracy varies exponentially with mode number j as already
seen above. There are two minor complications.

Legendre Eq.: 4 bound states N=90, 120

0 20 40 60 80
mode number j

Figure 7.7: The reciprocal eigenvalue drift ratios 1/δj,nearest (circles) and 1/δj,ordinal (x™s)
are plotted on a logarithmic scale versus mode number j for Legendre™s equation uyy +
{ 15.75 sech2 (y) + E}u = 0 where E is the eigenvalue (the energy). The modes are ordered
by real part of E with j = 1 the smallest. For this case, there are precisely four bound states
which are the four numerical eigenvalues with E < 0.

First, the obvious scaling for the j-th eigenvalue is |»j | itself. However, an eigenvalue
may accidentally be very close to zero. Furthermore, for a “nice” Sturm-Liouville problem,
the eigenvalues typically grow as O(j 2 ), but the difference between adjacent eigenvalues
grows only as fast as O(j). We can solve both these dif¬culties by de¬ning an “intermodal
separation” via

≡ |»1 ’ »2 |
≡ (|»j ’ »j’1 | + |»j+1 ’ »j |) , (7.19)
σj j>1
and then conservatively scale the eigenvalues by σj . (In the case of degeneracy, that is,
two or more eigenfunctions with the same eigenvalue, it is assumed that the differences in
Eq. (7.19) are taken between the non-degenerate modes nearest to mode j.)
Second, in many problems, the eigenvalue ordering is invariant when the resolution
is changed ” at least after the eigenvalues for each run have been sorted by magnitude
or real part ” and it is suf¬cient to plot the reciprocal of what we shall dub the (scaled)
“ordinal” difference

δj,ordinal ≡ |»j (N1 ) ’ »j (N2 )| / σj (7.20)

where the arguments denote the number of degrees of freedom N in the low and high
resolution computations. However, for many problems, such as Laplace™s Tidal Equation
(O™Connor, 1995, 1996, Boyd, 1996c), the well-resolved eigenvalues are not evenly spaced.
(Tidal oscillations fall into two main classes: Rossby waves, which have low frequency, and
gravity waves, which have high frequency. Thus, the gravest Rossby and gravity modes
(well-resolved) are separated by an in¬nite number of higher Rossby modes of intermedi-
ate frequency, which are not as accurately computed.) For such more intricate relationships
between eigenvalue magnitude and the spatial scale of the wave, one needs to compare the
j-th low resolution mode with whatever eigenvalue of the high resolution computation
agrees most closely with it. This “nearest” difference (scaled) is

δj,nearest ≡ min |»j (N1 ) ’ »k (N2 )| / σj (7.21)
k∈[1,N2 ]

Fig. 7.7 is a plot of the scaled differences or “drift-with-N ” for the equation solved by
the Associated Legendre equation using 90 and 120 rational Chebyshev functions on the
in¬nite interval. The exact solution has four discrete modes plus a continuous spectrum.
The discrete modes of the spectral matrix cannot converge to the continuous spectrum, so
that the best we can do is to resolve the four discrete modes. The graph shows that with
N = 90, these modes are very well-resolved in the sense that these eigenvalues, scaled by
the smaller of |»j | or σj , change by less than one part in a million when the basis size is
increased to N = 120. For this case, the “ordinal ratio” would have been suf¬cient be-
cause the discrete modes are the four smallest numerically computed eigenvalues at both
resolutions; the “ordinal” and “nearest” ratios are so close that the circles and x™s are super-
imposed. For the unresolved modes of the continuous spectrum (bottom of the graph), the
“nearest” ratios are consistently a little larger than the “ordinal” ratios. The reason is that
by coincidence, two inaccurate eigenvalues for different resolution may be rather close,
producing a nearest ratio which is small even though the spectral method is generating
random numbers for these modes. It is both quicker and less confusing to use the “or-
dinal” ratio wherever possible. As shown in Boyd(1996c), however, other eigenproblems
such as Laplace™s Tidal equation absolutely require the “nearest” ratio.
Fig. 7.7 is easy to interpret ” four discrete modes, all the rest nonsense ” because the
reciprocal ratios 1/δj are so many orders of magnitude larger for the “good” modes than

for the bad. A low precision ¬nite difference calculation would be much more ambiguous
“ is this eigenvalue real or just a numerical artifact that would jump around at higher

7.6 The Curse of “Spurious” Eigenvalues
Generalized eigenproblems may have some eigenvalues which are physically spurious rather
than merely numerically underresolved. The dif¬culty was ¬rst characterized by Gottlieb
and Orszag (1977, pg. 145), who noted “low modes are given accurately ... but there ap-
pear spurious unstable modes with large growth rates. Similar spurious unstable modes
appear in ¬nite-difference solution of the Orr-Sommerfeld equation.” (pg. 145).
Much time and energy has been expended in the invention of slight modi¬cations
to standard spectral methods to solve the “spurious eigenvalue” dif¬culty: Gottlieb and
Orszag (1977, pg. 143-146), Brenier, Roux and Bontoux(1986), Zebib(1987b), Gardner, Trog-
don and Douglass(1989), McFadden, Murray and Boisvert (1990), Huang and Sloan (1994),
and Dawkins, Dunbar, and Douglass (1998). However, a major theme of the chapter is that
convert-to-matrix methods always have lots of nonsense eigenvalues, so what is so special
about these “spurious eigenvalues”?
The short answer is that because the “spurious eigenvalues” are spurious due to mis-
represented physics, rather than mere underresolution of genuine eigenmodes, there are
some differences between them and the other unacceptable eigenvalues which are worth
discussing. Before we can discuss these differences, however, we must ¬rst note that the
bland label “spurious eigenvalues”, used in many previous studies, is a semantic atrocity
because it blurs the distinction between these eigenvalues and those which are in error
merely because N is too small.
To correct this semantic sloppiness, we offer the following.

PHYSICALLY SPURIOUS EIGENV ALUES are numerically-computed eigenvalues which
are in error because of misapplication of boundary conditions or some other misrepresentation of the
NUMERICALLY SPURIOUS EIGENV ALUES are poor approximations to exact eigenvalues
because the mode is oscillating too rapidly to be resolved by N degrees of freedom. A given numeri-
cally spurious eigenvalue can always be computed accurately by using suf¬ciently large N .

An example, ¬rst studied by Gottlieb and Orszag (1977, pg. 143-145), will help. The
equations for a viscous, incompressible ¬‚uid can, in the limit of small amplitude and one-
dimensional ¬‚ow, be written as the system

= νζxx
ζ = ψxx

where ν is the constant viscosity. The bounary conditions are

ψ(±1) = ψx (±1) = 0

Note that these are entirely on the streamfunction ψ with no boundary conditions on the
vorticity ζ. This system can be reduced to the single equation:

ψtxx = ν ψxxxx

The Galerkin method, called the “tau method” in Gottlieb and Orszag, is to write
ψ(t) ≈ (7.25)
an (t)Tn (x)

and then impose the N ’ 3 conditions that the residual should be orthogonal, in the usual
Chebyshev inner product, to Tk (x) for k = 0, 1, . . . , N ’ 4. These conditions are supple-
mented with the four boundary conditions (7.23).
Unfortunately, the resulting system of ODEs in time is a disaster wrapped in a catas-


. 31
( 136 .)