. 32
( 136 .)


trophe: all time-marching schemes for the system are unstable. Note that both (7.24) and
its Chebyshev discretization can be solved exactly by expansion in the solutions (exact or
Chebyshev-discretized, respectively) of the eigenvalue problem

ν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

ψt = ψxxxx

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-

ν uxxxx = »u

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

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
(1 ’ x2 )
hj (x) ≡ (7.29)
Cj (x)
(1 ’ x2 )

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

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,
dp T N 1
(±1) ∼ N 2p 1 + O(1/N 2 (7.31)
dx 2k + 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.


T2j - T 0 Basis



2 Heinrichs Basis

0 1 2 3
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.



Condition number

T - T Basis
2j 0


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

1 2
10 10

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)

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)
Ck (x) =
(1 ’ x2 )2 πx (xk )(x ’ xk )

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


. 32
( 136 .)