its Chebyshev discretization can be solved exactly by expansion in the solutions (exact or

Chebyshev-discretized, respectively) of the eigenvalue problem

(7.26)

νuxxxx = »uxx

where each mode depends upon time as exp(»t) where » is the eigenvalue. The exact

eigenvalues of the original PDE are » = ’νµ where either µ = nπ or µ is any nonzero

root of tan µ = µ; the eigenvalues are negative real so that the exact eigenmodes all de-

cay monotonically with time. The Chebyshev discretized problem, however, has two pos-

tive real eigenvalues whose magnitude increases as roughly O(N 4 ). Thus, the Chebyshev-

discretized system of ODEs in time has two eigenmodes which blow up very fast, and no

time-marching scheme can forestall disaster.

These two physically-spurious eigenvalues arise because (7.24) is an “explicitly-implicit”

problem, that is, in order to apply standard software for solving the system of ordinary dif-

ferential equations in time, we must rewrite the equation as

1

(7.27)

ψt = ψxxxx

2

‚xx

The problem is that we impose four boundary conditions on the streamfunction ψ, consis-

tent with the fact that the differential equation is fourth order, but the differential operator

that must be inverted to compute the right-hand side of (7.27) is only of second order. In

general, a second order differential equation with four boundary conditions is insoluble. A

similar inconsistency arises when converting (7.26) from a generalized to a standard eigen-

problem:

1

(7.28)

ν uxxxx = »u

2

‚xx

This suggests “ as McFadden, Murray and Boisvert (1990) have con¬rmed through

many numerical experiments “ that these physically spurious eigenvalues only arise for

generalized eigenproblems, and are absent for standard eigenproblems which do not in-

volve the inversion of a differential operator. Furthermore, this dif¬culty happens only for

non-periodic problems. In a periodic problem, one cannot have inconsistencies in the num-

ber of boundary conditions because the boundary condition of periodicity is the same for

differential operators of all orders. It is only generalized eigenproblems in non-periodic

geometry where physically spurious eigenvalues may arise.

Dawkins, Dunbar and Douglass (1998) have rigorously proved the existence of the large

positive eigenvalues for the eigenproblem (7.26). They show that the overspeci¬cation of

boundary conditions technically creates two L2 eigenmodes of the differential eigenprob-

lem with in¬nite eigenvalues, which the discretization approximates as positive eigenval-

ues of magnitude O(N 4 ).

7.6. “SPURIOUS” EIGENVALUES 141

Several remedies have been proposed, but all seem to be closely related, so we shall

describe only the variant due to Huang and Sloan (1994). Their spectral basis is composed

of Chebyshev-Lobatto cardinal functions on an N -point grid xk , k = 1, . . . , N . Dirichlet

homogeneous boundary conditions are imposed by omitting C1 (x) and CN (x) from the

basis. The twist is that for (7.26), a non-standard basis is used to represent the fourth

derivatives:

(1 ’ x2 )

hj (x) ≡ (7.29)

Cj (x)

(1 ’ x2 )

j

where the Cj (x) are the standard cardinal functions. Note that the factor of (1’x2 ) enforces

the homogeneous Neuman condition so that a sum of hj from j = 2, . . . , N ’ 1 must have

a double zero at both endpoints. The coef¬cients of both basis sets are the same: the values

of u(x) at the interior grid points. However, the modi¬ed functions hk are used only to

represent the fourth derivative while the standard basis functions with “lazy” imposition

of boundary conditions (just two) are used for the second derivatives. Thus, the discrete

generalized eigenvalue problem Au = » Bu for Eq.(7.26) has the matrix elements

(7.30)

Aij = hj,xxxx (xi ), Bij = Cj,xx (xi )

where u is the column vector of grid point values {u(x2 ), . . . , u(xN ’1 )}.

This works, but the question still remains: With so many numerically underresolved

eigenvalues to battle, do a couple of physically spurious eigenvalues matter? We shall

offer three positive answers, each followed by a rejoinder.

The ¬rst argument is that these physically spurious eigenvalues are very important

because they cause unconditional instability of time-marching schemes. However, this is

really an issue to be fought out in a discussion of time-marching. It is not a “spurious

eigenvalues” problem so much as it is a “screwing up the bounary conditions” problem.

Further, Gottlieb and Orszag note that “this version of the tau method [with physically

spurious eigenvalues] may be suitable for eigenvalue problems even though it is uncondi-

tionally unstable for the initial-value problem”.

Spurious modes are often generated when a system of equations is reduced to a single

equation, as often happens in setting up eigenproblems as illustrated by Eq.(7.26). The

most famous example of “spurious” modes is the generation of unreal pressure modes

when the equations of incompressible ¬‚uid ¬‚ow are reduced to a single equation for the

pressure. Although there are no physical boundary conditions for the pressure, the Poisson-

like equation for it requires boundary conditions “ and evaluating the Navier-Stokes equa-

tions at solid walls actually gives more boundary conditions than are needed.

The remedy is now well-understood: discretize the system, impose boundary conditions

on the system, and then reduce the system to a single equation after discretization. Another

important tactic is to use polynomials of different degrees where necessary. Again, the most

well-known illustration is the practice of approximating the pressure in incompressible

¬‚ow by a polynomial whose degree is two less than that of the velocity. Similarly, Huang

and Sloan™s strategy approximates the second derivative using basis functions whose de-

gree is two less than that of the modi¬ed cardinal functions hj which are used for the fourth

derivative.

So, physically spurious eigenvalues do have some signi¬cance because they imply tem-

poral instability. However, this is a well-understood phenomena now, and besides it is a

tale that belongs in another chapter. When the goal is only eigenvalues, temporal instabil-

ity does not matter.

The second argument in favor of modi¬ed eigenvalue-solving is that if the target is

unstable modes, then physically spurious modes of large positive real part might be con-

fused with genuine instabilities, and so waste a lot of time. However, as explained in the

CHAPTER 7. LINEAR EIGENVALUE PROBLEMS

142

previous section, this shouldn™t be a dif¬culty if one plots the change in eigenvalues with

resolution. The physically-spurious eigenvalues seem to scale as O(N 4 ) for most problems

and thus are among the very easiest “bad” eigenvalues to detect through sensitivity to N .

The third argument for purging physically spurious eigenvalues is that misrepresen-

tation of boundary conditions worsens the condition number of the matrices. This is per-

fectly true, and is particularly annoying since eigenproblems “ again because of the reduc-

tion of a system of equations to one differential equation “ often involve high order deriva-

tives. Since the derivatives of Chebyshev polynomials oscillate more and more wildly in

narrow boundary layers as the order of the derivative increases, Chebyshev differentiation

matrices become more and more ill-conditioned for problems with high order derivatives,

with or without physically spurious eigenvalues. We shall describe how to ameliorate this

ill-conditioning in the next section.

7.7 Reducing the Condition Number

One unpleasant dif¬culty with a Chebyshev or Legendre basis is that the derivatives of

the polynomials oscillate near the endpoints with increasing amplitude and decreasing

wavelength as the degree and order of differentiation increase. To be precise,

p

dp T N 1

(±1) ∼ N 2p 1 + O(1/N 2 (7.31)

p

dx 2k + 1

k=1

In turn, this implies that the condition number of the Chebyshev and Legendre matrices for

a differential equation are O(N 2p ). Thus, for a sixth order equation, the matrix is blowing

up as O(N 12 ), and the linear algebra subroutine is likely to crash spectacularly even for N

little larger than 10.

6

10

5

T2j - T 0 Basis

10

4

10

3

10

2 Heinrichs Basis

10

1

10

0 1 2 3

t=arccos(x)

Figure 7.8: Absolute values of the second derivative of the highest basis function versus

the trigonometric argument t = arccos(x) for two different basis sets (N = 50 collocation

points). Upper curve: Each basis function is the difference of a Chebyshev polynomial

with either 1 or x. Its maximum is 6.9E5. Bottom curve: Heinrichs™ basis: φ50 (x) = (1 ’

x2 ) T49 (x). The maximum is only 2.4E3.

7.7. REDUCING THE CONDITION NUMBER 143

7

10

6

10

5

Condition number

T - T Basis

10

2j 0

4

10

3

10

2

(1-x*x)*T j(x)

10

1

10

1 2

10 10

N

Figure 7.9: Solid: condition numbers for the matrices which discretize the second deriva-

tive operator in two different basis sets, plotted versus N , the size of the matrix. Top: the

“difference” basis in which each function is Tj (x) minus either 1 or x, depending on parity.

Bottom with circles: Heinrichs™ basis: φj (x) = (1 ’ x2 ) Tj (x). The condition number is

de¬ned to be the ratio of the largest to the smallest singular value of the matrix. The dotted

lines have slopes of N 4 (top) and N 2 (bottom), illustrating that the rate of growth of the

condition number closely follows these respective power laws.

Fortunately, there are a couple of remedies. Wilhelm Heinrichs (1989a, 1991b, 1991c)

observed that if we apply basis recombination to create basis functions that satisfy homo-

geneous boundary conditions, it is possible to choose particular sets that greatly reduce the

condition number. For example, suppose the target is to solve a second order differential

equation with the homogeneous Dirichlet boundary conditions u(±1) = 0. The simplest

choice of basis functions, and one that works very well except for N >> 100, is the “differ-

ence” basis

φ2j (x) ≡ T2j (x) ’ T0 , φ2j’1 (x) ≡ T2j’1 (x) ’ T1 (7.32)

but the second derivative of these is equal to that of a Chebyshev polynomial and thus this

basis gives poor condition number.

Heinrichs proposed instead

φj ≡ (1 ’ x2 )Tj , (7.33)

j = 0, 1, . . .

The endpoint and near-the-endpoint values of the second derivatives of these functions,

which we shall dub the “Heinrichs” basis, are much smaller than those of the Chebyshev

polynomials or the difference basis: O(N 2 ) basis versus O(N 4 ). His reasoning is that the

second derivative of his basis functions is

φj,xx = (1 ’ x2 )Tj,xx ’ 4xTj,x ’ 2Tj (7.34)

CHAPTER 7. LINEAR EIGENVALUE PROBLEMS

144

Now Tj,xx is O(N 4 ) near the endpoints, but this is precisely where the factor (1 ’ x2 ) tends

to zero. The contribution of the ¬rst derivative is not cancelled near the boundary, but since

the ¬rst derivative is only O(N 2 ), the second derivative of the Heinrichs basis is no larger.

Fig. 7.8 compares the second derivatives of two basis functions of the same degree in

the two different basis sets. It is necessary to plot the absolute values on a logarithmic

scale, and to plot them versus the trigonometric coordinate t (= arccos(x)) because the

oscillations grow so rapidly in magnitude and wavelength near the endpoints that the

obvious plot of the second derivatives on a linear scale versus x shows little but some

wiggly, indecipherable spikes near x = ±1. The lower plot shows that the (1 ’ x2 ) factor in

the Heinrichs basis turns the magnitude back towards zero so that the maximum amplitude

is roughly the square root of the maximum for the “difference” basis. Fig. 7.9 shows that this

translates into a great difference in the condition number of the Chebyshev discretization

matrices for the second derivative: O(N 2 ) for the Heinrichs basis versus O(N 4 ) (ouch!) for

the “difference” basis.

There are a couple of little white lies hidden in these graphs. First, the Heinrichs™ basis

makes the condition number of the matrix which discretizes the undifferentiated basis worse.

The result is that when we convert uxx + »u = 0 into the standard matrix eigenproblem

Gu = » u with G = B’1 A where A is the discretization of the second derivative and B that

of the undifferentiated term, it turns out the condition number of G is almost identical for

the two basis sets. The other white lie is that for a second order problem, ill-conditioning

is a problem only when N > 100, which is to say it is not a problem.

For fourth order and higher order differential equations, however, ill-conditioning is

more of a problem, and the Heinrichs™ basis comes increasingly into its own. The idea

generalizes in an obvious way to higher order equations. For example, for a fourth order

problem with the boundary conditions u(±1) = ux (±1) = 0, the Heinrichs basis is

φj (x) ≡ (1 ’ x2 )2 Tj (x) (7.35)

where the quartic multiplier enforces a double zero at both endpoints and thus lowers the

condition number of the fourth derivative from O(N 8 ) to O(N 4 ).

Huang and Sloan (1994) propose a related strategy which is more suitable for a cardinal

function basis. For a fourth order problem, their cardinal functions are

(1 ’ x2 )2 π(x)

(7.36)

Ck (x) =

(1 ’ x2 )2 πx (xk )(x ’ xk )

k

where π(x) is the product of monomial factors constructed from the interior points of the