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)

=0

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.

CHAPTER 7. LINEAR EIGENVALUE PROBLEMS

136

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

regularity.

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. CRITERIA FOR REJECTING EIGENVALUES 137

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

14

10

12

10

10

10

8

10

6

10

4

10

2

10

0

10

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.

CHAPTER 7. LINEAR EIGENVALUE PROBLEMS

138

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 |

σ1

1

≡ (|»j ’ »j’1 | + |»j+1 ’ »j |) , (7.19)

σj j>1

2

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

7.6. “SPURIOUS” EIGENVALUES 139

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

resolution?

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.

De¬nition 20 (SPURIOUS EIGENVALUES)

PHYSICALLY SPURIOUS EIGENV ALUES are numerically-computed eigenvalues which

are in error because of misapplication of boundary conditions or some other misrepresentation of the

physics.

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

ζt

(7.22)

ζ = ψxx

where ν is the constant viscosity. The bounary conditions are

(7.23)

ψ(±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:

(7.24)

ψtxx = ν ψxxxx

CHAPTER 7. LINEAR EIGENVALUE PROBLEMS

140

The Galerkin method, called the “tau method” in Gottlieb and Orszag, is to write

N

ψ(t) ≈ (7.25)

an (t)Tn (x)

j=0

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-