The Origin and Nature of Spurious Eigenvalues in the Spectral tau Method

Paul T. Dawkins Department of Mathematics and Statistics University of Nebraska-Lincoln Lincoln, Nebraska 68588-0323

             pdawkins@math.unl.edu

Steven R. Dunbar Department of Mathematics and Statistics University of Nebraska-Lincoln Lincoln, Nebraska 68588-0323

              sdunbar@math.unl.edu

Rod W. Douglass Computational Mechanics Idaho National Engineering Laboratory Idaho Falls, Idaho 83415-3808

                rxd@inel.gov ^{∗}^{†}^{‡}^{§}

June 24, 1997

Abstract

The Chebyshev-tau spectral method for approximating eigenvalues of boundary value problems may produce spurious eigenvalues with large positive real parts even when all true eigenvalues of the problem are known to have negative real part. We explain the origin and nature of the “spurious eigenvalues” in an example problem. The explanation will demonstrate that the large positive eigenvalues are in fact not spurious but instead are an approximation of infinite eigenvalues in a nearby generalized eigenvalue problem. The explanation also demonstrates that a range of related spectral methods generate large magnitude eigenvalues, both positive and negative, which are approximations of the infinite eigenvalue. We also discuss various ad hoc methods which other authors have proposed for eliminating the spurious eigenvalues from the persepctive of generalized eigenvalue problems.

∗This work was performed, in part, under the auspices of the U. S. Department of Energy, contract DE-AC07-94ID13223, through the support of the INEL Long Term Research Initiative in Computational Mechanics

†Keywords spectral methods, Chebyshev-tau method, spurious eigenvalues, generalized eigenvalue problems

‡PACS classification numbers: 02.10.Sp, 02.30.Mv, 02.60.Dc, 02.60.Lj, 02.70.Hm §1991 Mathematics Subject Classification. 65N35, 65F15, 41A10, 41A50, 42C99, 15A22, 47A75, 15A18

1 Introduction

Spectral methods of numerical solution for some differential eigenvalue problems including hydrodynamic stability problems may produce a set of “spurious eigenvalues” along with approximations to the true eigenvalues. An extensive literature that includes the early work of [11] and later in [3, 5, 9, 16, 21, 24, 25] documents this phenomenon. The intent of the present paper is to uncover the nature and origin of these “spurious eigenvalues” in the context of a typical model problem, proving that they exist at all orders of truncation in the Chebyshev-tau method, proving that “spurious eigenvalues” exist at all orders of truncation in a range of associated spectral methods, establishing how they behave under increases in truncation order for a range of spectral approximation methods, and explaining why they arise. Several authors have asserted that the spurious eigenvalues are a consequence of the “discretization” of the problem, while other authors have stated that the reason for the occurrence of the spurious eigenvalues is unknown. Our explanation will show that the spurious eigenvalues are an approximation of infinite eigenvalues in a generalized eigenvalue problem from a “nearby” Legendre spectral approximation of the boundary value problem. The infinite eigenvalues occur because the boundary conditions fit together with the Legendre spectral approximation, allowing an exceptional solution corresponding to infinite eigenvalues. Our explanation of the origin of the spurious eigenvalues gives a reason ad hoc methods for eliminating the spurious eigenvalues work. Thus, this paper should give users of the Chebyshev-tau method, one of several spectral methods available, guidance in the method’s use for approximating eigenvalues.

Previous authors usually define “spurious eigenvalues” as the positive or unstable eigenvalues when it is obvious that positive eigenvalues are not eigenvalues of the model problem under investigation. Here we show that some negative eigenvalues arising from the ad hoc methods for eliminating spurious eigenvalues share perturbation and growth characteristics with the spurious eigenvalues from the direct Chebyshev-tau method. This provides a clue about the nature and origin of the “spurious eigenvalues” as approximations of infinite eigenvalues in a generalized eigenvalue problem and permits a generalized definition of spurious eigenvalues.

The Chebyshev-tau method and other spectral methods applied to the model problem lead to a generalized eigenvalue problem of the form Ax = λBx. A more symmetric form of the generalized eigenvalue problem is βAx = αBx, where pairs of complex numbers α, β with ratio λ = α/β become a generalized eigenvalue. Stewart and Sun [20] have a complete theory for generalized

Figure 1: Schematic diagram of perturbations of generalized eigenvalues

eigenvalues that we use extensively later in the paper. The important point to notice now is that if B has a non-trivial null-space while A is non-singular, then α = 1, β = 0 is a generalized eigenvalue. This pair (1, 0) is an “infinite eigenvalue”. More precisely, since pairs of complex numbers α, β with common ratio λ represent the same eigenvalue, lines through the origin in the complex plane ℂ^{2} represent generalized eigenvalues. More simply, lines through the origin correspond to points on the unit circle. For the model problem we consider it is simpler yet, since we can take α, β to be real numbers, and so we can picture the generalized eigenvalues as points on the unit circle in ℝ^{2}. The point (1, 0) now acts as a “point at infinity”. Perturbations can then move the eigenvalue into either the upper half-plane, or the lower half-plane, resulting in respectively a large-magnitude positive or negative ratio, see the schematic digrams in Figure 1. This picture describes our generalized definition of spurious eigenvalues as perturbations of infinite eigenvalues, including the large magnitude positive eigenvalues previously observed.

Aside from the spurious eigenvalues, the Chebyshev-tau method computes the remaining eigenvalues accurately and efficiently, accounting for the method’s popularity and utility. The practical effect of the existence of “spurious eigenvalues” is to cast doubt on results obtained from the application of the Chebyshevtau method to a stability analysis of a flow or structure. If positive eigenvalues, indicating exponential growth of a disturbance, which are really spurious result then the analyst may conclude incorrectly that the system is unstable when it is not. It is important to either alter the method to not compute spurious modes or to identify spurious modes so that they can be ignored in the analysis. Works addressing elimination of spurious modes include [9, 13, 16, 25] and those that simply state that spurious modes were computed and then ignored include [3, 22, 24]. The perturbation results in this paper give a more effective means of identifying spurious eigenvalues through their growth in magnitude as a function of truncation order.

We begin in Section 2 by defining the model problem that guides the explanation of the spurious eigenvalues in the spectral-tau method. In Section 3 we then discuss the Chebyshev-tau method applied to the model problem including the form of the resulting linear algebra problem. The next Section 4 describes a parity reduction of the problem that simplifies later work. In Section 5 we then derive the characteristic polynomial, first in the Chebyshev-tau case, then in the more general Gegenbauer-tau case. This allows us to determine the sign and growth of the largest magnitude eigenvalues for the entire range of tau methods. The next Section 6 starts with a reduction to an equivalent basis of monomials instead of Gegenbauer polynomials. The monomial basis is convenient for theory because it reduces the linear algebra to a standard triangular-Hessenberg basis. We can then apply the theory for generalized eigenvalues to show the model problem has an infinite generalized eigenvalue. For the related Gegenbauer-tau spectral methods, we show there is a generalized eigenvalue that is a perturbation of the infinite eigenvalue. We determine the size of the perturbation and find that it matches the growth of the magnitude of the eigenvalues from the characteristic polynomial method derived previously. Section 7 gives the explanation of the origin of infinite generalized eigenvalue in the Legendre-tau method. The boundary conditions and the form of the Legendre polynomials coincide permitting the existence of a special solution corresponding to the infinite eigenvalue. We then briefly discuss the ad hoc methods for the elimination of positive eigenvalues in Section 8. Finally, we present some conclusions based on this paper’s results. Some necessary background facts and mathematical proofs are in appendices.

2 The Model Problem

We will consider the model problem

\begin{aligned}u^{(4)} & = su′′,  −1 < x < 1, \\ u(−1) & = u(1) = u′(−1) = u′(1) = 0 \end{aligned}  (1)

throughout the paper as an example of the occurrence of spurious eigenvalues when applying the Chebyshev-tau method to boundary value problems. The model problem occurs in [11, pages 143-144] from a separation of variables applied to a one-dimensional model of the vorticity-streamfunction equations for low Reynolds number incompressible flow. This model problem is also considered in [9, 13, 16].

A similar boundary value problem occurs in hydrodynamic stability analysis in spherical coordinates:

D^{2}u = sDu,  r_{1} < r < r_{2},
u(r_{1}) = u(r_{2}) = u′(r_{1}) = u′(r_{2}) = 0

where the second-order differential operator

Du = u′′(r) + \frac{2}{r}u′ − \frac{l(l + 1)u}{r^{2}}

occurs instead of the second-derivative operator. The operator D arises from the spectral decomposition of the Laplacian in spherical coordinates. This problem also has spurious eigenvalues when applying the Chebyshev-tau method. See [9, 10, 23] for details. The structural similarity to the model problem (1) is obvious.

Another example is the Orr-Sommerfeld stability equation for plane Poiseuille flow

[u^{(4)} − 2α^{2}u′′ + α^{4}u]/(−iαR) + [(U − s)(u′′ − α^{2}u) − U′′u] = 0, −1 < x < 1,

with boundary conditions

u(−1) = u(1) = u′(−1) = u′(1) = 0

where u is the amplitude of the velocity disturbance; α is the wavenumber; R is the Reynolds number, the stability parameter for this problem; and U(x) = 1 − x^{2} is the known steady base flow whose stability is being examined. It is important to know the value of R where at least one eigenvalue first has a positive real part as R is increased from zero, because for this eigenvalue the disturbance will grow exponentially in time instead of being damped out. For more details see [9, 17, 24]. In this problem, [u^{(4)} − 2α^{2}u′′ + α^{4}u]/(−iαR) is a fourth-order differential operator and [(U − s)(u′′ − α^{2}u) − U′′u] is a secondorder differential operator. The structural similarity to the model problem (1) is also obvious. Using the Chebyshev-tau spectral method produces two spurious eigenvalues with large imaginary parts, these spurious eigenvalues are clearly recognized by the large magnitudes of the values of at least one of the tau coefficients, see [9] for details.

The eigenvalues of the model problem (1) are all negative and satisfy either s = −n^{2}π^{2} or tan(\sqrt{−s}) = \sqrt{−s}. The first five eigenvalues are (numerically)

\begin{aligned}s_{1} & = −π^{2} ≈ −9.869604401 \\ s_{2} & ≈ −20.19072856 \\ s_{3} & = −4π^{2} ≈ −39.47841760 \\ s_{4} & ≈ −59.67951594 \\ s_{5} & = −9π^{2} ≈ −88.826439612. \end{aligned}  (2)

The corresponding eigenfunctions are respectively

1 + (−1)^{n+1} cos(nπx)

and

−\sqrt{−s} cos(\sqrt{−s}x) + sin(\sqrt{−s}x).

3 The tau method and the model problem

Lanczos [14, 15] first proposed the tau method as a means of solving boundary value problems without requiring basis functions to satisfy the boundary conditions. Fox [7, 8] later extensively developed its use with Chebyshev polynomials. Ortiz [18] extended the method with the use of canonical polynomials. Orszag [17] applied and advocated the Chebyshev-tau method for a wide variety of problems. The tau method uses a truncated series expansion in a complete set of orthogonal functions as an approximation for the solution of an ordinary differential equation. Chebyshev polynomials work well for this approximation technique because of their nearly optimal uniform approximation of continuous functions, their orthogonality and completeness, and other extremal properties, [19]. Boyd [4, Chapter 18] has a useful and interesting overview of the history and philosophy of the tau method.

3.1 The tau method

To illustrate the use of the tau method in its usual form, consider an approximate truncated series solution to the boundary value problem

L{u} = 0,  −1 < x < 1  (3)
B_{i}{u(−1)} = B_{j}{u(1)} = 0,  i + j = N_{b}.  (4)

Here L is a linear ordinary differential operator of order N_{b}, u is an unknown function in the variable x, and B_{i} and B_{j} are linear operators such that B_{i}{u(−1)} = 0 and B_{j}{u(1)} = 0 represent the boundary conditions. The subscripts i and j denote the number of boundary conditions applied at x = −1 and x = 1 respectively. Let u be approximated by u(x, N),

u(x, N) = ∑_{k=0}^{N+N_{b}} a_{k}f_{k}(x)  (5)

where the a_{k}’s are unknown coefficients and the f_{k}’s are functions from the complete orthogonal set corresponding to the inner product symbolized by ⟨·, ·⟩. We will use the family of Gegenbauer (or ultraspherical) polynomials as the complete orthonormal set of functions. The family of Gegenbauer polynomials include the Legendre polynomials and Chebyshev polynomials of the first and second kind as special cases, see Appendix B for more details. Substitute (5) into (3), and expand the derivatives in terms of the orthogonal functions. Using the inner product with the functions f_{l}, l = 0, 1, 2, …, N yields a system of N + 1 algebraic equations for the unknown coefficients a_{k}, k = 0, 1, 2, …, N + N_{b}. The remaining N_{b} equations result from substituting (5) into (4) so that the following system of N + N_{b} + 1 equations, linear in the a_{k}’s, is produced:

∑_{k=0}^{N+N_{b}} a_{k}⟨L{f_{k}(x)}, f_{l}⟩ = 0,  l = 0, 1, 2, …, N  (6)

∑_{k=0}^{N+N_{b}} B_{i}{f_{k}(−1)}a_{k} = ∑_{k=0}^{N+N_{b}} B_{j}{f_{k}(1)}a_{k} = 0,  i + j = N_{b}. (7)

If the problem is an eigenvalue problem, then the unknown eigenvalue appears in the corresponding matrix equation. Standard matrix and eigenvalue solver packages will then solve equations (6), (7).

The philosophy of the tau-method is to solve a nearby or approximate differential equation exactly in a polynomial. The method can also be described as “Galerkin’s method with boundary bordering”, [4, Chapter 18]. The tau method applied to a constant coefficient linear differential equation solves the following residual problem exactly:

L{u} = ∑_{k=1}^{N_{b}} τ_{k}f_{N+k}(x)  (8)

∑_{k=0}^{N+N_{b}} B_{i}{f_{k}(−1)}a_{k} = ∑_{k=0}^{N+N_{b}} B_{j}{f_{k}(1)}a_{k} = 0,  i + j = N_{b}.  (9)

where the coefficients τ_{k} for k = 1, …, N_{b} are unknowns, called the tau coefficients. The residual is a measure of ability of u(x, N) to satisfy the original differential equation. The coefficients of the residual, traditionally using the notation τ_{k}, account for the name of the method.

Using the inner product for the orthogonal functions f_{k}, equation (6) reduces to the system of equations

⟨L{u(x, N)}, f_{k}⟩ = 0,  k = 0, 1, … N  (10)
⟨L{u(x, N)}, f_{N+k}⟩ = τ_{k}⟨f_{N+k}, f_{N+k}⟩,  k = 1, … N_{b}.  (11)

After solving (6) and (7), the tau coefficients can be computed from (11). Fox [7] has shown that the tau coefficients provide error bounds on the approximate solution u(x, N). Applying the algorithm by taking inner products with the orthogonal polynomials does not explicitly involve the tau coefficients, so often the linear system is solved directly without reference to the perturbation. For the Chebyshev-tau method the error introduced by the tau method will be small if the magnitudes of the tau coefficients are small, since the Chebyshev polynomials are bounded in absolute value by one on the interval [−1, 1]. In common practice the tau coefficients are not computed, although Gardner et al.

[9] used the tau coefficients to identify spurious eigenvalues.

3.2 Application of the Chebyshev-tau method to the model problem

To further illustrate the application of the usual tau method using Chebyshev polynomials, consider the model eigenvalue problem (1). We will follow the notation of [19] for the Chebyshev polynomials. Let T_{k}(x) = cos(k arccos x) for −1 ≤ x ≤ 1 be the Chebyshev polynomial of degree k. The Chebyshev polynomials form a complete orthogonal set with respect to the inner product

⟨f(x), g(x)⟩ = \frac{2}{π} ∫_{−1}^{1} f(x)g(x) \frac{dx}{\sqrt{1 − x^{2}}} .

The Chebyshev polynomials are a special case of the Gegenbauer polynomials G_{n}^{ν}(x) with ν = 0, so this example illustrates a case of the Gegenbauer-tau method. To each function f(x) integrable on [−1, 1], there is an associated Chebyshev expansion

f(x) ∼ \begin{matrix}∞ \\ ∑ & ′ \\ k=0 \end{matrix} A_{k}T_{k}(x) =^{def} A_{0}/2 + A_{1}T_{1}(x) + A_{2}T_{2}(x) + …

where

A_{k} = ⟨f(x), T_{k}(x)⟩ = \frac{2}{π} ∫_{−1}^{1} f(x)T_{k}(x) \frac{dx}{\sqrt{1 − x^{2}}} .

In certain cases, the Chebyshev expansion is optimal among all expansions in terms of Gegenbauer polynomials, see Appendix B for a precise statement of the details. This optimality accounts for the common use of Chebyshev polynomials for numerical approximations of all kinds.

The differential equation (1) becomes

a_{k}^{(4)} − sa_{k}^{(2)} = 0  k = 0, 1, 2, …, N  (12)

where a_{k}^{(2)} and a_{k}^{(4)} are linear combinations of the unknown coefficients. These linear combinations come from well-known formulas used to expand derivatives of Chebyshev polynomials as expansions in Chebyshev polynomials, see [9, 11, 19]. The formulas are given and applied to the model problem in Appendix A. The boundary conditions are also linear combinations of the unknown coefficients, using the values of the Chebyshev polynomials and their derivatives at the endpoints −1 and +1.

The system of equations in matrix-vector form is

Aa = sBa.  (13)

For example, if N = 1, then A and B are

A = (\begin{matrix}0 & 0 & 0 & 0 & 192 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1920 \\ 1 & −1 & 1 & −1 & 1 & −1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & −4 & 9 & −16 & 25 \\ 0 & 1 & 4 & 9 & 16 & 25 \end{matrix}) , B = (\begin{matrix}0 & 0 & 4 & 0 & 32 & 0 \\ 0 & 0 & 0 & 24 & 0 & 120 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}) . (14)

N  1  2  3  4  5  6
eigenvalues 40 251.453 529.388 1331.132 2284.474 4272.17
12 40.  251.453 529.388 1331.132 2284.47
−11.453 −11.453 −9.892 −9.892 −9.870
−25.388 −25.388 −20.338 −20.338
−61.239 −61.239 −40.663
−104.136 −104.136
−189.638

Table 1: Eigenvalues, including spurious eigenvalues, for the model problem (1) computed with the Chebyshev-tau method

Table 1 gives the eigenvalues from the various orders of Chebyshev-tau approximation of the boundary value problem (1). Note that 2 large positive eigenvalues result from the calculations, and that they increase as the order of approximation increases. The negative eigenvalues appear to be converging to the values (2) of the boundary value problem (1). Note also that the eigenvalues “leap-frog”, that is, the largest eigenvalue in column N becomes the second largest positive eigenvalue in column N + 1. Likewise the least eigenvalue in column N becomes the next-to-least eigenvalue in column N + 1. An abridged table of eigenvalues for higher truncation orders is given in [11, 16]. The entry for N = 6 common to the table above and to [11, 16] agree. However, the interest here is not in extending the table to high orders of approximation, but rather to explain the origin and nature of the two large positive eigenvalues that are spurious for the boundary value problem (1).

4 Parity reduction of the problem

The parity of the model problem (1) allows a reduction of the problem for theoretical purposes that exposes the more important and interesting phenomenon of the spurious eigenvalues in a clearer fashion. Recall that odd index Gegenbauer polynomials are odd polynomials and even index Gegenbauer polynomials are even, see Appendix B, for all values of the parameter ν. Since the derivatives in (1) are of even order, any odd index Gegenbauer polynomial inserted into the differential equation will remain odd after differentiations. The re-expansion of the derivatives in Gegenbauer polynomials will be in terms of odd indices of lower order. This has the effect of further partitioning the problem (6), (7) or equivalently (13) into odd and even portions. Such a parity reduction is common in spectral methods, see [4, Chapter 7]. It will suffice to examine only one subportion, say the even portion of the problem.

A specific application of this reduction to the Chebyshev-tau method explains the “leap-frogging” of eigenvalues observed in Section 3.2. For say M odd, M = 2N + 1, the problem (6), (7) factors into an “odd” equation of size (N +1)×(N +1) and an “even” equation of size (N +1)×(N +1). The odd portion comes from the odd Chebyshev polynomials T_{1}(x), T_{3}(x), …, T_{2N+1}(x). Increasing the size of the approximation (5) to M +1, we would partition the problem (6), (7) or equivalently (13) into an “odd” portion of size (N + 1) × (N + 1) and an “even” portion of size (N +2)×(N +2). The odd portion still comes from the odd Chebyshev polynomials T_{1}(x), T_{3}(x), …, T_{2N+1}(x) as in the previous problem. This accounts for the leap-frogging of the eigenvalues in the previous table. Those eigenvalues occurring in problem M + 1 which also previously occurred in problem M are the recalculation of the eigenvalues from the subportion of the problem that was unchanged in passing from M to M + 1. The leapfrogging occurs only for the model problem (1) and is not true for more general eigenvalue problems containing mixed even and odd orders of differentiation, see [9].

5 The characteristic polynomial and the lower bound for the positive eigenvalue

We derive a general expression for the characteristic polynomial of the eigenvalue problem arising from a range of spectral-tau methods. In particular, we obtain the characteristic polynomials from the Chebyshev-tau method and the Legendre-tau method. The results will show that a range of spectral methods generate positive (or spurious) eigenvalues, increasing to an infinite eigenvalue from the Legendre-tau method. Prof. Hans Weinberger suggested this derivation of the characteristic polynomial to us and we thank him for his assistance. We can also derive the same characteristic polynomial directly from determinants of the triangular-Hessenberg form of the matrix generalized eigenvalue problem through an intricate sequence of recursion arguments and reductions.

As an example, we first derive the characteristic polynomial for the Chebyshevtau method. This has the advantage of being simple and explicit since the properties of the Chebyshev polynomials that we use are well known. Then for full generality we derive the characteristic polynomial for the Gegenbauer-tau method.

5.1 The characteristic polynomial for the Chebyshev-tau method

Consider the residual problem (8) of the spectral-tau method in Section 3.1 reduced by parity as in Section 4 applied to the model problem (1).

L{u} = u^{(4)} − su′ = τ_{1}f_{2N+2}(x) + τ_{2}f_{2N+4}(x).

Use the Chebyshev polynomials as the set of orthogonal polynomials. Specifically, we want to choose the coefficients of the even polynomial u of degree 2N + 4 in such a way that the residual of L{u} = u^{(4)} − su′′ (which is an even polynomial) is orthogonal to T_{0}, T_{2}, T_{4}, … T_{2N}. Then the residual will only be a multiple of T_{2N+2}.

We write

L{u} = (D^{2} − s)D^{2}{u} = (−D^{2} + s){−u′′}  (15)

where D is the derivative operator Because −u′′ is a polynomial, we use the Neumann series for the inverse of the operator −D^{2} + s to find that

−u′′ = τ_{2N+2} ∑_{k=0}^{N+1} s^{−k−1}D^{2k}T_{2(N+1)}.  (16)

This equation has a polynomial solution that satisfies the boundary conditions if and only if the integrals of 1 and of x times the right-hand side are zero. The integral of x times the right-hand side is zero by symmetry. Multiplying by s^{N+2} to create a polynomial, the condition for s to be an eigenvalue is that

s^{N+1} ∫_{−1}^{1} T_{2(N+1)}(x) dx + 2 ∑_{k=1}^{N+1} s^{N+1−k}T_{2(N+1)}^{(2k−1)}(1) = 0.  (17)

This is the characteristic equation. Using [19, page 5, Exercise 1.1.4]

∫_{−1}^{+1} T_{2(N+1)}(x) dx = −2/((2(N + 1))^{2} − 1).

Thus the leading coefficient is always negative. From [19, page 38, Exercise 1.5.6]

T_{n}^{(j)}(1) = (∏_{i=0}^{j−1}(n^{2} − i^{2})) / (∏_{i=0}^{j−1}(2i + 1)) .

The characteristic polynomial is then

p(s) =^{def} \frac{−2}{(2(N + 1))^{2} − 1}s^{N+1} + 2 ∑_{k=1}^{N+1} \frac{∏_{i=0}^{2k−2}((2(N + 1))^{2} − i^{2})}{∏_{i=0}^{2k−2}(2i + 1)}s^{N+1−k}.

Notice that all coefficients except the leading coefficient are positive. We will use this information to find a lower bound on the unique positive root of the characteristic polynomial.

Since all coefficients except the first are positive, we can then say for s > 0,

p(s) > −\frac{2}{(2(N + 1))^{2} − 1}s^{N+1}+2(2N+2)^{2}s^{N} +2\frac{∏_{i=0}^{2N}((2(N + 1))^{2} − i^{2})}{∏_{i=0}^{2N}(2i + 1)} ≡ q(s)

Taking the derivative of q(s) with respect to s and solving for s gives a single positive critical point

s_{cr} = 2N(2N + 1)(2N + 2)(2N + 3) > (2N + 1)^{4}.

The characteristic polynomial p(s) is always greater than q(s) for all s > 0 since all coefficients except the leading coefficient are positive. In particular, p(s) ≥ q(s) > 0 for 0 ≤ s ≤ s_{cr}. We also know that due to the negative leading coefficient, the characteristic polynomial p(s) will eventually be negative and so by the intermediate value theorem the characteristic polynomial must have a root greater than (2N+1)^{4}. So we have a lower bound on the positive eigenvalue, i.e. s > (2N + 1)^{4}.

5.2 The characteristic polynomial for the general Gegenbauertau method

Again consider the general residual problem (8) in Section 3.1 reduced by parity as in Section 4 applied to the model problem (1). Use the Gegenbauer polynomials as the set of orthogonal polynomials for the tau method.

Proceed just as in equations (15), (16), (17) to obtain the characteristic polynomial. The resulting characteristic polynomial is

s^{N+1} ∫_{−1}^{1} G_{2(N+1)}^{ν}(x) dx + 2 ∑_{k=1}^{N+1} s^{N+1−k} (G_{2(N+1)}^{ν})^{(2k−1)} (1) = 0.

We will deduce the special limiting case of ν = 0 reducing to the previously derived Chebyshev-tau characteristic polynomial later.

The following lemmas derive the coefficients of the characteristic polynomial explicitly.

Lemma 1 Let n ≥ 0 be an integer and ν ∈ (−\frac{1}{2}, ∞), ν ≠ 0. Then

∫_{−1}^{1} G_{n}^{ν}(x) dx = [(−1)^{n} + 1] \frac{(2ν − 1) Γ(n + 2ν − 1)}{(n + 1)! Γ(2ν)}

A complete proof is in [6], or the reader can check the lemma by examining special cases. The proof proceeds by writing the integral in terms of hypergeometric functions and using hypergeometric function identities.

Lemma 2 Let n, k ≥ 0 be integers and ν ∈ (−\frac{1}{2}, ∞), ν ≠ 0, then,

\frac{d^{k}}{dx^{k}} G_{n}^{ν}(x)|_{x=1} = 2^{k}\binom{n + k + 2ν − 1}{n − k} ∏_{j=0}^{k−1}(ν + j)

A complete proof is in [6].

Notice that if the order of the differentiation, k, is greater than the order of the Gegenbauer polynomial, n, Lemma 2 will give the correct answer. In this case n − k < 0, and since this term is in the bottom of the binomial coefficient, the binomial coefficient will be zero.

Lemmas 1 and 2 give the characteristic polynomial:

p(s) = \frac{2}{Γ(2ν)} (\frac{Γ(2N + 2ν + 1)(2ν − 1)}{(2N + 3)!}) s^{N+1}

+ 2 ∑_{k=1}^{N+1} s^{N+1−k}2^{2k−1}\binom{2N + 2k + 2ν}{2N − 2k + 3} ∏_{j=0}^{2k−2}(ν + j). (18)

Note that these lemmas do not cover the ν = 0 case, when the Gegenbauer polynomials are the Chebyshev polynomials. In Lemma 1 there is a Γ(2ν) in the denominator of the formula. For ν = 0, Γ(2ν) has an pole singularity. An intuitive interpretation would then make the leading coefficient of the characteristic polynomial zero for all n. Similarly, for ν = 0 the product in the formula from Lemma 2 yields 0 because of the j = 0 factor. All coefficients of the characteristic polynomial would then be 0. To recover the Chebyshev-tau characteristic polynomial at ν = 0 from the general Gegenbauer-tau characteristic polynomial we must use a limiting argument.

5.3 Spurious eigenvalues from the Gegenbauer-tau method

We can now use the characteristic polynomial (18) to prove the existence of a spurious eigenvalue for a range of ν’s. We will examine the ν = 0 case corresponding to the Chebyshev polynomials of the first kind later (see Section 5.4).

First examine the leading coefficient. The limits on N and ν, namely, N ≥ 0 and ν > −\frac{1}{2}, imply Γ(2N + 2ν + 1) > 0. Table 2 completes the sign analysis for the leading coefficient.

−\frac{1}{2} < ν < 0 0 < ν < \frac{1}{2} ν = \frac{1}{2}  ν > \frac{1}{2}
Γ(2ν)  negative  positive positive positive
2ν − 1  negative  negative  zero  positive
leading coefficient  positive  negative  zero  positive

Table 2: Sign Analysis on the Leading Coefficient of the Characteristic Polynomial

The sign analysis on the non-leading coefficients is much simpler. Since N ≥ 0, k ≥ 1, and ν > −\frac{1}{2} we have 2N + 2k + 2ν > 0. Therefore, the binomial coefficient will always be positive. The sign of the non-leading coefficient then will be determined by the sign of the product, which is in turn determined by the j = 0 term. The non-leading coefficients will be negative for −\frac{1}{2} < ν < 0 and positive for ν > 0.

With these facts we have the following theorem:

Theorem 1 For −\frac{1}{2} < ν < \frac{1}{2}, ν ≠ 0, the model problem (1) solved using the Gegenbauer-tau method to (8), (9) and reduced by parity to the even portion will have a single positive eigenvalue, λ_{+} where

λ_{+} > \frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}.

Proof: First look at the case of 0 < ν < \frac{1}{2}. For this case the leading coefficient is negative and the non-leading coefficients are all positive.

Dropping all of the terms in the characteristic polynomial except for the k = 1 and k = N + 1 terms will give a polynomial

q(s) = \frac{2}{Γ(2ν)} (\frac{Γ(2N + 2ν + 1) (2ν − 1)}{(2N + 3)!}) s^{N+1}

+ 4ν\binom{2N + 2ν + 2}{2N + 1}s^{N} + 2^{2N+2} (4N + 2ν + 2) ∏_{j=0}^{2N}(ν + j).

Note that q(0) = p(0) and that q(s) < p(s) for s > 0.

The polynomial q(s) has only one critical point, s_{cr}.

s_{cr} = −4νN\binom{2N + 2ν + 2}{2N + 1}\frac{Γ(2ν) (2N + 3)!}{(2N + 2) Γ(2N + 2ν + 1) (2ν − 1)}.

Since (2ν−1) < 0, s_{cr} > 0. Now, eliminating the gamma functions and factoring gives

s_{cr} = \frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}.

The polynomial q(s) is positive at s = 0 and the leading coefficient is negative and so the q(s) will eventually be negative. Therefore, since s_{cr} is the only positive critical point it must be a maximum. This along with the fact that q(s) will eventually be negative implies that q(s) must have a single root in the right half line which is greater than s_{cr}.

Like q(s) the characteristic polynomial p(s) is positive at s = 0 and for sufficiently large s will be negative. Also, by definition p(s) > q(s). Therefore, the characteristic polynomial will have a root greater than the root of q(s), which is in turn greater than s_{cr}. Hence, the Gegenbauer tau approximation to the solution of the model problem will have an eigenvalue λ_{+} such that

λ_{+} > s_{cr} = \frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}  (19)

Also, once the characteristic polynomial becomes negative it will remain negative and so there will be no other roots of the characteristic polynomial in the right half line.

For the case −\frac{1}{2} < ν < 0 the leading coefficient is positive and all other coefficients are negative. In this case the characteristic polynomial is simply the negative of the characteristic polynomial from the first case. Therefore, when −\frac{1}{2} < ν < 0, the Gegenbauer-tau approximation to the solution of the model problem will also have an eigenvalue λ_{+} such that

λ_{+} > s_{cr} = \frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}  (20)

□

Note that as ν → ±\frac{1}{2} we have s_{cr} → ∞ and so λ_{+} → ∞.

5.4 Reduction of the Gegenbauer-tau characteristic polynomial to the Chebyshev-tau characteristic polynomial (ν = 0)

Recall the definition of the Chebyshev polynomials in terms of Gegenbauer polynomials [1, 2]

T_{n}(x) = \frac{n}{2} lim_{ν→0} \frac{G_{n}^{ν}(x)}{ν}

The main difficulty when examining the ν = 0 case is that G_{n}^{0}(x) = 0 for all n and x. The limit definition takes care of this problem, but introduces other problems.

In order to handle this limit we need two facts:

lim_{ν→0} ν Γ(aν) = \frac{1}{a} for a > 0

and for positive integers m

\frac{Γ(ν + m)}{ν Γ(ν)} → Γ(m)

uniformly as ν → 0.

These facts imply that the convergence of the limit in the Gegenbauer polynomial definition of the Chebyshev polynomials is uniform [6],

\frac{n G_{n}^{ν}(x)}{2ν} →^{unif} T_{n}(x) as ν → 0

Uniform convergence will allow this limit to be interchanged with integrals and derivatives, for example

\frac{d}{dx} (\frac{nG_{n}^{ν}(x)}{2ν}) →^{unif} T′_{n}(x) as ν → 0.

First calculate the leading coefficient. The uniformity of the limit will first be used to interchange the limit and integration.

∫_{−1}^{1} T_{2N+2}(x) dx = ∫_{−1}^{1} \frac{2N + 2}{2} lim_{ν→0} \frac{G_{2N+2}^{ν}(x)}{ν} dx

= (N + 1) lim_{ν→0} \frac{1}{ν} ∫_{−1}^{1} G_{2N+2}^{ν} dx

= −\frac{2}{(2N + 3)(2N + 1)}.

Likewise calculate the remaining coefficients.

\frac{d^{2k−1}}{dx^{2k−1}} T_{2N+2}(x)|_{x=1} = \frac{d^{2k−1}}{dx^{2k−1}} (\frac{2N + 2}{2} lim_{ν→0} \frac{G_{2N+2}^{ν}(x)}{ν})|_{x=1}

= \frac{2N + 2}{2} lim_{ν→0} \frac{1}{ν} (\frac{d^{2k−1}}{dx^{2k−1}} G_{2N+2}^{ν}(x))|_{x=1}

= \frac{∏_{j=0}^{2k−2}((2N + 2)^{2} − j^{2})}{∏_{j=0}^{2k−2}(2j + 1)}

Thus, we can recover the Chebyshev-tau characteristic polynomial from the Gegenbauer-tau characteristic polynomial.

As in Section 5 the Chebyshev-tau approximation to the solution of the model problem has a spurious positive eigenvalue, λ_{+} such that

λ_{+} > 2N(2N + 1)(2N + 2)(2N + 3)

Notice that this is precisely the lower bound found in Theorem 1 with ν = 0. We have the following theorem:

Theorem 2 For −\frac{1}{2} < ν < \frac{1}{2} the model problem (1) reduced by parity and solved using the Gegenbauer-tau method will have a single positive eigenvalue, λ_{+} where

λ_{+} > \frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}

6 Mathematical reductions and reformulations

The matrix formulation of the tau method in Section 3.2 is inconvenient for theoretical analysis for two reasons. The matrices A and B are dense, with a partial “checkerboard” structure. Furthermore, because of the “boundary bordering” (see [4]) the matrix B is singular. The singularity of B requires a partitioning process for the determination of the eigenvalues, see [9].

In this section, we describe an alternative but equivalent spectral formulation of the problem (1). The alternative formulation results in a non-singular matrix system with a triangular-Hessenberg structure that is convenient for theoretical analysis. The triangular-Hessenberg form helps to explain the nature and origin of the spurious eigenvalues.

In the following subsections we briefly describe the steps that lead to the triangular-Hessenberg form. We impose additional modifications of the procedure that we outlined in Sections 3.1 and 3.2 in arriving at (13). The modifications do not alter the problem or its eigenvalues, but do take advantage of various properties to eliminate unnecessary or obscuring information from the problem. First, we have already shown in Section 4 how to factor the original procedure into odd and even parts, reducing by half the size of the problem to be solved. Next, we show how to automatically incorporate the boundary conditions into the process to eliminate the singular rows of the matrix eigenvalue problem. Finally, we show that the original problem leads to a special triangular-Hessenberg form by using monomial basis functions instead of the (equivalent) Chebyshev or Gegenbauer polynomials. It is in this final context that the origin of the spurious eigenvalues is understood.

6.1 Galerkin-Gegenbauer bases, monomial bases, Galerkinmonomial bases

Finding a polynomial of degree N + 4 which satisfies the differential equation and boundary conditions in the sense of having residual orthogonal to G_{0}^{ν}(x), …, G_{N}^{ν}(x) approximately solves the model problem (1). Of course, users most commonly apply this approximation with ν = 0, the Chebyshev polynomials, for reasons discussed in Appendix B. However, investigation of the spurious eigenvalues requires increased generality with the full range of Gegenbauer polynomials.

The boundary conditions imply that the approximating polynomial will have a root of multiplicity 2 at x = 1 and also at x = −1. This means that it is possible to factor (1 − x)^{2}(1 + x)^{2} = (1 − x^{2})^{2} from the polynomial solution. Therefore, the final solution will have the form

u(x, N) = (1 − x^{2})^{2}(d_{0}G_{0}^{ν}(x) + ⋯ + d_{N}G_{N}^{ν}(x))
= d_{0}(1 − x^{2})^{2}G_{0}^{ν}(x) + ⋯ + d_{N}(1 − x^{2})^{2}G_{N}^{ν}(x)

for some constants d_{0}, …, d_{N}. Thus we may solve for the approximation in terms of the basis functions

(1 − x^{2})^{2}G_{k}^{ν}(x).  (21)

This is a “basis recombination” as in [4, Section 6.5]. That is, for the model problem (1) we are able to recast the Gegenbauer-tau method (or more commonly, the Chebyshev-tau method) into a Galerkin-Gegenbauer-tau method (see also the discussion of nomenclature in [4, Chapter 18]). For example, in the simplest case N = 1, the matrix-eigenvalue problem derived by using the Galerkin-Chebyshev basis elements (21) is

(\begin{matrix}24 & 0 \\ 0 & 120 \end{matrix}) (\begin{matrix}d_{0} \\ d_{1} \end{matrix}) = s (\begin{matrix}2 & 0 \\ 0 & 3 \end{matrix}) (\begin{matrix}d_{0} \\ d_{1} \end{matrix})

The eigenvalues of the system using the Galerkin-Chebyshev basis functions (21) are the same as the eigenvalues of the system (13) using the Chebyshev polynomials and the subsequent reduction (see [9]) with the boundary conditions.

One final theoretical consideration will also simplify the problem. We will restrict our attention to even basis functions. For the theoretical analysis we will use basis functions of the form (1−x^{2})^{2}x^{2j}. That is, we use an even polynomial incorporating the boundary conditions (1−x^{2})^{2}p_{N}(x) = (1−x^{2})^{2} ∑_{k=0}^{N} a_{k}x^{2k} as the approximation to the solution of the boundary value problem instead of (5) with the Gegenbauer polynomials. Inserting this into the differential equation, then taking the inner product with each of the even Gegenbauer polynomials G_{0}^{ν}(x), …, G_{2N}^{ν}(x) yields N + 1 equations in the N + 1 unknowns a_{k}. This is a “Petrov-Galerkin method”, see [4, Chapter 18, page 598]. This amounts to a change of variables in the formulation of the problem. The change of basis can be expressed by using the change of basis from the Gegenbauer polynomials to the monomials. We will denote the resulting matrix equation as

L_{N+1}^{ν}a = sR_{N+1}^{ν}a.  (22)

The superscript indicates the index of the Gegenbauer polynomials used, and the subscript indicates the order of approximation. When it is clear from context, or unnecessary, we will omit the indices. Formulas for the entries of L and R are in Appendix C for the Chebyshev case, and Appendix D for the general Gegenbauer case. Note that L is an upper triangular matrix, and R is an upper Hessenberg matrix. This is another simplification in the problem, since the standard QZ algorithm for the solution of generalized eigenvalue problems first reduces the problem to the triangular-Hessenberg form (which is always possible) and then proceeds to solve the resulting generalized eigenvalue problem, [20]. Note that this reduction is specific to the special form of the model problem and its boundary conditions. Although the monomial basis is ill-conditioned for numerical computation, using the monomial basis for reduction to triangular-Hessenberg form does help expose the central point of this theoretical investigation which is the nature and origin of the spurious eigenvalues.

Increasing the order of approximation from 2N to 2(N + 1) by using the polynomial (1−x^{2})^{2} ∑_{k=0}^{N+1} a_{k}x^{2k} adds a column to both L and R corresponding to the inner products of the derivatives of (1 − x^{2})^{2}x^{2(N+1)} with G_{2j}^{ν}(x) and a row corresponding to the inner products of derivatives of (1 − x^{2})^{2}x^{2i} with G_{2(N+1)}^{ν}(x). This means that as the size of the approximation increases, the matrices L and R need not be completely recalculated. It also means that for a given degree 2N of approximation, the principal submatrices of L and R contain the lesser degrees of approximation. Both properties are useful for the theoretical analysis in this paper.

6.2 The Gegenbauer-tau method yields a regular generalized eigenvalue problem

In discussing generalized eigenvalue problems, we will generally follow the discussion and notation of [20, Chapter VI]. For α, β ∈ ℂ let (α, β) ≠ (0, 0). Then for any complex scalar γ, ⟨α, β⟩ =^{def} {γ(α, β)^{T} : γ ∈ ℂ}. Write ⟨λ⟩ = ⟨λ, 1⟩ and define ⟨∞⟩ = ⟨1, 0⟩.

For vector spaces V_{1} and V_{2}, x ∈ V_{1}, and operators A, B with A, B : V_{1} → V_{2}, if

βAx = αBx

for (α, β) ≠ (0, 0) and x ≠ 0, then ⟨α, β⟩ is a (generalized) eigenvalue of the pair (A, B) with (right) (generalized) eigenvector x.

The matrix pair (A, B), with both matrices square of the same dimensions, is singular if for all (α, β), det(βA − αB) = 0. Otherwise the pair is regular. If either det(A) ≠ 0 or det(B) ≠ 0, then the pair (A, B) is regular. If the pair of N × N matrices (A, B) is regular, then there are precisely N generalized eigenvalues.

In triangular-Hessenberg form it is simple to show that the matrix generalized eigenvalue problem (22) resulting from the Gegenbauer-tau method is regular. Since L^{(ν)} is upper triangular, it suffices to show that the diagonal entries are non-zero. From Appendix D the diagonal entries of L^{(ν)} are

(2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{−2j} \frac{(2j + ν)Γ(ν)Γ(2j + 1)}{Γ(2j + ν + 1)}.  (23)

Thus, there are N + 1 (generalized) eigenvalues for the matrix pair. We may find the eigenvalues as the roots of the characteristic polynomial

det(L^{(ν)} − λR^{(ν)}).

From Appendix D we can show that for Legendre polynomials, when ν = 1/2, the first row of the matrix R_{N+1}^{(1/2)} is zero, see also Section 6.4. Hence, R_{N+1}^{(1/2)} is singular, with a rank of N. Therefore there is only one linearly independent right null vector for R_{N+1}^{(1/2)}. Then we have

Theorem 3 The matrix pair (L_{N+1}^{(1/2)}, R_{N+1}^{(1/2)}) has one infinite generalized eigenvalue ⟨1, 0⟩ = ⟨∞⟩.

6.3 Application of generalized eigenvalue perturbation theory

We use generalized eigenvalue perturbation theory to show that for 0 ≤ ν < 1/2 the system (22) has an eigenvalue that approaches ⟨∞⟩ as the order of approximation increases. We first need more background from [20]. Let (A, B) be a regular matrix pair with simple generalized eigenvalue ⟨α, β⟩ and left and right eigenvectors y and x. Let (Ã, B̃) = (A+E, B+F) be a regular matrix pair and a perturbation of (A, B) with corresponding generalized eigenvalue ⟨α̃, β̃⟩. Let ϵ = \sqrt{∥E∥^{2} + ∥F ∥^{2}} where we will define the matrix norm ∥ · ∥ shortly.

The chordal metric χ(⟨·, ·⟩, ⟨·, ·⟩) for generalized eigenvalues has the form

χ(⟨α, β⟩, ⟨γ, δ⟩) = \frac{|αδ − βγ|}{\sqrt{|α|^{2} + |β|^{2}}\sqrt{|γ|^{2} + |δ|^{2}}} .

Stewart and Sun [20, Theorem 2.1, page 293] show that for small perturbations of a matrix pair, for each eigenvalue ⟨α, β⟩ of the original matrix pair (A, B), there will be a corresponding eigenvalue ⟨α̃, β̃⟩ for the matrix pair (Ã, B̃). More precisely, let (A, B) be a regular pair, and let its eigenvalues

be ⟨λ_{1}⟩, …, ⟨λ_{n}⟩. Then there is an ordering ⟨λ̃_{1}⟩, …, ⟨λ̃_{n}⟩ of the eigenvalues of (Ã, B̃) such that

lim_{ϵ→0} χ (⟨λ_{i}⟩, ⟨λ̃_{i}⟩) = 0 i = 0, …, n

In other words, no eigenvalues are lost in the process of perturbing the original matrix pair for small enough perturbations.

From [20, Theorem 2.2, page 293] we get

⟨α̃, β̃⟩ = ⟨y^{H} Ãx, y^{H} B̃x⟩ + O(ϵ^{2}).

Actually, in the proof Stewart and Sun show that

χ (⟨α̃, β̃⟩, ⟨y^{H} Ãx, y^{H} B̃x⟩) < O(ϵ^{2}).

For convenience, we will indicate the size of the matrices with a subscript. We will find a matrix pair, (L̃_{(N+1)}, R̃_{(N+1)}) = (L_{(N+1)} − E_{(N+1)}, R_{(N+1)} − F_{(N+1)}) which has an infinite generalized eigenvalue ⟨∞⟩ for each N. Defining ϵ = \sqrt{∥E∥_{2}^{2} + ∥F ∥_{2}^{2}} as the size of the perturbation as in [20], the matrix pair (L_{(N+1)}, R_{(N+1)}) will be an O(ϵ) perturbation of (L̃_{(N+1)}, R̃_{(N+1)}). If ⟨α, β⟩ is the generalized eigenvalue of (L_{(N+1)}, R_{(N+1)}) corresponding to the infinite eigenvalue of (L̃_{(N+1)}, R̃_{(N+1)}) the perturbation Theorem 2.2 of Stewart and Sun gives the asymptotic behavior of ⟨α, β⟩.

Start by defining E_{(N+1)} to be the (N +1)×(N +1) zero matrix so L̃_{(N+1)} = L_{(N+1)}. Define

F_{(N+1)} = (\begin{matrix}0 & ⋯ & 0 & 0 & 0 \\ ⋮ & ... & ⋮ & ⋮ & ⋮ \\ 0 & ⋯ & 0 & 0 & 0 \\ 0 & ⋯ & 0 & R_{N,N−1} & R_{N,N} \end{matrix} )

so R̃_{(N+1)} = R_{(N+1)} − F_{N+1} and all of the entries in the last row of R̃_{(N+1)} are zero for every N. Hence R̃_{(N+1)} is singular for every N.

Now it must be checked that the new pair (L̃_{(N+1)}, R̃_{(N+1)}) satisfies the hypotheses of Theorem 2.2. First,

det (1 · L̃_{(N+1)} − 0 · R̃_{(N+1)}) = det L̃_{(N+1)} = det L_{(N+1)} ≠ 0

and so the matrix pair is regular. Then if x̃ is any null vector of R̃_{(N+1)} consider

0 · L̃_{(N+1)} x̃ − 1 · R̃_{(N+1)} x̃ = 0

so the matrix pair (L̃_{(N+1)}, R̃_{(N+1)}) has an infinite eigenvalue.

Now det (L̃_{(N+1)} − λR̃_{(N+1)}) = 0 is a polynomial of degree N and so has N roots (possibly complex), say λ_{i} for i = 0, …, N−1. Therefore, (L̃_{(N+1)}, R̃_{(N+1)}) has N eigenvalues of the form ⟨λ_{i}, 1⟩ for i = 0, …, N −1. None of these generalized eigenvalues are infinite eigenvalues. The regular matrix pair (L̃_{(N+1)}, R̃_{(N+1)}) has exactly N +1 eigenvalues, counting multiplicity, hence the infinite eigenvalue must be simple.

To investigate the asymptotic behavior of the corresponding generalized eigenvalue ⟨α, β⟩ of (L_{(N+1)}, R_{(N+1)}) we must verify that ϵ → 0 as N → ∞ where ϵ^{2} = ∥F_{N+1}∥_{2}^{2}. Now ϵ^{2} is the largest eigenvalue of F_{N+1}^{T}F_{N+1}. This matrix has two eigenvalues: 0 of multiplicity N and R_{N,N−1}^{2} +R_{N,N}^{2} of multiplicity 1. Hence

∥F ∥_{2}^{2} = R_{N,N−1}^{2} + R_{N,N}^{2} = O(N^{6}/2^{4N}) as N → ∞

by using the expressions for the entries R_{N,N−1} and R_{N,N} in Appendix D. Using the continuity of generalized eigenvalues ([20, Theorem 2.1, page 291]) ⟨α, β⟩ does indeed approach ⟨∞⟩ as N → ∞.

Note that the left and right eigenvectors ỹ and x̃ corresponding to ⟨1, 0⟩ are the left and right null vectors of R̃_{(N+1)}. The left eigenvector is ỹ^{T} = (0, 0, ⋯, 0, 1). The right eigenvector, x̃^{T} = (x_{0}, …, x_{N}) solves

R̃_{(N+1)}x̃ = 0

where there is one more variable than equations because the last row of R̃_{(N+1)} is zero. With this in mind choose x_{N} to be the independent value. Rewriting the matrix equation,

R_{(N)} ( x_{0} ; ⋮ ; x_{N−1}) = ( −R_{0,N}x_{N} ; ⋮ ; −R_{N−1,N}x_{N} ) .

The solution for x_{N−1} is x_{N−1} = −(2N^{2} + N + 2)x_{N}/(4N + 4), see [6, Section 4.3]. We can normalize the null vector by taking x_{N} = 1/L_{N,N} which is sufficient to make the generalized eigenvalue ⟨ỹ^{H} L̃_{(N+1)}x̃, ỹ^{H} R̃_{(N+1)}x̃⟩ = ⟨1, 0⟩.

Using the known generalized eigenvalue ⟨ỹ^{H} L̃_{(N+1)}x̃, ỹ^{H} R̃_{(N+1)}x̃⟩ = ⟨1, 0⟩ of the pair (L̃_{(N+1)}, R̃_{(N+1)}), the perturbations −E_{(N+1)} = 0 and −F_{(N+1)} with nonzero entries given in Appendix D, and the eigenvectors ỹ and x̃ in the approximation formula on page 293 of [20] we obtain

χ (⟨1, 0⟩, ⟨α_{(N+1)}, β_{(N+1)}⟩) ≅ |ỹ^{H}F x̃| = \frac{2N − 1}{8(2N + 3)(2N + 1)(N + 1)}.

Theorem 4 There is a generalized eigenvalue of the Chebyshev-tau approximation system (22) to the model problem (1) which is O(N^{−2})-close to ⟨∞⟩.

6.4 Infinite eigenvalues in the Legendre-tau method (ν = 1/2)

The first step in proving the existence of spurious eigenvalues for the Gegenbauertau method is to examine the case of ν = 1/2. Recall that ν = 1/2 corresponds to using Legendre polynomials. This section will show that for ν = 1/2 the Gegenbauer-tau method yields an infinite generalized eigenvalue.

Theorem 1 Let N > 0 be an integer. The matrix pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) is regular and has a simple infinite generalized eigenvalue.

Proof: As in (23) in Section 6.2 we know that L_{(N+1)}^{(1/2)} is upper triangular and nonsingular. Recall from Section 6.2 that if one of the matrices in the matrix pair in non-singular then the matrix pair is regular.

We also saw in Section 6.1 that R_{(N+1)}^{(1/2)} is upper Hessenberg. The elements of the subdiagonal are

R_{i,i+1}^{(1/2)} = \sqrt{π} (2i + 2)(2i + 1)2^{−2i} \frac{(2i + 1/2) Γ(2i + 1)}{Γ(2i + \frac{3}{2})}

which are nonzero for 1 ≤ i ≤ N. Therefore, rows 1 through N of R_{(N+1)}^{(1/2)} are

linearly independent and so R_{(N+1)}^{(1/2)} has rank of at least N.

The elements of the first row of R_{(N+1)}^{(1/2)} are (see Appendix D)

R_{0,j}^{(1/2)} = \sqrt{π}(2j)(2j − 1)2^{−2j+1}\binom{2j − 2}{j − 1}\frac{Γ(j)}{Γ(j + 1/2)}

− 2\sqrt{π}(2j + 2)(2j + 1)2^{−2j−1}\binom{2j}{j}\frac{Γ(j + 1)}{Γ(j + \frac{3}{2})}

+ \sqrt{π}(2j + 4)(2j + 3)2^{−2j−3}\binom{2j + 2}{j + 1}\frac{Γ(j + 2)}{Γ(j + \frac{5}{2})}.

Using

\frac{Γ(x + 1/2)}{Γ(2x)} = \frac{π}{2^{2x−1}Γ(x)}

to eliminate the fractions from the Γ functions in the denominator of each term gives R_{0,j}^{(1/2)} = 0. Therefore, the first row of R_{(N+1)}^{(1/2)} is a row of zeroes and

R_{(N+1)}^{(1/2)} is singular and has rank N.

The matrix pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) has at least one infinite generalized eigenvalue ⟨1, 0⟩. It is not hard to see that this is a simple generalized eigenvalue. The associated right eigenvector of this infinite generalized eigenvalue must be a right null vector for R_{(N+1)}^{(1/2)} and the null space of R_{(N+1)}^{(1/2)} is of dimension one.

Therefore, there is only one linearly independent right null vector for R_{(N+1)}^{(1/2)} and so the infinite generalized eigenvalue must be simple.

□

In Section 6.5 we will perform a perturbation analysis on the Legendre matrices and so we will need the right and left eigenvectors of the matrix pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) associated with the infinite generalized eigenvalue.

Theorem 2 The left eigenvector of (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) is

y_{(N+1)} = (1, 0, …, 0)

where there are N + 1 entries in the vector. The right eigenvector, x_{(N+1)}, is

x_{j} = \frac{(−1)^{N−j}}{2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}

for 0 ≤ j ≤ N

The proof is technical, and can be found in [6, Section 5.3]. The reader can verify the formulas in specific cases.

Now that we know that the matrix pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) is regular and has a simple infinite generalized eigenvalue we can use this to show that for a range of ν’s around ν = 1/2 the Gegenbauer tau approximation to the solution of the model problem will give a spurious eigenvalue.

6.5 Perturbation analysis

With the technical lemmas from Appendix E we can prove that a spurious eigenvalue will arise from the application of the Gegenbauer-tau approximation method to the model problem for 0 ≤ ν < \frac{1}{2} and ν > \frac{1}{2}. Recall from Theorem 1

that the matrix pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) which arises from the application of

the Legendre-tau method (ν = \frac{1}{2}) to the model problem has a simple infinite generalized eigenvalue. Also, recall from Section 6.3 that if another matrix pair (A, B) of (N + 1) × (N + 1) matrices is “close” to (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) then there will be a generalized eigenvalue ⟨α_{(N+1)}, β_{(N+1)}⟩ of (A, B) that corresponds to the infinite generalized eigenvalue of (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}). The corresponding eigenvalue is approximately

⟨α_{(N+1)}, β_{(N+1)}⟩ = ⟨y^{H}Ax, y^{H}Bx⟩ + O(ϵ^{2})

where y and x are the left and right eigenvectors of the pair (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}) and ϵ is

ϵ = \sqrt{||A − L_{(N+1)}^{(1/2)}||^{2} + ||B − R_{(N+1)}^{(1/2)}||^{2}}.

Let 0 ≤ ν < \frac{1}{2} or ν > \frac{1}{2} be given. We then have the following theorem proving the existence of a spurious eigenvalue for these values of ν. This follows because under the norm defined in Appendix E and using Lemma 6 we can make, for sufficiently large N, R_{(N+1)}^{(ν)} as close to R_{(N+1)}^{(1/2)} as we want and similarly for

L_{(N+1)}^{(ν)} and L_{(N+1)}^{(1/2)}.

Theorem 3 Let η > 0 and ν ≥ 0, ν ≠ 1/2 be given. Let ⟨α^{(ν)}, β^{(ν)}⟩ be the generalized eigenvalue of (L_{(N+1)}^{(ν)}, R_{(N+1)}^{(ν)}) that corresponds to the infinite gen-

eralized eigenvalue langle∞⟩ = ⟨1, 0⟩ of (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}). Then there exists an M such that for N ≥ M,

⟨α^{(ν)}, β^{(ν)}⟩ = ⟨1, \frac{2(2ν − 1) (2ν − 3)}{(2N + 3)(2N + 2) (2N − 2ν + 3) (2N + 2ν + 2)}⟩ + O(ϵ^{2})

where

ϵ^{2} < 16 (A_{N,ν} + A_{N, \frac{1}{2}} ) N^{2} (51N^{2} + 7) e^{2J_{min}−2N}(1 + η)

J_{min} is defined in Appendix E.

Proof: First, establish M and the upper bound on ϵ^{2}. For the given η we get an M from Lemma 6 such that

||L_{(N+1)}^{(ν)}|| < 816A_{N,ν} N^{4}e^{2J_{min}−2N}(1 + η)
||L_{(N+1)}^{(1/2)}|| < 816A_{N,1/2} N^{4}e^{2J_{min}−2N}(1 + η)

and

||R_{(N+1)}^{(ν)}|| < 112A_{N,ν} N^{2}e^{2J_{min}−2N}(1 + η)
||R_{(N+1)}^{(1/2)}|| < 112A_{N,1/2} N^{2}e^{2J_{min}−2N}(1 + η)

for N ≥ M. Hence,

ϵ^{2} = ||L_{(N+1)}^{(ν)} − L_{(N+1)}^{(1/2)}||^{2} + ||R_{(N+1)}^{(ν)} − R_{(N+1)}^{(1/2)}||^{2}

≤ ||L_{(N+1)}^{(ν)}|| + ||L_{(N+1)}^{(1/2)}|| + ||R_{(N+1)}^{(ν)}|| + ||R_{(N+1)}^{(1/2)}||

< 816 (A_{N,ν} + A_{N, \frac{1}{2}} ) N^{4}e^{2J_{min}−2N}(1 + η) + 112 (A_{N,ν} + A_{N, \frac{1}{2}} ) N^{2}e^{2J_{min}−2N}(1 + η)

= 16 (A_{N,ν} + A_{N, \frac{1}{2}} ) N^{2} (51N^{2} + 7) e^{2J_{min}−2N}(1 + η).

Notice that ϵ^{2} → 0 at about the same rate as N^{5}e^{−2N} → 0.

Now we can calculate the approximation to ⟨α_{(N+1)}, β_{(N+1)}⟩. Recall that

⟨α^{(ν)}, β^{(ν)}⟩ = ⟨y^{H}L_{(N+1)}^{(ν)}x, y^{H}R_{(N+1)}^{(ν)}x⟩ + O(ϵ^{2})
≈ ⟨y^{H}L_{(N+1)}^{(ν)}x, y^{H}R_{(N+1)}^{(ν)}x⟩

where the right and left eigenvectors of (L_{(N+1)}^{(1/2)}, R_{(N+1)}^{(1/2)}), x and y respectively, are given by Theorem 2. Using the fact that only the first entry in the left eigenvector is non-zero we see that,

y^{H}L_{(N+1)}^{(ν)}x = ∑_{j=0}^{N} L_{0,j}^{(ν)}x_{j}

y^{H}R_{(N+1)}^{(ν)}x = ∑_{j=0}^{N} R_{0,j}^{(ν)}x_{j}.

Using Lemmas 7 and 8 from Appendix E to evaluate the summations and simplifying gives (see [6] for details)

y^{H}L_{(N+1)}^{(ν)}x = \frac{AB_{ν} (\frac{5}{2} − ν)_{N} (2N + 4)(2N + 3)(2N + 2)(2N + 1)}{Γ(ν + 1) (ν + 1)_{N}}
y^{H}R_{(N+1)}^{(ν)}x = −\frac{A B_{ν} (2ν − 1) (\frac{3}{2} − ν)_{N} (N + 2)(2N + 1)}{Γ(ν + 2) (ν + 2)_{N}}

Recall from the definition of generalized eigenvalues that if ⟨α_{(N+1)}, β_{(N+1)}⟩ is a generalized eigenvalue for a matrix pair then so is

⟨1, \frac{β}{α}⟩

provided α ≠ 0. Therefore we get

⟨α^{(ν)}, β^{(ν)}⟩ ≈ ⟨y^{H}L_{(N+1)}^{(ν)}x, y^{H}R_{(N+1)}^{(ν)}x⟩ = ⟨1, \frac{y^{H}R_{(N+1)}^{(ν)}x}{y^{H}L_{(N+1)}^{(ν)}x}⟩

Substituting in and canceling wherever possible we get,

⟨α^{(ν)}, β^{(ν)}⟩ ≈ ⟨1, −\frac{(2ν − 1) (\frac{3}{2} − ν)_{N} Γ(ν + 1) (ν + 1)_{N}}{2 (\frac{5}{2} − ν)_{N} (2N + 3)(2N + 2) Γ(ν + 2) (ν + 2)_{N}} ⟩ (24)

Notice that

\frac{Γ(ν + 1) (ν + 1)_{N}}{Γ(ν + 2) (ν + 2)_{N}} = \frac{1}{(N + ν + 1)}

Also,

\begin{matrix}(\frac{3}{2} − ν) & N \\ (\frac{5}{2} − ν) & N \end{matrix} = \frac{(\frac{3}{2} − ν)}{(\frac{3}{2} − ν + N)}

With these

⟨α^{(ν)}, β^{(ν)}⟩ ≈= ⟨1, \frac{2(2ν − 1) (2ν − 3)}{(2N + 3)(2N + 2) (2N − 2ν + 3) (2N + 2ν + 2)}⟩.

□

This theorem leads to the following results.

Theorem 4 The Chebyshev-tau method applied to the model problem, (1), will yield a spurious eigenvalue.

Theorem 5 Let 0 ≤ ν ≤ 1 be given. The Gegenbauer-tau approximation to the model problem (1) will yield

1. a positive spurious eigenvalue for 0 ≤ ν < 1/2,

2. an infinite generalized eigenvalue for ν = 1/2,

3. a negative spurious eigenvalue for 1/2 < ν ≤ 1

for all truncation orders.

7 Origin of the spurious eigenvalues in the Legendretau method

This section explains how the boundary conditions act together with the Legendretau method to allow the existence of an infinite eigenvalue. The Gegenbauer-tau method for values of ν near the Legendre case ν = 1/2 will then have an eigenvalue that is a perturbation of the infinite eigenvalue. For 0 ≤ ν < 1/2 we have already shown that the perturbation moves the infinite eigenvalue into the regime where it is positive. For 1/2 ≤ ν < 1 we have already shown that the perturbation moves the eigenvalue into the regime where it is negative. The origin of the infinite eigenvalue in the Legendre-tau method together with the perturbation theory for generalized eigenvalues explains the origin of the spurious eigenvalues.

We need to show that the Legendre-tau method of all orders applied to the model problem (1) yields an infinite generalized eigenvalue. For a given order N + 4 an infinite generalized eigenvalue would be the pair ⟨α_{(N+1)}, β_{(N+1)}⟩ = ⟨1, 0⟩ for which there is a polynomial of degree N + 4 whose residual for the model problem

βD^{4}u = αD^{2}u
u(−1) = u(1) = u′(−1) = u′(1) = 0

is orthogonal to the Legendre polynomials P_{0}(x), … P_{N}(x). That is, the goal is to find a polynomial solution, expanded in Legendre polynomials, for the problem

(0 · D^{4} − 1 · D^{2}) ( ∑_{j=0}^{N+4} a_{j}P_{j}(x)) = τP_{N+2}(x)

together with the boundary conditions.

As explained previously in Section 6.1, because we seek polynomial solutions we may incorporate the boundary conditions directly into the equation. That is, we want to show that there is a solution of the problem

−D^{2} (∑_{j=0}^{N} a_{j}(1 − x^{2})^{2}P_{j}(x)) = τP_{N+2}(x).

We can rescale and slightly simplify the problem by dividing through by −τ and incorporating it into the unknown coefficients a_{j} to obtain

D^{2} (∑_{j=0}^{N} a_{j}(1 − x^{2})^{2}P_{j}(x)) = P_{N+2}(x).  (25)

To simplify and eliminate the consideration of cases, consider only when N is even, say N = 2M. We will solve this equation by expanding it in the linearly independent even Legendre polynomials as basis functions. Equation (25) can be rewritten as

∑_{j=0}^{M} a_{j}D^{2}[(1 − x^{2})^{2}P_{2j}(x)] = P_{N+2}(x).

Now expand the polynomial D^{2}[(1 − x^{2})^{2}P_{j}(x)] in even Legendre polynomials as

D^{2}[(1 − x^{2})^{2}P_{2j}(x)] = ∑_{k=0}^{j+1} K_{j,k}P_{2k}(x).

From the orthogonality relation, the coefficients K_{j,k} for k = 0, …, j + 1 are

K_{j,k} = \frac{2k + 1}{2} ∫_{−1}^{+1} D^{2}[(1 − x^{2})^{2}P_{2j}(x)]P_{2k}(x) dx.

Integrating by parts twice and rearranging, the coefficients are

K_{j,k} = \frac{2k + 1}{2} ∫_{−1}^{+1} P_{2j}(x)(1 − x^{2})^{2}P′′_{2k}(x) dx.

For a given j, by orthogonality only the even Legendre polynomials P_{2k}(x) of degrees 2j + 2, 2j and 2j − 2 can contribute non-zero terms. That is,

D^{2}[(1 − x^{2})^{2}P_{2j}(x)] = K_{j,j−1}P_{2j−2}(x) + K_{j,j}P_{2j}(x) + K_{j,j+1}P_{2j+2}(x)

Using this information, we must show the M+2 equations in the M+1 unknowns a_{0} … a_{M} generated from (25) by equating coefficients of P_{k}(x), k = 0, …, M +1 on left and right are consistent.

The first of these equations, equating coefficients of P_{0}(x), would have all coefficients on the left 0, since

K_{j,0} = \frac{1}{2} ∫_{−1}^{+1} P_{2j}(x)(1 − x^{2})^{2}P′′_{0}(x) dx = 0.

The equation from the coefficients for P_{2k}, k = 1, …, M − 1 will involve a_{k−1}, a_{k}, a_{k+1}. Finally, the equation from the coefficients of P_{2M+2} will only involve a_{M}. For example, when M = 4, so we are expanding the solution in P_{0}(x), P_{2}(x), P_{4}(x), P_{6}(x), P_{8}(x) and the residual is P_{10}(x), the matrix is

(\begin{matrix}0 & 0 & 0 & 0 & 0 \\ 8 & −\frac{16}{7} & \frac{8}{21} & 0 & 0 \\ 0 & \frac{72}{7} & −\frac{720}{77} & \frac{360}{143} & 0 \\ 0 & 0 & \frac{560}{33} & −\frac{224}{11} & \frac{112}{17} \\ 0 & 0 & 0 & \frac{336}{13} & −\frac{672}{19} \\ 0 & 0 & 0 & 0 & \frac{11880}{323} \end{matrix} ) .

Clearly, the system is tridiagonal. The rank of the coefficient matrix is M + 1, and the rank of the augmented matrix is also M + 1. The system is consistent and there is a solution. The implication is that the Legendre-tau method allows an infinite eigenvalue.

The Legendre polynomials can be expressed by the Rodriguez formula [1]

P_{j}(x) = \frac{(−1)^{j}}{j! 2^{j}} \frac{d^{j}}{dx^{j}} (1 − x^{2})^{j}.

The Rodriguez formula for the Gegenbauer polynomials is, for comparison,

G_{j}^{(ν)}(x) = \frac{1}{e_{j}(1 − x^{2})^{ν−1/2}} \frac{d^{j}}{dx^{j}} [(1 − x^{2})^{ν−1/2}(1 − x^{2})^{j}]

where the normalizing constant e_{j} is

e_{j} = (−1)^{j}2^{j}j! \frac{Γ(2ν)Γ(ν + j + 1/2)}{Γ(ν + 1/2)Γ(j + 2ν)}.

See [1, Equation 22.11.2, page, 785].

Using the Rodriguez formula in (25), we want to show that

\frac{d^{2}}{dx^{2}} (∑_{j=0}^{N} a_{j}(1 − x^{2})^{2} \frac{d^{j}}{dx^{j}} \frac{(−1)^{j}}{j! 2^{j}} (1 − x^{2})^{j}) = \frac{d^{N+2}}{dx^{N+2}} \frac{(−1)^{N+2}}{(N + 2)! 2^{N+2}} (1 − x^{2})^{N+2}. (26)

has a solution. Of course, from the argument in the previous paragraph, we know that (26) has a solution. The expression in terms of the Rodriguez expansion makes it plausible that a solution exists.

The expansion in terms of the Rodriguez formula explains why the Legendretau spectral method has an infinite eigenvalue for the boundary-value problem (1). The boundary conditions, when incorporated into a polynomial solution, has exactly the form of the basis of the polynomial. The similarity of the boundary conditions to the form of the polynomial solution permits a solution where none would ordinarily be expected.

8 Discussion of ad hoc Methods for Eliminating the Spurious Eigenvalue

8.1 Factoring the differential equation

Gottlieb and Orszag suggest a simple method for eliminating the spurious unstable (positive) eigenvalues in the model problem, and also in the Orr-Sommerfeld equation, [11, page 145]. Gardner, Trogdon and Douglass generalized the method to fourth-order eigenvalue problems [9]. For the model problem, the suggestion is to solve a factored problem

\begin{aligned}v′′ & = sv \\ u′′ & = v \\ u(−1) & = u′(−1) = u(1) = u′(1) = 0 \end{aligned}  (27)

for the functions u and v. The Chebyshev-tau method sets u(x) = ∑_{k=0}^{N+2} a_{k}T_{k}(x) and v(x) = ∑_{k=0}^{N+2} b_{k}T_{k}(x). Gottlieb and Orszag prove that the solution of (27) with the Chebyshev-tau method has only negative eigenvalues, at least in the special case with all odd coefficients a_{2n+1} = b_{2n+1} = 0. Gardner et al. show with examples that the factorization method has no positive eigenvalues even in more complicated problems, see [9].

Following [9], since the boundary conditions apply to the function u, therefore involving only the coefficients a_{n}, we eliminate the b-coefficients from the resulting pair of matrix equations. All eigenvalues of the resulting matrix equation are negative, and the small magnitude eigenvalues approach the eigenvalues (2) of the original problem (1). Although Gardner et al. did not do so, we find the dominant eigenvalue by the power method. The results are in Table 3.

The absence of positive eigenvalues might suggest that there are no spurious eigenvalues, but the sequence of dominant eigenvalues in the table shows that the dominant eigenvalue (the eigenvalue largest in magnitude) is rapidly growing in magnitude. The dominant eigenvalue grows approximately as N^{3.4}.

8.2 Elimination of spurious eigenvalues by the ad hoc method of McFadden, Murray, and Boisvert

In [16] McFadden, Murray, and Boisvert propose an extremely simple modification of the Chebyshev-tau method to eliminate the spurious eigenvalues. For the model problem (1) instead of using the approximating equation (12) from u(x, N) = ∑_{k=0}^{N+4} a_{k}T_{k}(x), McFadden, Murray and Boisvert propose instead using

a_{k}^{(4)} − sā_{k}^{(2)} = 0  k = 0, 1, 2, …, N  (28)

where they define

c_{k}ā_{k}^{(2)} = ∑_{\begin{matrix}p=k+2 \\ p+k even \end{matrix}}^{N+2} p(p^{2} − k^{2})a_{p}

N Dominant Eigenvalue
3  −26.6667
5  −110.3716
7  −305.9584
9  −686.5736
11  −1343.2847
13  −2385.2583
15  −3939.8243
17  −6152.4924
21  −13225.0986
25  −25130.8741
29  −43688.6149
33  −71008.1217
37  −109490.2008
41  −161826.6649

Table 3: The dominant eigenvalue computed using the factorization method of Gardner, et al.

and c_{k} is the normalizing factor arising from the orthogonality of the Chebyshev polynomials, compare (31) in Appendix A.

Note that ā_{k}^{(2)} is obtained from a_{k}^{(2)} in (31) by simply dropping the last two terms. Thus, the modified matrix eigenvalue problem corresponding to (28) is obtained from the original eigenvalue problem (13) by setting the last two columns of B to zero. Using this modification of McFadden et al., all of the eigenvalues are non-positive. For example, using N = 5 the eigenvalues are −9.8701, −20.1942, −40.5258, −62.6474, −201.604, −305.958. The absence of positive eigenvalues might suggest that there are no spurious eigenvalues in the suggested method, but closer inspection indicates that the dominant eigenvalue (the eigenvalue largest in magnitude) is rapidly growing in magnitude. In fact, Table 4 indicates that the magnitude of the dominant eigenvalue is growing as approximately N^{2.8}. Note that the values of the dominant eigenvalue for N in the table are precisely the same as the values for N + 2 in Table 3 from the factorization method of Gardner et al.. This is not accidental, since for the model problem (1) for given N the equations (28) given by McFadden et al. in [16] are exactly the same as the single reduced matrix equation (3.9b) on page 147 of Gardner et al. in [9]. That is, for the model problem the ad hoc method of McFadden et al. is the same as the factorization method of Gardner et al.

8.3 Elimination of spurious eigenvalues by the pseudospectral method of Huang and Sloan

In [13], Huang and Sloan present a pseudospectral method for eigenvalue problems. They present several numerical examples that show that the method produces only negative eigenvalues and that it is comparable in accuracy to the

N Dominant Eigenvalue
5  −305.9584
9  −1343.2847
13  −3939.8243
17  −9186.9569
21  −18466.9859
25  −33453.2068
29  −56109.9177
33  −88692.4219
37  −133747.0290
41  −177784.9199

Table 4: The dominant eigenvalue computed using the method of McFadden et al.

Chebyshev-tau method. They use interpolating polynomials u^{N}(x) and ũ^{N}(x) satisfying the boundary conditions u^{N}(±1) = (u^{N})′(±1) = 0 and ũ^{N}(±1) = 0 respectively. For the collocation points

−1 = x_{1} < x_{2} < ⋯ < x_{N−1} < x_{N} = 1

the collocation equations are

(u^{N})^{(4)}(x_{j}) = s(ũ^{N})′′(x_{j}),  2 ≤ j ≤ N − 1.  (29)

The key feature of the scheme is that it uses two distinct interpolation polynomials. The degree of the polynomial which approximates the second-order derivative term on the right-hand side of (29) is lower by two than the degree of the polynomial that approximates the fourth-order derivative on the left. This has characteristics in common with the approach of McFadden et al. described above.

Using the pseudospectral method of Huang and Sloan, all of the eigenvalues are nonpositive. For example, using N = 5 the eigenvalues are −9.8531,−80.0000, −146.1469. The absence of positive eigenvalues might again suggest that there are no spurious eigenvalues in the suggested method, but closer inspection indicates that the dominant eigenvalue (the eigenvalue largest in magnitude) is rapidly growing in magnitude. In fact, Table 5 indicates that the dominant eigenvalue is growing as approximately N^{4.7}.

Even with the preconditioning matrix suggested by Huang and Sloan, [13, equation 2.15, and the Appendix], the matrix-eigenvalue problem (29) is numerically very sensitive and requires very high precision to obtain convergence of the eigenvalues to the known eigenvalues. The sensitivity accounts for our inability to extend the table to higher orders of truncation as was done for the other spectral methods.

N Dominant eigenvalue
5  −146.1469
9  −2632.3385
13  −14771.6090
17  −49706.5247
21  −126419.7558
25  −269734.6733

Table 5: The dominant eigenvalue computed with the pseudospectral method of Huang and Sloan

9 Conclusions

In prior sections, we have an explanation for the origin and a description of the nature of the spurious eigenvalues in spectral methods for a differential eigenvalue problem. In particular, we have a explanation for the widely reported positive spurious eigenvalues in the Chebyshev-tau method. The analysis of the spurious eigenvalues was in terms of a simplified model problem; a problem that is simple enough to have an explicit solution to validate the spectral methods, yet still contains the essence of more general problems. Many other investigators have used the model problem as an example of the use of spectral methods. Finally, the model problem occurs naturally in the investigation of some fluid dynamics problems. From that analysis we observe that:

1. The model problem factors into odd and even problems each generating a spurious eigenvalue. The factoring simplifies the analysis, but is not essential to the analysis.

2. The popular Chebyshev-tau and Legendre-tau methods are each instances in a family of spectral methods using the Gegenbauer polynomials.

3. One member of the family of Gegenbauer-tau methods, namely the Legendretau method when ν = 1/2, applied to the model problem has an infinite eigenvalue in the sense of generalized eigenvalue theory.

4. Positive spurious eigenvalues occurring for 0 ≤ ν < 1/2 are approximations to infinite eigenvalues in a generalized eigenvalue problem for ν = 1/2.

5. Large magnitude negative eigenvalues occur when 1/2 ≤ ν ≤ 1, prompting us to enlarge the definition of spurious eigenvalues to include all largemagnitude eigenvalues that are perturbations of an infinite eigenvalue, regardless of sign.

6. The boundary conditions, when incorporated into a polynomial solution, match the form of the Legendre polynomials. The similarity of the boundary conditions to the form of the Legendre polynomials permits the infinite eigenvalue.

Figure 2: Schematic diagram of spectral methods and perturbations of infinite eigenvalues

7. Spurious eigenvalues from the Gegenbauer-tau method grow at least as

fast as

\frac{(2N)(2N + 2ν + 1)(2N + 2ν + 2)(2N + 3)}{1 − 4ν^{2}}

in the truncation order N. In particular, the spurious eigenvalues from the Chebyshev-tau method are larger than (2N)(2N +1)(2N +2)(2N +3).

8. Ad hoc methods intended to eliminate positive eigenvalues but compute large magnitude eigenvalues nonetheless.

For a schematic diagram of the results, consider the Legendre-tau method on the model problem with infinite eigenvalues (along with good approximations to the true eigenvalues) as a point in a space of approximation methods. Other approximation methods are perturbations away from this point. These perturbations change the infinite eigenvalue by perturbing it away from infinity into the positive eigenvalue regime. Some other spectral methods that are seldom used perturb the eigenvalue into the negative eigenvalue regime. Ad hoc methods change either the problem or the method for the Chebyshev-tau method, perturbing the large eigenvalues into the negative eigenvalue regime. We intend further research to make this precise with perturbation results from generalized eigenvalue theory.

Figure 3: Schematic diagram of differential eigenvalue problems and perturbations of infinite eigenvalues

The results suggest another way of visualizing the results is to think of perturbing the differential eigenvalue problem from the model problem. Other more complicated differential eigenvalue problems of a similar form generate positive spurious eigenvalues from application of the Chebyshev-tau method, [9]. One can think of this as being a sequence of perturbations from the Legendretau method applied to the model problem. First the change to the Chebyshevtau method perturbs the infinite eigenvalue into the positive eigenvalue regime. Then the approximation equations change to accommodate the change in the differential equation, changing the eigenvalue approximations, but keeping the positive spurious eigenvalue in the positive eigenvalue regime. We intend further research to also make this precise with perturbation results from generalized eigenvalue theory.

Finally, we observe that a reliable diagnostic tool for use in showing that a problem has eigenvalues that are approximations to infinite eigenvalues (i.e., are spurious) is to observe the growth rate for dominant eigenvalues computed with truncation order, N in the spectral series. If the magnitude of the eigenvalue grows like N^{n}, n ≈ 4, it is probably spurious in the sense we have explained in this paper.

A Formulae for derivatives of Chebyshev polynomials and their application to the model problem

The coefficients of the expansions of the second- and fourth-order derivatives in the Chebyshev-tau expansion (12) are (cf. [11, 9])

c_{k}a_{k}^{(4)} = \frac{1}{24} ∑_{\begin{matrix}p=k+4 \\ p+k even \end{matrix}}^{N+4} p[p^{2}(p^{2} − 4)^{2} − 3p^{4}k^{2} + 3p^{2}k^{4} − k^{2}(k^{2} − 4)^{2}]a_{p} (30)
c_{k}a_{k}^{(2)} = ∑_{\begin{matrix}p=k+2 \\ p+k even \end{matrix}}^{N+4} p(p^{2} − k^{2})a_{p}  (31)

where c_{k} is a normalizing factor arising from the orthogonality of the Chebyshev polynomials:

c_{k} = \begin{cases}2,  if k = 0 \\ 1,  if k > 0 \end{cases} .

The boundary conditions from (1) become

∑_{j=0}^{N+4}(−1)^{j}a_{j} = ∑_{j=0}^{N+4} a_{j} = ∑_{j=0}^{N+4}(−1)^{j+1}j^{2}a_{j} = ∑_{j=0}^{N+4} j^{2}a_{j} = 0.

The system of equations in matrix-vector form is

Aa = sBa

where A has elements

A_{ni} = \frac{χ_{4}(n, i)}{24c_{n}}i[i^{2}(i^{2} − 4)^{2} − 3i^{4}n^{2} + 3i^{2}n^{4} − n^{2}(n^{2} − 4)^{2}]

n = 0, …, N; i = 0, …, N + 4 (32)

where

χ_{M}(n, i) = \begin{cases}0, & if i < n + M or n + i + M odd, \\ 1, & if i ≥ n + M and n + i + M even \end{cases}

and

A_{N+1,i} = (−1)^{i}, A_{N+2,i} = 1, A_{N+3,i} = (−1)^{i+1}i^{2}, A_{N+4,i} = i^{2}

i = 0, …, N + 4. (33)

On the right side, B has elements

B_{ni} = \frac{χ_{2}(n, i)}{c_{n}}i(i^{2} − n^{2})  n = 0, …, N; i = 0, …, N + 4
B_{ni} = 0  n = N + 1, …, N + 4; i = 0, …, N + 4.

B Gegenbauer Polynomials

The Gegenbauer polynomials (sometimes also called the ultraspherical polynomials) are a parametrized family of orthogonal polynomials. The Gegenbauer polynomials contain the Chebyshev polynomials of first and second kind and Legendre polynomials as special cases. The Gegenbauer polynomials {G_{n}^{ν}(x)}_{n=0}^{∞} consist of polynomials G_{n}^{ν}(x) of degree n orthogonal on [−1, 1] with respect to the weight function w_{ν}(x) = (1 − x^{2})^{ν−1/2} where ν > −1/2. Normalize the ultraspherical polynomials so that G_{n}^{ν}(1) = 1, then lim_{νto0} G_{n}^{0}(x) = T_{n}(x), the Chebyshev polynomials of the first kind; G_{n}^{1/2}(x) = P_{n}(x), the Legendre polynomials; and G_{n}^{1}(x) = U_{n}(x), the Chebyshev polynomials of the second kind.

The Gegenbauer polynomials of odd index n are odd polynomials for all ν. Likewise, the Gegenbauer polynomials of even index n are even polynomials for all ν.

The Rodriguez formula for the Gegenbauer polynomials is, [1, Equation 22.11.2, page, 785],

G_{j}^{(ν)}(x) = \frac{1}{e_{j}(1 − x^{2})^{ν−1/2}} \frac{d^{j}}{dx^{j}} [(1 − x^{2})^{ν−1/2}(1 − x^{2})^{j}]

where the normalizing constant e_{j} is

e_{n} = (−1)^{j}2^{j}j! \frac{Γ(2ν)Γ(ν + j + 1/2)}{Γ(ν + 1/2)Γ(j + 2ν)}.

Let

s_{n,ν}(f; x) = ∑_{k=0}^{n} f_{k}(ν)G_{k}^{ν}(x) n = 0, 1, 2, …

where

f_{k}(ν) = ∫_{−1}^{1} f(x)G_{k}^{ν}(x) \frac{dx}{(1 − x^{2})^{k−1/2}}

be the partial sums of the expansion of f in terms of the G_{k}^{ν}(x). Let

S_{n}(f; ν) = ∥f − s_{n,ν}(f)∥ = max_{−1≤x≤1} |f − s_{n,ν}(f)|

In certain cases, the Chebyshev expansion is optimal among all expansions in terms of Gegenbauer polynomials, [19, 3.6, pages 188-191]. In fact

Theorem 5 If ν > 0, either f_{k}(ν) > 0 for k > n or f_{k}(ν) alternate in sign for k > n, and ∑_{k=0}^{∞} f_{k}(ν) converges, then

S_{n}(f, ν) ≥ S_{n}(f, 0)

C Triangular-Hessenberg form of the monomial-Galerkin-Chebyshev method

Lemma 3 1. The matrix L^{(0)} is upper triangular with entries given by

L_{i,j}^{(0)} = (2j)(2j − 1)(2j − 2)(2j − 3)2^{5−2j}\binom{2j − 4}{j − 2 − i}

− 2(2j + 2)(2j + 1)(2j)(2j − 1)2^{3−2j}\binom{2j − 2}{j − 1 − i}

+ (2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{1−2j}\binom{2j}{j − i}.

2. The matrix R^{(0)} is upper Hessenberg with entries given by

R_{i,j}^{(0)} = (2j)(2j − 1)2^{3−2j}\binom{2j − 2}{j − 1 − i}

− 2(2j + 2)(2j + 1)2^{1−2j}\binom{2j}{j − i}

+ (2j + 4)(2j + 3)2^{−1−2j}\binom{2j + 2}{j + 1 − i}.

Proof: On the right side of the matrix generalized eigenvalue problem is an (N + 1) × (N + 1) matrix R^{(0)} indexed from 0 to N. The entry in the (i, j) position on the right is the inner product of

[(1−x^{2})^{2}x^{2j}]′′ = (2j)(2j −1)x^{2j−2} −2(2j +2)(2j +1)x^{2j} +(2j +4)(2j +3)x^{2j+2}

with T_{2i}(x). By [19, equations (1.147), (1.148), page 54-55], the (i, j) entry of R^{(0)} is

\begin{aligned}R_{i,j}^{(0)} & = (2j)(2j − 1)2^{1−2j+2}\binom{2j − 2}{j − 1 − i} − 2(2j + 2)(2j + 1)2^{1−2j}\binom{2j}{j − i} \\ & + (2j + 4)(2j + 3)2^{1−2j−2}\binom{2j + 2}{j + 1 − i}. \end{aligned}

Interpreting binomial coefficients with negative entries to be 0 as usual recovers the upper Hessenberg nature of the matrix.

The entry in the (i, j) position of L^{(0)} is the inner product of

[(1 − x^{2})^{2}x^{2j}]′′′′ = (2j)(2j − 1)(2j − 2)(2j − 3)x^{2j−4}
− 2(2j + 2)(2j + 1)(2j)(2j − 1)x^{2j−2}
+ (2j + 4)(2j + 3)(2j + 2)(2j + 1)x^{2j}

with T_{2i}(x). The (i, j) entry of L^{(0)} is

L_{i,j}^{(0)} = (2j)(2j − 1)(2j − 2)(2j − 3)2^{1−2j+4}\binom{2j − 4}{j − 2 − i}

− 2(2j + 2)(2j + 1)(2j)(2j − 1)2^{1−2j+2}\binom{2j − 2}{j − 1 − i}

+ (2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{1−2j}\binom{2j}{j − i}.

□

D Triangular-Hessenberg form of the monomial-Galerkin-Gegenbauer method

Let −\frac{1}{2} < ν < ∞, ν ≠ 0 be given. Let 2N or 2N + 1 be the order of approximation in the monomial-Galerkin-Gegenbauer method applied to the model problem (1).

Lemma 4 1. The matrix L^{(ν)} is upper triangular with entries given by

L_{i,j}^{(ν)} = (2j)(2j − 1)(2j − 2)(2j − 3)2^{−2j+4}\binom{2j − 4}{j − i − 2}\frac{(2i + ν) Γ(ν) Γ(j + i − 1)}{Γ(j + i + ν − 1)}

− 2(2j + 2)(2j + 1)(2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − i − 1}\frac{(2i + ν) Γ(ν) Γ(j + i)}{Γ(j + i + ν)}

+ (2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j − i}\frac{(2i + ν) Γ(ν) Γ(j + i + 1)}{Γ(j + i + ν + 1)}

2. The matrix R^{(ν)} is upper Hessenberg with entries given by

R_{i,j}^{(ν)} = (2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − i − 1}\frac{Γ(ν)Γ(j + i)}{Γ(j + i + ν)}

− 2(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j − i}\frac{Γ(ν)Γ(j + i + 1)}{Γ(j + i + ν + 1)}

+ (2j + 4)(2j + 3)2^{−2j−2}\binom{2j + 2}{j − i + 1}\frac{Γ(ν)Γ(i + i + 2)}{Γ(j + i + ν + 2)}

Proof: On the right side of the matrix generalized eigenvalue problem is an (N + 1) × (N + 1) matrix R^{(ν)} indexed from 0 to N. The entry in the (i, j) position on the right is the inner product of

[(1−x^{2})^{2}x^{2j}]′′ = (2j)(2j −1)x^{2j−2} −2(2j +2)(2j +1)x^{2j} +(2j +4)(2j +3)x^{2j+2} with G_{2i}^{ν}(x). By [12, equation 7.311.2, page 826], the (i, j) entry of R^{(ν)} is

h_{2i}[(2j)(2j − 1)\frac{π 2^{−2j−2ν+3} Γ(2i + 2ν) Γ(2j − 1)}{Γ(ν) Γ(j − i) Γ(j + i + ν) (2i)!}

− 2(2j + 2)(2j + 1)\frac{π 2^{−2j−2ν+1} Γ(2i + 2ν) Γ(2j + 1)}{Γ(ν) Γ(j − i + 1) Γ(j + i + ν + 1) (2i)!}

+ (2j + 4)(2j + 3)\frac{π 2^{−2j−2ν−1} Γ(2i + 2ν) Γ(2j + 3)}{Γ(ν) Γ(j − i + 2) Γ(j + i + ν + 2) (2i)!}]. (34)

Here is the normalizing factor associated with the inner product with G_{2i}^{(ν)}. Cancelling common factors, and simplifying, the inner product is

(2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − i − 1}\frac{(2i + ν) Γ(ν)Γ(j + i)}{Γ(j + i + ν)}

− 2(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j − i}\frac{(2i + ν) Γ(ν)Γ(j + i + 1)}{Γ(j + i + ν + 1)}

+ (2j + 4)(2j + 3)2^{−2j−2}\binom{2j + 2}{j − i + 1}\frac{(2i + ν) Γ(ν)Γ(j + i + 2)}{Γ(j + i + ν + 2)}. (35)

See [6] for more details.

Likewise, the matrix entries for L^{(ν)} are

L_{i,j}^{(ν)} = ⟨\frac{d^{4}}{dx^{4}} [(1 − x^{2})^{2}x^{2j}] , G_{2i}^{ν}(x)⟩ .

Expanding, using [12, equation 7.311.2, page 826], introducing the normalizing factor h_{2i} for the inner product, cancelling common factors, and simplifying we obtain

L_{i,j}^{(ν)} = (2j)(2j−1)(2j−2)(2j−3)2^{−2j+4}\binom{2j − 4}{j − i − 2}\frac{(2i + ν) Γ(ν) Γ(j + i − 1)}{Γ(j + i + ν − 1)}

− 2(2j + 2)(2j + 1)(2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − i − 1}\frac{(2i + ν) Γ(ν) Γ(j + i)}{Γ(j + i + ν)}

+ (2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j − i}\frac{(2i + ν) Γ(ν) Γ(j + i + 1)}{Γ(j + i + ν + 1)}. (36)

E Definition of the matrix norm

This section contains several facts for the generalized eigenvalue analysis. All proofs are technical and are omitted. Detailed proofs are in [6].

The first detail is the definition of the matrix norm defining the size of the perturbation ϵ in Section 6.3. Let A_{(N+1)} be an (N +1)×(N +1) matrix. Then define

||A_{(N+1)}|| = sup_{0≤i,j≤N} e^{i+j−2N} |a_{i,j}|

It is easy to show that this is indeed a norm.

The rest of this section consists of technical lemmas that will help with the analysis in Section 6.5. These lemmas, while important to the analysis, detract from the actual analysis.

The first two lemmas will allow an estimation on ϵ to be made in Section 6.5. The first lemma will gives an estimation on |L_{i,j}^{(ν)}| and |R_{i,j}^{(ν)}| for 0 ≤ i, j ≤ N and 0 ≤ ν ≤ 1.

Lemma 5 Let ν ≥ 0 and N > 0 be given. Then for any integers i and j such that 0 ≤ i, j ≤ N we have,

|L_{i,j}^{(ν)}| < 2A_{i,ν} [(2j)(2j − 1)(2j − 2)(2j − 3)2^{−2j+4}\binom{2j − 4}{j − 2}+

2(2j + 2)(2j + 1)(2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − 1} +

(2j + 4)(2j + 3)(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j}]

|R_{i,j}^{(ν)}| ≤ 2A_{i,ν} [(2j)(2j − 1)2^{−2j+2}\binom{2j − 2}{j − 1}+

2(2j + 2)(2j + 1)2^{−2j}\binom{2j}{j} +

(2j + 4)(2j + 3)2^{−2j−2}\binom{2j + 2}{j + 1}]

where

A_{i,ν} = \begin{cases}1 & ν = 0 \\ (2i + ν) Γ(ν) & ν ≠ 0 \end{cases} .

The next lemma gives an upper bound on the norm of the matrices L_{(N+1)}^{(ν)}

and R_{(N+1)}^{(ν)} under the norm defined earlier in this section.

Lemma 6 Let ν ≥ 0 be given. Then given ϵ > 0 there exists an M such that if N ≥ M we have,

||L_{(N+1)}^{(ν)}|| < 816A_{N,ν} N^{4}e^{2J_{min}−2N}(1 + ϵ)

and

||R_{(N+1)}^{(ν)}|| < 112A_{N,ν} N^{2}e^{2J_{min}−2N}(1 + ϵ)

where J_{min} is defined below.

For a given η > 0 we know that we can find J such that for j ≥ J,

1 − η < \binom{2j}{j} (\frac{2^{2j}}{\sqrt{πj}} )^{−1} < 1 + η

So, choose J_{min} > 2 such that for j ≥ J_{min} we have

1 − ϵ < \binom{2j − 4}{j − 2} (\frac{2^{2j−4}}{\sqrt{π(j − 2)}})^{−1} < 1 + ϵ

1 − ϵ < \binom{2j − 2}{j − 1} (\frac{2^{2j−2}}{\sqrt{π(j − 1)}})^{−1} < 1 + ϵ

1 − ϵ < \binom{2j}{j} (\frac{2^{2j}}{\sqrt{πj}} )^{−1} < 1 + ϵ

1 − ϵ < \binom{2j + 2}{j + 1} (\frac{2^{2j+2}}{\sqrt{π(j + 1)}})^{−1} < 1 + ϵ

Notice that the choice of J_{min} is independent of the value of N.

The last two lemmas in this section will be used in Section 6.5 to estimate perturbed eigenvalues. The first lemma gives the evaluations of several summations.

Lemma 7 Let N > 0 be an integer. Then

1.

∑_{j=0}^{N} [\frac{2^{−2j}(2j)! (−1)^{N−j}}{j! Γ(j + ν + 2) 2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}]

= \frac{A (\frac{3}{2} − ν)_{N}}{Γ(ν + 2) (ν + 2)_{N}}

2.

∑_{j=0}^{N} [\frac{j 2^{−2j}(2j)! (−1)^{N−j}}{j! Γ(j + ν + 2) 2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}]

= \frac{A (N^{2} + \frac{5}{2}N) (\frac{3}{2} − ν)_{N}}{(\frac{3}{2} − ν) Γ(ν + 2) (ν + 2)_{N}}

3.

∑_{j=0}^{N} [\frac{2^{−2j}(2j)! (−1)^{N−j}}{j! Γ(j + ν + 1) 2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}]

= \frac{A (\frac{5}{2} − ν)_{N}}{Γ(ν + 1) (ν + 1)_{N}}

∑_{j=0}^{N} [\frac{j 2^{−2j}(2j)! (−1)^{N−j}}{j! Γ(j + ν + 1) 2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}]

= \frac{A (N^{2} + \frac{5}{2}N) (\frac{5}{2} − ν)_{N}}{(\frac{5}{2} − ν) Γ(ν + 1) (ν + 1)_{N}}

5.

∑_{j=0}^{N} [\frac{j^{2} 2^{−2j}(2j)! (−1)^{N−j}}{j! Γ(j + ν + 1) 2^{N−j} (N − j)!} \frac{∏_{k=0}^{2N−2j−1}(2N − k)}{∏_{k=0}^{N−j−1}(4N + 3 − 2k)}]

= \frac{A (N^{2} + \frac{5}{2}N) (\frac{5}{2} − ν)_{N}}{(\frac{5}{2} − ν) Γ(ν + 1)(ν + 1)_{N}} (1 + \frac{(N − 1) (N + \frac{7}{2})}{(\frac{7}{2} − ν)})

where

A = \frac{(2N)!}{2^{N} N! ∏_{k=0}^{N−1}(4N + 3 − 2k)}

The last lemma in this section gives simplified formulas for R_{0,j}^{(ν)} and L_{0,j}^{(ν)} for ν ≥ 0

Lemma 8 Let N > 0 be an integer, 0 ≤ j ≤ N and ν ≥ 0 then,

R_{0,j}^{(ν)} = B_{ν}(2ν − 1) \frac{2^{−2j} (2j)!}{j! Γ(j + ν + 2)} (((2ν − 3)j − 2))
L_{0,j}^{(ν)} = B_{ν} \frac{2^{−2j} (2j)!}{j! Γ(j + ν + 1)} ((16ν^{2} − 96ν + 140)j^{2} − (16ν^{2} − 100)j + 24)

where

B_{v} = \begin{cases}2 & ν = 0 \\ νΓ(ν) & ν ≠ 0 \end{cases}

References

[1] M. Abramowitz and I. Stegun. Handbook of Mathematical Functions. Dover, New York, 1965.

[2] L. C. Andrews. Special Functions for Engineers and Applied Mathematicians. Macmillan, New York, 1985.

[3] John P. Boyd. Complex coordinate methods for hydrodynamic instabilities and Sturm-Liouville eigenproblems with an infinite singularity. Journal of Computational Physics, 57:454–471, 1985.

[4] John P. Boyd. Chebyshev & Fourier Spectral Methods, volume 49 of Lecture Notes in Engineering. Springer-Verlag, 1989.

[5] B. Brenier, B. Roux, and P. Bontoux. Comparaison des méthodes tau-Chebyshev et Galerkin dans l’étude de stabilité des mouvements de convection naturelle. Problème des valeurs propres parasites. Journal de Méchanique théorique et appliquée, 5:95–119, 1986.

[6] P. Dawkins. Spurious Eigenvalues in Spectral-tau Methods. PhD thesis, University of Nebraska-Lincoln, 1997.

[7] L. Fox. Chebyshev methods for ordinary differential equations. Computer Journal, 4:318, 1962.

[8] L. Fox and I. Parker. Chebyshev Polynomials in Numerical Analysis. Oxford University Press, London, 1968.

[9] David B. Gardner, Steven A. Trogdon, and Rod W. Douglass. A modified tau spectral method that eliminates spurious eigenvalues. Journal of Computational Physics, 80(1):137–167, January 1989.

[10] David R. Gardner. Linear Stability and Bifurcation of Natural Convection Flows in Narrow-Gap Concentric Spherical Annulus Enclosures. PhD thesis, University of Nebraska, 1988.

[11] David Gottlieb and Steven A. Orszag. Numerical Anlaysis of Spectral Methods. SIAM, Philadelphia, 1977.

[12] I. S. Gradshteyn and I. M. Ryzhik. Tables of Integrals, Series, and Products. Academic Press, 4th edition edition, 1969.

[13] Weizhang Huaung and David M. Sloan. The pseudospectral method for solving differential eigenvalue problems. Journal of Computational Physics, pages 399–410, 1994.

[14] C. Lanczos. Trigonometric interpolation of empirical and analytical functions. J. Math. Phys., 7:123–199, 1938.

[15] C. Lanczos. Applied Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1964.

[16] G. B. McFadden, B. T. Murray, and R. F. Boisvert. Elimination of spurious eigenvalues in the Chebyshev tau spectral method. Journal of Computational Physics, 91:228–239, 1990.

[17] Steven A. Orszag. Accurate solution of the Orr-Sommerfeld stability equation. Journal of Fluid Mechanics, 50:689–703, 1971.

[18] E. L. Ortiz. The Tau method. SIAM Journal on Numerical Analysis, 6:480–492, 1969.

[19] Theodore J. Rivlin. Chebyshev Polynomials. John Wiley & Sons, New York, second edition, 1990.

[20] G. W. Stewart and Ji-Guang Sun. Matrix Pertubation Theory. Academic Press, Boston, 1990.

[21] B. Straughan and D. W. Walker. Two very accurate and efficient methods for computing eigenvalues and eigenfunctions in porous convection problems. Journal of Computational Physics, 127:128–141, 1996.

[22] Y. Y. Su and B. Khomami. Numerical solution of eigenvalue problems using spectral techniques. Journal of Computational Physics, 100:297–305, 1992.

[23] Kevin Tebeest. Linear Stability and Bifurcation of Natural Convection Flows within Arbitrary-Gap Spherical Annuli. PhD thesis, University of Nebraska, 1992.

[24] Abdelfattah Zebib. A Chebyshev method for the solution of boundary value problems. Journal of Computational Physics, 53:443–455, 1984.

[25] Abdelfattah Zebib. Removal of spurious modes encountered in solving stability problems by spectral methods. Journal of Computational Physics, 70:521–525, 1987.