Continuum Mechanics and Thermodynamics manuscript No. (will be inserted by the editor)

Hydrodynamic Stability, the Chebyshev Tau Method and Spurious Eigenvalues

David Bourne

Department of Mathematical Sciences, University of Durham, South Road, Durham, DH1 3LE, U.K.

e-mail: David.Bourne@dunelm.org.uk

July 4, 2002

Key words Bénard convection, Orr-Sommerfeld equation, Chebyshev tau method, Spurious eigenvalues.

Abstract The Chebyshev tau method is examined; a numerical technique which in recent years has been successfully applied to many hydrodynamic stability problems. The orthogonality of Chebyshev functions is used to rewrite the differential equations as a generalized eigenvalue problem. Although a very efficient technique, the occurrence of spurious eigenvalues, which are not always easy to identify, may lead one to believe that a system is unstable when it is not. Thus, the elimination of spurious eigenvalues is of great importance. Boundary conditions are included as rows in the matrices of the generalized eigenvalue problem and these have been observed to be one cause of spurious eigenvalues. Removing boundary condition rows can be difficult. This problem is addressed here, in application to the Bénard convection problem, and to the Orr-Sommerfeld equation which describes parallel flow. The procedure given here can be applied to a wide range of hydrodynamic stability problems.

1 Introduction

The Chebyshev tau method is an efficient numerical technique for finding eigenvalues and eigenfunctions of systems in continuum mechanics. One strength of the algorithm is that it calculates as many eigenvalues as required, not just the dominant one. Dongarra et al. (1996) highlight work done in viscous shear flow and viscous channel flow, by Butler and Farrell (1992) and Reddy and Henningson (1993) respectively, to show the importance of examining more than just the leading mode. The main disadvantage of the Chebyshev tau method is the occurrence of spurious eigenvalues.

These are values returned by the algorithm which do not satisfy the eigenvalue problem.

The Chebyshev tau method and the spurious eigenvalue problem have been the attention of much study recently. Gardner et al. (1989) and Mc-Fadden et al. (1990) describe modified (Chebyshev) tau methods to avoid spurious eigenvalues. Straughan (2001) observes the problems in detecting spurious eigenvalues. They need not be very large, e.g. of O(1), and the number of them often changes with the number of Chebyshev polynomials used. Gardner et al. (1989) calculate the tau coefficients to search for spurious eigenvalues. An alternative approach, as discussed in Dongarra et al 1996, is to look at the eigenfunctions and residuals. For a fascinating insight into the origin of spurious eigenvalues see Dawkins et al 1998.

Dongarra et al.(1996), Payne and Straughan (2000) and Straughan (2001) discuss the merits of the D^{4}, D^{2} and D versions of the Chebyshev tau method (in which the systems of equations are written as order four, order two and order one, respectively). The matrix coefficients in the D^{4} method have large growth rates which can lead to serious round off errors. The D^{2} method has a smaller growth rate and the D method is even better. However, the D method requires matrices of twice the width of those used in the D^{2} method and so computing times are longer (although this is not such a big problem due to modern computer speeds). The choice between the D and D^{2} methods seems to depend mostly on the boundary conditions and on which Chebyshev matrices have to be calculated.

McFadden et al. (1990) noticed that rows of zeros, correponding to the boundary conditions, in the generalized eigenvalue problem are one source of spurious eigenvalues. For many problems it is necessary to use the D or D^{2} method to reduce the order of the equations. This introduces new variables for which we have no boundary conditions. Thus it is difficult to remove zero boundary condition rows (in an attempt to suppress spurious eigenvalues). This problem was first addressed by Straughan (2001), for porous convection situations. We present a method which we believe to be very well suited for programming into a computer.

In this note we illustrate our ideas by looking at Bénard convection and the Orr-Sommerfeld equation, in particular plane Poiseuille flow and plane Couette flow.

2 Bénard Convection

2.1 The Chebyshev tau method

When a horizontal layer of fluid is heated from below the fluid near the bottom of the layer expands and, under the right conditions, will rise forming a regular pattern of convection cells. It turns out that for this instability to occur the temperature gradient or layer depth must exceed some critical value. This is known as Bénard convection.

One can show that the following equations govern the linear instability of the system (see Chandrasekhar, 1981, Chap. II):

\begin{aligned}(D^{2} − a^{2})^{2}W & = & Ra^{2}Θ \\ (D^{2} − a^{2})Θ & = & −RW, \end{aligned}  (1)

for z ∈ (0, 1). D is the differential operator d/dz and a is the wavenumber. Θ(z) and W (z) are obtained by separating off the z dependence from the temperature perturbation and the vertical component of the velocity perturbation, respectively. R^{2} is the Rayleigh number and, for a given fluid, is proportional to the temperature gradient multiplied by the cube of the layer depth.

System (1) is of order six and so we require six boundary conditions. The following four conditions are frequently enforced:

W (0) = W (1) = Θ(0) = Θ(1) = 0.  (2)

The final two boundary conditions depend on what type of surfaces bound the fluid layer. The three cases usually considered are

(i) free-free: bottom and top surfaces free,

D^{2}W (0) = D^{2}W (1) = 0,  (3)

(ii) rigid-rigid: bottom and top surfaces fixed,

DW (0) = DW (1) = 0,  (4)

(iii) rigid-free: bottom surface fixed, top surface free,

DW (0) = D^{2}W (1) = 0.  (5)

See Chandrasekhar (1981 pp.21-22) for a derivation of these boundary conditions. The free surface condition assumes that the interface remains flat and ignores surface tension.

System (1) is an eigenvalue problem for R, for a fixed. Let R_{crit} be the least positive eigenvalue of this system. We have linear stability if R < R_{crit} and linear instability if R > R_{crit}. As an interesting aside, we mention that the nonlinear stability boundary also satisfies (1), Straughan (1992, pp.51-55).

We will implement the D^{2} version of the Chebyshev tau method to solve this eigenvalue problem. Let

A(z) = (D^{2} − a^{2})W  (6)

Using this, we can rewrite system (1) as

\begin{aligned}L_{1}(W, A, Θ)^{T} & ≡ & (D^{2} − a^{2})W − A & = & 0, \\ L_{2}(W, A, Θ)^{T} & ≡ & (D^{2} − a^{2})A − Ra^{2}Θ & = & 0, \\ L_{3}(W, A, Θ)^{T} & ≡ & (D^{2} − a^{2})Θ + RW & = & 0. \end{aligned}  (7)

Note that, following the D^{2} method, we have changed from a system of two equations, of order four and order two, to a system of three order two equations.

Transform the equations in system (7) to the spatial domain (−1, 1), decompose W, A and Θ as an infinite sum of Chebyshev polynomials and then trancate each sum as follows:

W = ∑_{i=0}^{N+2} w_{i}T_{i}(z), A = ∑_{i=0}^{N+2} a_{i}T_{i}(z), Θ = ∑_{i=0}^{N+2} θ_{i}T_{i}(z).  (8)

The Chebyshev tau method asserts that, because of the truncation, instead of solving system (7) we solve

L_{1}(W, A, Θ)^{T} = τ_{11}T_{N+1} + τ_{12}T_{N+2},
L_{2}(W, A, Θ)^{T} = τ_{21}T_{N+1} + τ_{22}T_{N+2},  (9)
L_{3}(W, A, Θ)^{T} = τ_{31}T_{N+1} + τ_{32}T_{N+2}

where τ_{11}, …, τ_{32} are the tau coefficients.

We have 3(N + 3) = 3(N + 1) + 6 coefficients to determine: w_{i}, a_{i} and θ_{i} for i = 0, …, N + 2. We take the inner product of T_{i} with each equation in system (9), for i = 0, …, N. This yields 3(N + 1) equations governing our coefficients. The remaining six equations that we need come from the boundary conditions. By doing this, we obtain the generalized eigenvalue problem

Ax = R Bx  (10)

where

x = (w_{0}, …, w_{N+2}, a_{0}, …, a_{N+2}, θ_{0}, …, θ_{N+2})^{T}  (11)

and

A = ( \begin{matrix}D^{2} − a^{2}I & −I & 0 \\ 0 & D^{2} − a^{2}I & 0 \\ 0 & 0 & D^{2} − a^{2}I \\ boundary & condition & rows \end{matrix} )  (12)

B = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & a^{2}I \\ −I & 0 & 0 \\ boundary & condition & rows \end{matrix} )_{.}  (13)

A and B are matrices of dimension 3(N + 3) × 3(N + 3). The matrix blocks, such as D^{2} −a^{2}I, have dimension (N +1)×(N +3). I is the identity matrix and there are six boundary condition rows. In Dongarra et al. (1996) the elements of the Chebyshev tau differentiation matrices D and D^{2} (which represent d/dz and d^{2}/dz^{2} respectively) are given as

\begin{aligned}D_{0,2j−1} & = & 2j − 1, & j ≥ 1, \\ D_{i,i+2j−1} & = & 2(i + 2j − 1), & i ≥ 1, j ≥ 1, \\ D_{0,2j}^{2} & = & \frac{1}{2}(2j)^{3}, & j ≥ 1, \\ D_{i,i+2j}^{2} & = & (i + 2j)4j(i + j), & i ≥ 1, j ≥ 1. \end{aligned}  (14)

Note that D^{2} = D.D, where . denotes the usual matrix multiplication. The generalized eigenvalue problem, equation (10), can now be solved using the QZ algorithm of Moler and Stewart (1973).

The six boundary condition rows in B are all zero. When using the Chebyshev tau method to analyse stability problems in hydodynamics, zero boundary condition rows, as we have here, often lead to the occurrence of spurious eigenvalues. To try to avoid this problem one would like to use the boundary conditions to rewrite w_{N+1}, w_{N+2}, a_{N+1}, a_{N+2}, θ_{N+1} and θ_{N+2} in terms of w_{0}, …, w_{N}, a_{0}, …, a_{N}, θ_{0}, …, θ_{N}. The zero rows can then be removed along with the columns corresponding to w_{N+1}, w_{N+2}, a_{N+1}, a_{N+2}, θ_{N+1} and θ_{N+2}. The difficulty, which we resolve, is how to find expressions for a_{N+1} and a_{N+2} when all the boundary conditions are on W and Θ.

The expressions for w_{N+1}, w_{N+2}, θ_{N+1} and θ_{N+2} are common to all three sets of boundary conditions. Using the boundary conditions given in equation (2) and the fact that T_{n}(1) = 1 and

T_{n}(0) = \{ \begin{matrix}0 & n odd, \\ 1 & n even, 4 | n, \\ −1 & n even, 4 ∤ n, \end{matrix}  (15)

one easily sees that

w_{N+1} = (−1, 0, 1, 0, …, −1, 0, 1, 0) w,  (16)
w_{N+2} = (0, −1, −2, −1, …, 0, −1, −2, −1) w,  (17)
θ_{N+1} = (−1, 0, 1, 0, …, −1, 0, 1, 0) θ,  (18)
θ_{N+2} = (0, −1, −2, −1, …, 0, −1, −2, −1) θ,  (19)

where w := (w_{0}, …, w_{N})^{T}, θ := (θ_{0}, …, θ_{N})^{T} and we have taken N to be odd and N + 1 divisible by four.

Now we would like to find expressions for a_{N+1} and a_{N+2}. We examine each boundary value problem in turn. Note that there is no need to solve the free-free problem numerically; it was solved analytically in 1916 by Lord Rayleigh (Chandrasekhar, 1981, pp.35-36). However, we include this case for completeness.

2.2 The free-free problem

The free-free boundary conditions are given in equations (2) and (3). In this case there is no challenge in finding a_{N+1} and a_{N+2}. From equation (6) we see that

A(0) = D^{2}W (0) − a^{2}W (0) = 0,  (20)

using our boundary conditions. Similarly A(1) = 0 and so, in exactly the same way as above, we obtain

a_{N+1} = (−1, 0, 1, 0, …, −1, 0, 1, 0) a  (21)
a_{N+2} = (0, −1, −2, −1, …, 0, −1, −2, −1) a  (22)

where a := (a_{0}, …, a_{N})^{T}.

2.3 The rigid-rigid problem

The rigid-rigid boundary conditions are given in equations (2) and (4). Differentiating equation (6) with respect to z we obtain

DA = D^{3}W − a^{2}DW.  (23)

Therefore

\begin{aligned}DA(0) & = & D^{3}W (0), \\ DA(1) & = & D^{3}W (1). \end{aligned}  (24)

Now

DA = ∑_{s=0}^{N+2} a_{s}T′_{s}(z)

= ∑_{s=0}^{N+2} ( ∑_{r=0}^{N+2} D_{rs}T_{r}(z)) a_{s}.  (25)

Therefore

DA(1) = ∑_{s=0}^{N+2} ( ∑_{r=0}^{N+2} D_{rs}T_{r}(1)) a_{s}

= ∑_{s=0}^{N+2} ( ∑_{r=0}^{N+2} D_{rs}) a_{s}
= (D_{1}, …, D_{N+3})(a_{0}, …, a_{N+2})^{T}  (26)

where D_{i} is the sum of the elements in the ith column of D. Let A_{1} = (D_{1}, …, D_{N+3}).

DA(0) = ∑_{s=0}^{N+2} ( ∑_{r=0}^{N+2} D_{rs}T_{r}(0)) a_{s}

= ∑_{s=0}^{N+2} ( ∑_{r=0}^{N+2} D′_{rs}) a_{s}  (27)

where

D′ = diag(1, 0, −1, 0, …, 1, 0, −1, 0) × D.  (28)

Thus

DA(0) = (D′_{1}, …, D′_{N+3})(a_{0}, …, a_{N+2})^{T}  (29)

where D′_{i} is the sum of the elements in the ith column of D′. Let A_{0} = (D′_{1}, …, D′_{N+3}).

In a similar manner we find that

D^{3}W (1) = W_{1}^{3}(w_{0}, …, w_{N})^{T}  (30)

where

W_{1}^{3} = (D_{3,1}, …, D_{3,N+1})  (31)

and D_{3,i} is the sum of the elements in the ith column of D^{3} after the columns relating to w_{N+1} and w_{N+2} have been removed. These columns are removed using the expressions found in the previous section for w_{N+1} and w_{N+2} in terms of w_{0}, …, w_{N}.

Also

D^{3}W (0) = W_{0}^{3}(w_{0}, …, w_{N})^{T}  (32)

where

W_{0}^{3} = (D′_{3,1}, …, D′_{3,N+1})  (33)

and D′_{3,i} is the sum of the elements in the ith column of the matrix

diag(1, 0, −1, 0, …, 1, 0, −1, 0) × D^{3}.

(by D^{3} we really mean D^{3} with the w_{N+1} and w_{N+2} columns removed). Let w = (w_{0}, …, w_{N})^{T} and a = (a_{0}, …, a_{N})^{T}. Our aim is to solve the following equations for a_{N+1} and a_{N+2} in terms of a_{0}, …, a_{N}, w_{0}, …, w_{N}:

A_{0}(a^{T}, a_{N+1}, a_{N+2})^{T} = W_{0}^{3} w
A_{1}(a^{T}, a_{N+1}, a_{N+2})^{T} = W_{1}^{3} w (34)

Rewrite this system as

Ã_{0} a + D′_{N+2}a_{N+1} + D′_{N+3}a_{N+2} = W_{0}^{3} w
Ã_{1} a + D_{N+2}a_{N+1} + D_{N+3}a_{N+2} = W_{1}^{3} w (35)

where

Ã_{0} = (D′_{1}, …, D′_{N+1}) and Ã_{1} = (D_{1}, …, D_{N+1}).  (36)

Rearranging gives

( \begin{matrix}D′_{N+2} & D′_{N+3} \\ D_{N+2} & D_{N+3} \end{matrix} ) ( \begin{matrix}a_{N+1} \\ a_{N+2} \end{matrix} ) = ( \begin{matrix}W_{0}^{3} & −Ã_{0} \\ W_{1}^{3} & −Ã_{1} \end{matrix} ) ( \begin{matrix}w \\ a \end{matrix} )_{.}  (37)

The matrix on the left-hand side is invertible and so we have successfully written a_{N+1} and a_{N+2} in terms of a_{0}, …, a_{N}, w_{0}, …, w_{N}.

2.4 The rigid-free problem

The rigid-free boundary conditions are given in equations (2) and (5). After all the work done in the previous section it is an easy exercise to check that a_{N+1} and a_{N+2} are given by:

( \begin{matrix}D′_{N+2} & D′_{N+3} \\ 1 & 1 \end{matrix} ) ( \begin{matrix}a_{N+1} \\ a_{N+2} \end{matrix} ) = ( \begin{matrix}W_{0}^{3} & −Ã_{0} \\ 0 & −1 \end{matrix} ) ( \begin{matrix}w \\ a \end{matrix} )_{.}

3 The Orr-Sommerfeld Stability Equation

The Orr-Sommerfeld equation, given below, describes the linear stability of viscous parallel flow, i.e. the transition from laminar flow to instability.

(D^{2} − a^{2})^{2}Φ − iaRe(U − c)(D^{2} − a^{2})Φ + iaReU′′Φ = 0  (38)

for z ∈ (−1, 1), subject to boundary conditions

Φ(±1) = DΦ(±1) = 0.  (39)

See Drazin and Reid (1981, Sect. 25). D = d/dz, Re is the Reynolds number, a the wavenumber and c the growth rate. Φ(z) is the amplitude of the stream function Ψ:

Ψ = Φ(z)e^{ia(x−ct)}.  (40)

U(z) is the steady base flow whose stability we are interested in. Two special cases are particularly important:

Plane Poiseuille flow: U(z) = 1 − z^{2}, flow forced by a pressure gradient in the x-direction,

Plane Couette flow: U(z) = z, flow forced by the upper boundary being sheared against the lower one.

The stream function has time dependence exp(−iact) whose modulus is exp(ac_{i}t), in the obvious notation. Therefore we have linear stability when c_{i} < 0 and linear instability when c_{i} > 0.

Once again, we choose to use the D^{2} variety of the Chebyshev tau method to find marginal stability, as recommended for the Orr-Sommerfeld equation by Dongarra et al. (1996). Orszag (1971) applied the Chebyshev tau method to the Orr-Sommerfeld problem to obtain accurate values of the spectrum. He used a D^{4} method and inserted boundary condition rows into the A matrix.

Introduce

A(z) = (D^{2} − a^{2})Φ  (41)

and use this to rewrite equation (38) as the system

\begin{aligned}(D^{2} − a^{2})Φ − A & = & 0, \\ (D^{2} − a^{2})A − iaRe(U − c)A + iaReU′′Φ & = & 0. \end{aligned}  (42)

Expand Φ and A as

Φ = ∑_{i=0}^{N+2} ϕ_{i}T_{i}(z), A = ∑_{i=0}^{N+2} a_{i}T_{i}(z).  (43)

Substitute for Φ and A from (43) into (42) and take the inner product of T_{i} with each equation in system (42), for i = 0, …, N. This yields 2(N + 1) equations for the 2(N + 1) + 4 coefficients ϕ_{0}, …, ϕ_{N+2}, a_{0}, …, a_{N+2}. The final four equations come from the boundary conditions, equation (39). In this way we obtain the generalized eigenvalue problem

Ax = c Bx  (44)

where

x = (ϕ_{0}, …, ϕ_{N+2}, a_{0}, …, a_{N+2})^{T}  (45)

and

A = ( \begin{matrix}D^{2} − a^{2}I & −I \\ ia Re U′′ & D^{2} − a^{2}I − ia Re U \\ boundary & condition rows \end{matrix} )  (46)

B = ( \begin{matrix}0 & 0 \\ 0 & −iaReI \\ boundary & condition rows \end{matrix} )_{.}  (47)

A and B are matrices of dimension 2(N +3)×2(N +3). The matrix blocks have dimension (N + 1) × (N + 3) and there are four boundary condition rows. U and U′′ are the Chebyshev matrices representing U(z) and U′′(z). We examine the two special cases mentioned earlier:

Plane Poiseuille flow: U(z) = 1 − z^{2}

In this case U = I − P, where P is the Chebyshev matrix of z^{2}. In Dongarra et al. (1996) it is suggested that one constructs P by writing z^{2} = \frac{1}{2}(1 + T_{2}(z)), taking the inner product < T_{i}, z^{2}A > and using the fact that

2 T_{m}(z) T_{n}(z) = T_{m+n}(z) + T_{|m−n|}(z).  (48)

Alternatively, and requiring no more or no less effort, one can take the inner product < T_{i}, z^{2}A > and use the following (derived from the recurrence relation for Chebyshev polynomials):

4 z^{2} T_{n}(z) = T_{n−2}(z) + 2T_{n}(z) + T_{n+2}(z).  (49)

Note that U′′ = 2I.

Plane Couette flow: U(z) = z

In this case U = C, where C is the Chebyshev matrix of z. This can be obtained by taking the inner product < T_{i}, zA > and using the recurrence relation

2 zT_{n}(z) = T_{n−1}(z) + T_{n+1}(z).  (50)

U′′ is the (N + 1) × (N + 3) zero matrix.

The four boundary condition rows in B are zero and are thought to be one cause of spurious eigenvalues. Use of the standard Chebyshev tau method to solve the plane Poiseuille problem gives two spurious eigenvalues with positive imaginary part. See McFadden et al. (1989). As in Sect. 2, to avoid this problem we would like to find expressions for ϕ_{N+1}, ϕ_{N+2}, a_{N+1} and a_{N_{2}} in terms of ϕ_{0}, …, ϕ_{N}, a_{0}, …, a_{N}. The zero rows can then be removed along with the columns corresponding to ϕ_{N+1}, ϕ_{N+2}, a_{N+1} and a_{N+2}. The challenge here, as observed by McFadden et al. (1990), is to find expressions for a_{N+1} and a_{N+2} when all the boundary conditions are on Φ and its derivative.

If we take N to be odd and use the fact that T_{n}(±1) = (±1)^{n}, it is easy to show from the boundary conditions Φ(±1) = 0 that

ϕ_{N+1} = (−1, 0, −1, 0, …, −1, 0)ϕ,  (51)
ϕ_{N+2} = (0, −1, 0, −1, …, 0, −1)ϕ,  (52)

where ϕ = (ϕ_{0}, …, ϕ_{N})^{T}.

Differentiate equation (41) with respect to z, evaluate at z = ±1 and use the boundary conditions DΦ(±1) = 0 to find that

\begin{aligned}DA(1) & = & D^{3}Φ(1), \\ DA(−1) & = & D^{3}Φ(−1). \end{aligned}  (53)

As in Sect. 2, one can show that:

DA(1) = (D_{1}, …, D_{N+3})(a_{0}, …, a_{N+2})^{T} =: A_{1}(a_{0}, …, a_{N+2})^{T}, (54)

where D_{i} is the sum of elements in the ith column of D;

DA(−1) = (D′_{1}, …, D′_{N+3})(a_{0}, …, a_{N+2})^{T} =: A_{−1}(a_{0}, …, a_{N+2})^{T},

(55)

where D′_{i} is the sum of the elements in the ith column of the matrix diag(1, −1, …, 1, −1) × D;

D^{3}Φ(1) = (D_{3,1}, …, D_{3,N+1})(ϕ_{0}, …, ϕ_{N})^{T} =: Φ_{1}^{3}(ϕ_{0}, …, ϕ_{N})^{T}, (56)

where D_{3,i} is the sum of the elements in the ith column of D^{3} after the columns relating to ϕ_{N+1} and ϕ_{N+2} have been removed (using expressions 51 and 52);

D^{3}Φ(−1) = (D′_{3,1}, …, D′_{3,N+1})(ϕ_{0}, …, ϕ_{N})^{T} =: Φ_{−1}^{3}(ϕ_{0}, …, ϕ_{N})^{T},

(57)

where D′_{3,i} is the sum of the elements in the ith column of the matrix diag(1, −1, …, 1, −1) × D^{3} (by D^{3} we really mean D^{3} with the ϕ_{N+1} and ϕ_{N+2} columns removed).

Setting

a = (a_{0}, …, a_{N})^{T},  (58)
Ã_{1} = (D_{1}, …, D_{N+1}),  (59)
Ã_{−1} = (D′_{1}, …, D′_{N+1})  (60)

one can rewrite system (53) in the following way:

( \begin{matrix}D_{N+2} & D_{N+3} \\ D′_{N+2} & D′_{N+3} \end{matrix} ) ( \begin{matrix}a_{N+1} \\ a_{N+2} \end{matrix} ) = ( \begin{matrix}Φ_{1}^{3} & −Ã_{1} \\ Φ_{−1}^{3} & −Ã_{−1} \end{matrix} ) ( \begin{matrix}ϕ \\ a \end{matrix} )_{.}

We can solve this for a_{N+1} and a_{N+2} and so eliminate the zero boundary condition rows in B, as desired.

4 Conclusions

In this paper we have described a new method to remove boundary condition rows from Chebyshev matrices. Zero boundary condition rows are certainly not the only source of spurious eigenvalues, see Dawkins et al. (1998). Also, Straughan (2000) found that removing such rows in an inclined temperature gradient convection problem did not eliminate the spurious eigenvalue. However, there are still many practical problems where removing the boundary condition rows is fruitful, e.g. brine motion in thawing subsea permafrost, Staughan (2000), or the stress-free Orr-Sommerfeld problem, Dongarra et al. (1996). The technique we have presented may be applied to many more stability problems in hydrodymaics.

Acknowledgements The author would like to thank Prof. B. Straughan for helpful suggestions.

References

1. Butler KM, Farrell BF (1992) Three-dimensional optimal perturbations in viscous shear flow. Phys. Fluids A 4: 1637-1650

2. Chandrasekhar S (1981) Hydrodynamic and hydromagnetic stability. Dover, New York

3. Dawkins PT, Dunbar SR, Douglass RW (1998) The origin and nature of spurious eigenvalues in the spectral tau method. J. Comput. Phys. 147: 441-462

4. Dongarra JJ, Straughan B, Walker DW (1996) Chebyshev tau-QZ algorithm methods for calculating spectra of hydrodynamic stability problems. App. Num. Math. 22: 399-434

5. Drazin PG, Reid WH (1981) Hydrodynamic stability. Cambridge University Press, London New York Sydney

6. Gardner DR, Trogdon SA, Douglas RW (1989) A modified tau spectral method that eliminates spurious eigenvalues. J. Comput. Phys. 80: 137-167

7. McFadden GB, Murray BT, Boisvert RF (1990) Elimination of spurious eigenvalues in the Chebyshev tau spectral method. J. Comput. Phys. 91: 228-239

8. Moler CB, Stewart GW (1973) An algorithm for generalized matrix eigenproblems. SIAM J. Num. Anal. 10: 241-256

9. Orszag SA (1971) Accurate solution of the Orr-Sommerfeld stability equation. J. Fluid Mech. 50: 689-703 10. Payne LE, Straughan B (2000) A naturally efficient numerical technique for

porous convection stability with non-trivial boundary conditions. Int. J. Numer. Anal. Meth. Geomech. 24: 815-836

11. Reddy SC, Henningson DS (1993) Energy growth in viscous channel flows. J. Fluid Mech. 252: 209-238

12. Straughan B (1992) The energy method, stability, and nonlinear convection. Springer-Verlag, Berlin Heidelberg New York

13. Straughan B (2001) Porous convection, the Chebyshev tau method, and spurious eigenvalues. In: Straughan B, Greve R, Ehrentraut H, Wang Y (eds) Continuum Mechanics and Applications in Geophysics and the Environment. Springer, Berlin Heidelberg New York, pp 140-152