It is a matter of evidence that the hydrodynamic stability problems become more complex due to incorporation of technologically important effects. Consequently, the order of the differential equations increases rapidly.

Generally the discretization methods applied to such problems lead to matrices which become less conditioned as both the order of approximation and the order of differentiation increase. This fact is particularly true for spectral methods. Roughly, for $N=2^9$ ($N$ stands for the spectral parameter), the condition number for the second order Chebyshev differentiation matrix is of order $O\left({10}^{18}\right),$ that of the fourth order differentiation matrix attains something of order $O\left({10}^{24}\right)$ and eventually that of the eighth order differentiation matrix raises to something of order $O\left({10}^{30}\right).$

In this context, the aim of this paper is two-fold. First, we introduce a new method to solve numerically even and high order two-point boundary value problem with Dirichlet and hinged boundary conditions.  It reduces the problem to a system of second order differential equations supplied only with Dirichlet boundary conditions. Further, accurate spectral methods are used in order to solve that Dirichlet boundary value problem numerically. Second, the method is applied to solve a hydrodynamic stability problem of Bénard type.

It is well known that, the classical Orr-Sommerfeld eigenvalue problem, which is encountered in the linear hydrodynamic stability analysis of some basic flows, was first solved using a direct (tau) spectral method by S. Orszag in 1971 in his seminal paper [21]. He determined the critical Reynolds number associated with the most unstable shear mode corresponding to Poiseuille flow, i.e., $Re=5772.22.$

Two other important papers, of Dongarra, Straughan and Walker [5] and that of Melenk, Kirchner and Schwab [20], consider in a general mathematical framework spectral methods for hydrodynamic stability problems. They analyze the Orr-Sommerfeld equation supplied with homogeneous boundary conditions which contain derivative up to the first order. The pseudospectra and the numerical range of this Orr-Sommerfeld operator is also investigated by Reddy, Schmid and Henningson in their paper [25].

The high accuracy of tau, Galerkin and collocation spectral methods were subsequently, intensively exploited in solving second order as well as higher order (especially fourth order) differential eigenvalue problems. The main difficulty with the later problems consists in imposing boundary conditions which involve derivatives of order higher than one. From the mathematical point of view this reduces to a Hermite interpolation problem. In  a quite general setting it was solved by Huang and Sloan in their paper [18]. Unfortunately, in practical computations, the resulting trial and test functions lead to tedious computations with significant rounding errors.

Up to our knowledge there are two main strategies to avoid these difficulties. The first one requires a special attention in choosing the spaces of test and trial functions. It is important to design spaces of functions which satisfy as many as possible boundary conditions. In our previous paper [11] we built up a Chebyshev tau method in order to solve a hydrodynamic stability problem connected with the Marangoni-Plateau-Gibbs convection. The eigenvalue problem was a non-standard one, i.e., two out of the four boundary conditions attached to the Orr-Sommerfeld equation contained derivatives of order two and three as well as the eigenparameter. Our numerical results were lately confirmed by Greenberg and Marleta [13]. A more recent paper of Blyth and Pozrikidis [1], on the same topic, deals with analogue eigenproblems. The main drawback of this approach consists in the fact that it depends strictly on the problem at hand.

A much more flexible and general strategy is considered in this paper. Mainly it transforms the original problem into one containing lower order derivatives in the differential equation as well as in the boundary conditions. Sometimes this means to rewrite a problem as a first order system of differential equations. The shooting method and the compound matrix method are examples of such method. Instead we suggest a method which takes into account the particularities of the problem (parity of order of differentiation, symmetries, etc.). It reduces a fourth or sixth order eigenvalue problem (generally, an even order problem) with Dirichlet and hinged boundary conditions to a system of second order equations supplied exclusively with Dirichlet boundary conditions. These boundary conditions are much more simple to handle than the hinged one. Further, this procedure simplifies considerably the construction of test and trial functions in tau and Galerkin methods as well as the formulation of differentiation matrices in the collocation (pseudospectral) method. Thus, it is applicable to any differential boundary value problem which contains only even order of differentiation supplied with Dirichlet and hinged boundary conditions. This is our main theoretical result.

All methods provide a fairly accurate approximation of the lower part of the spectrum without any scale resolution restriction. The Chebyshev collocation is the fastest and the most flexible with respect to the implementation process.

The Viola’s eigenvalue problem is considered as a benchmark one. This is a fourth order singularly perturbed differential eigenvalue problem. We use two strategies in order to solve it. First, we solve the problem as a second order one supplied only with Dirichlet boundary conditions (”$D^2$” strategy) and observe that the differential singularity deteriorates considerably the conditioning of the discretization matrices. In spite of this inconvenience, our Chebyshev collocation method provides numerical results which improves those reported by Straughan in the Appendix of his monograph [29]. Second, we handle the problem as a fourth order one using in our collocation method shape functions which fulfill all boundary conditions. Unfortunately even for moderate values of $N,$ i.e., $N=128,$ the algorithm becomes unstable showing that the ”$D^2$” strategy is plainly superior with respect to accuracy. One drawback of this strategy is the superior cost of the computation of algebraic eigenvalues. As we work with the qz algorithm in solving the algebraic eigenvalue problem, the cost is eight times larger than that of solving the fourth order problem. In spite of this aspect, even for the hydrodynamic stability problem, where the cost is 27 times larger, the computing time and the memory requirements remain reasonable.

The above ”$D^2$” strategy as well as some weighted residuals (Chebyshev Galerkin type) methods are eventually used in order to solve a hydrodynamic stability problem of Bénard type involving variable gravity field.

Variations in the gravity field occurring in, on and above the Earth’s surface and due to the atmosphere and fluid dynamics, have been intensively studied. On a laboratory scale the variations are ignored and the gravity field can be considered as constant. However, their presence in crystal growth, e.g. [31], in convection problems from porous media, e.g. [22] and other convective motions from biomechanics, chemistry or geophysics make them attractive from the applications point of view. Consequently, their effects on the stability bounds must be investigated.

In the present paper we are concerned with a problem governing the convection-conduction in a horizontal layer bounded by two impermeable walls of a fluid heated from below. Across the layer the gravity field will vary according to some specific laws. The heat conducting viscous fluid is contained in a layer situated between the planes $z=0$ and $z=h$. For $t>0$ the conduction and convective motion is governed by the conservation equations of momentum, mass and internal energy ([29], P. 93)

where, as usual, $\nu$ is the coefficient of kinematic viscosity, $\rho$ is the density, $\alpha$ the thermal expansion coefficient, $k$ the thermal diffusivity, $p$ the pressure, $T$ the temperature, $𝐯$ the velocity field and $𝐠(z):=g(1+\epsilon h(z))𝐤$ is the variable gravity field, with $g$ constant, $𝐤$ the unit vector in the $z$-direction, $\epsilon$ is a scale parameter and $\epsilon h(z)$ represents the gravity variation.

The boundary conditions appropriate to fixed surfaces and specified constant temperature, $T_U$ on $z=h$, $T_L$ on $z=0$ are (see [29])

We assume that $T_L>T_U$.

The stationary solution of the equation (1) has the form

with $\beta$ the pressure gradient and $T_0$ considered as a reference temperature. In order to investigate the linear stability of (3), perturbations are introduced, the obtained perturbations equations are non-dimensionalized and then the convective terms are removed [29]. Assuming normal modes perturbations, the following two-point problem governing the linear stability of this conduction stationary solution results [29]

Here $D:={\displaystyle \frac{d}{dz}}$, $Ra=R^2={\displaystyle \frac{\alpha g\beta h^4}{k\nu }}$ is the Rayleigh number and it represents the eigenparameter of the problem (4)-(5), $a$ is the wavenumber, $W$ and $\mathrm{\Theta }$ are the amplitudes of the vertical velocity and temperature perturbation. They form the eigenfunction of the eigenvalue problem (4)-(5).This case is considered in detail in the monograph of Straughan [29].

The boundary conditions on a free surface are that the perturbations of the components of the stress are zero which implies that a free surface behaves as a rigid one with tangential slip but without any tangential stress [10]. In our case, for normal modes perturbations, the boundary conditions reduce to

We will treat here only the free boundaries case so the corresponding eigenvalue problem is (4)-(6). In [15] it was proven that in this situation the principle of exchange of stabilities is valid.

Straughan provided in [29] numerical evaluation of the Rayleigh number in the case of rigid boundaries by using the energy method. In [7], for linearly decreasing gravity fields, we obtained two stability criteria using isoperimetric inequalities and showing that in this case the domain of stability enlarges. Stability bounds were found by us in the case of both free and rigid boundaries [6] using spectral methods based on trigonometric series. In [9] the eigenvalue problem (4)-(6) was studied using expansion series based on polynomials (Legendre and Chebyshev polynomials).

Taking into account the order of differentiation in (4)-(6) we introduce the new variable

and thus, the two-point boundary value problem (4)-(6) can be rewritten as the second order system

with the boundary conditions

This means a generalized eigenvalue problem.

It is worth noting at this moment that the eigenvalue problem (4)-(6) can be rewritten more compactly as a sixth order differential equation

plus Dirichlet and hinged boundary conditions, namely

The new variable $\mathrm{\Psi }$ introduced by the transformation (7) reduced this problem to a Dirichlet problem, i.e., (8)-(9). The mathematical formalism requires one more boundary condition to be added to the boundary conditions (11) in order to close the problem. We introduce the condition $DW=0$ at $z=0,1$ and observe that with usual variational arguments it is simply to show that $R^2$ is indeed positive.

It is fairly important to underline why we integrate numerically the problem (8)-(9) instead of the problem (10)-(11).

The main reason is the rapid worsening (increasing) of the condition number of the discretization matrices as both the order of differentiation and the spectral parameter $N$ (number of terms in the spectral series) get larger. For the Chebyshev differentiation matrices, the fact is illustrated in the Figure 1.

As we do not have theoretical estimations for the spectrum of our hydrodynamic stability problem we will use the following eigenvalue problem as a benchmark one. More than that we try to make more explicit the reason we prefer to transform high order eigenvalue problems in systems of differential equations containing only second order derivatives, i.e., ”$D^2$” strategy.

The Viola eigenvalue problem reads (see Straughan [29] P. 218 and the works quoted there)

where the parameter $\theta$ satisfies The problem becomes singular as $\theta \to 1^{-}.$

With a change of variable analogous to (7), i.e., $v:={(1-\theta x)}^3{D}^{2}u,$ the problem (12) reduces to an eigenvalue problem for the pair ${(uv)}^T$ containing only Dirichlet boundary conditions, namely

As we make use of Chebyshev polynomials, the linear transformation of independent variables $z:=2x-1$ shifts our problem on the close interval $[-1,1].$ To go further with our method we solve this eigenvalue problem by the Chebyshev collocation method. This method transforms the problem into the following generalized algebraic eigenvalue problem

where the matrices $𝐀$ and $𝐁$ are

and the submatrices $\mathrm{𝐃𝟐},$ $\mathrm{𝐃𝐱},$ $𝐈$ and $𝐙$ stand respectively for the second order differentiation matrix on the Chebyshev nodes of the second kind $𝐱:=\{x_1,\mathrm{\dots },x_{N-1}\},$ $\mathrm{𝐝𝐢𝐚𝐠}\left(1-\theta 𝐱\right),$ the unitary matrix of order $N-1$ and the matrix with all zero entries of order $N-1.$ It is important to underline that in the differentiation matrix $\mathrm{𝐃𝟐}$ the boundary conditions are enforced, thus the rows and columns corresponding to the first node $x_0$ and the last node $x_N$ are deleted. The vector $𝐰$ contains the values of the approximations of $u$ and $v$ in the nodes of $𝐱$, i.e., $𝐰:={(u_1,\mathrm{\dots },u_{N-1},v_1,\mathrm{\dots },v_{N-1})}^T$. All the differentiation matrices we have used (chebdif,cheb2bc) come from the work of Weideman and Reddy [32]. The problem (13) was solved using the MATLAB code eig. Our results confirm up to some point those of Straughan [29]. Apart from those we conjecture that the first eigenvalue ${\lambda }_1$ approaches $1$ as $\theta \to 1^{-}$! It is sustained by the numerical results reported in Figure 2. They were obtained when $N=512$ and corresponding to $\theta =1-\left(1.0e-06\right)$ we got ${\lambda }_1=1.054343.$ More than that, for $N=600$ and the same $\theta$ we obtained ${\lambda }_1=1.035803.$ We have to observe that all the computations involved by the collocation method were performed in full precision (64 bit arithmetic).

An inconvenience of our approach consists in the fact that corresponding to all zero rows in the matrix $𝐁$ the eigenvalues equal infinity (the rank of matrix $𝐁$ equals $N-1)$. The MATLAB code eig, which implements the qz algorithm and gives access to intermediate results in the computation of generalized eigenvalues, explains this phenomenon. In fact the generalized eigenvalues for those zero rows are zero.

As a further check on accuracy, we computed two types of sensitivities defined as

for the $i-th$ eigenvalue where ${𝐱}_i$ and ${𝐲}_i$ are the right and left eigenvectors of (13) respectively. They were introduced by Stewart and Sun in [28] and used in an analogous problem by Dongarra et al. in [5]. Corresponding to Chebyshev collocation method applied to problem (13) with $N=512$ we summarize the situation in the Table 1.

It is fairly clear from this table that both measures of sensitivities increase drastically with the decreasing of parameter $\theta .$ In other words, the more singular is the differential problem the more sensitive are the eigenvalues. Since the quantity ${\log }_{10}{S}^{\left(i\right)}$ is a measure of the number of digits one lose in accuracy this table provide useful information on the accuracy of the actual computations. A similar conclusion provides the Figure 3. As the curves ${\log }_{10}cond(A)$ versus $\theta ,$ for $N=256$ and $N=512$ respectively, are quite parallel and horizontal for they sharply increase for $\theta \to 1^{-}.$

An attempt was made to solve the problem (12) as a sixth order one. In order to fulfill the boundary conditions we use this time the more complicated weighted functions

and in order to obtain the corresponding differentiation matrices we make use of the code poldif from the same paper of Weideman and Reddy [32]. In this situation the right hand side matrix of the generalized eigenvalue problem is a diagonal one and consequently it is easy invertible. This way we are in the situation to solve an usual eigenvalue problem. For values of $N$ around $128,$ the MATLAB code qz as well as the MATLAB code invshift (see Quarteroni and Saleri [24] P. 175) which evaluates the eigenvalue of minimum modulus of a matrix by the inverse power method, provide positive values for ${\lambda }_1.$ However, they are far larger than those obtained above. For $N$ beyond $128$ the situation is even worse. The eigenvalue of minimum modulus becomes negative and also complex eigenvalues appear in the spectrum. This means numerical instability.

Actually we look for the smallest eigenvalue $Ra$ (the Rayleigh number) in the problem (8)-(9) defining the neutral manifold.

In order to accomplish that, we reduce the problem to generalized algebraic eigenvalue first by a weighted residual (Galerkin) type method. To this end we formulate the eigenvalue problem in an operatorial framework by introducing the vectorial unknown $𝐔$ as $𝐔:={\left(W\mathrm{\Psi }\mathrm{\Theta }\right)}^T$, the operator $𝐋:X\to X,$ and the functional space $X:=H_w^2(-1,1)\times H_w^2(-1,1)\times H_w^2(-1,1)$ where

and

$I$ and $O$ stand for the identity and null operator. The space $H_w^p(-1,1),p\ge 0$ is a subspace of the Sobolev space $H^p(-1,1).$ Endowed with the inner product ${(\cdot ,\cdot )}_{p,w},$ defined as usual by ${(u,v)}_{p,w}:={\displaystyle \sum _{k=0}^p}{\displaystyle {\int }_{-1}^1}{\displaystyle \frac{d^{k}u\left(x\right)}{d{x}^p}}{\displaystyle \frac{d^{k}v\left(x\right)}{d{x}^p}}w(x)𝑑x,$ with the weight $w\left(x\right):=1/\sqrt{1-x^2},$ it is a Hilbert space. When $p=0$ we suppress this index from the definition of the scalar product.

It is clear that ${𝐋}_{\mathrm{𝟏}}$ is invertible and ${𝐋}_{\mathrm{𝟐}}$ does not have this property.

As in the case of Viola’s problem we first translate the problem on the interval $[-1,1].$

Consequently, our problem reads: find $Ra$ and $𝐔\in X$ $\backslash \left\{0\right\}$ such that

Various numerical methods can be used in order to solve the problem (16). In the monographs [2], [3], [19] or [30] to quote but a few, details on the spectral methods are provided including the approximation results in the weighted Sobolev spaces $H_w^p(-1,1),p\ge 0.$ In our recent monograph [12] we tried to collect some classical and at the same time difficult eigenvalue problems and compare the results obtained by spectral methods with their exact solutions or at least other numerical results (see [23] for such a collection).

Let us consider the following two finite dimensional spaces $X_M$, $Y_M$, both subspaces of $X$, made up by functions which satisfy the Dirichlet boundary conditions $𝐔=0$  at  $x=\pm 1.$ In our case an approximation ${𝐔}_M\in X_M\setminus \left\{0\right\}$ of the solution $𝐔\in X$ has the form ${𝐔}_M={(W_N{\mathrm{\Psi }}_N{\mathrm{\Theta }}_N)}^T$, with the $N$ index related to the discretization process. As usual, $W_N$, ${\mathrm{\Psi }}_N$, ${\mathrm{\Theta }}_N$ are defined as truncated series of orthogonal functions (trial functions)

In (17) the coefficients ${\widehat{W}}_k$, ${\widehat{\mathrm{\Psi }}}_k$, ${\widehat{\mathrm{\Theta }}}_k$, $k=0,1,2,\mathrm{\dots },N$ are unknown and the sets ${\{e_k^W(x)\}}_{k=0,1,\mathrm{\dots },N}$, ${\{e_k^{\mathrm{\Psi }}(x)\}}_{k=0,1,\mathrm{\dots },N}$, ${\{e_k^{\mathrm{\Theta }}(x)\}}_{k=0,1,\mathrm{\dots },N}$ are the sets of trial functions. This imply $M=3(N+1)$. The spectral approximation ${𝐔}_M$ is obtained by imposing the vanishing of the projection of the residual on the finite dimensional space $Y_M$, subspace of $X$, so a general weighted residual formulation for the problem (16) reads: find $Ra$ and ${𝐔}_M$ $\in$ $X_M$ $\backslash \left\{0\right\}$ such that

Whenever we perform an integration by parts of $D^2$ in the above formulation, in order to lower the order of differentiation, the formulation is called week and the method becomes a Galerkin-Petrov one. When the space of trial functions $X_M$ and that of test functions $Y_M$ are identical we come up with the well known Galerkin method.

In [9] the analytical study of the eigenvalue problem (8)-(9) was performed by using complete sets of orthogonal functions directly on $[0,1]$, i.e. shifted Chebyshev and Legendre polynomials.

For the particular varying gravity field $h(x)=-(2x-1)$ we solve the problem (18) using several basis in the test and trial spaces.

i) Weighted residuals (Galerkin) methods (WRG). We considered two family of functions as basis of the finite dimensional spaces $X_M$, $Y_M$. They have the form ${\{F_j\}}_{j=1,2}={\{S_{jk}\times S_{jk}\times S_{jk}\}}_{j=1,2;k=0,1,\mathrm{\dots },N}$ with the Shen basis ${\{S_{1k}\}}_{k=0,1,\mathrm{\dots },N}$ (SB1) [27] defined by

and ${\{S_{2k}\}}_{k=0,1,\mathrm{\dots },N}$, which is another Shen basis (SB2)[27], having the form

Then, in each case, the unknowns $W$, $\mathrm{\Psi }$, $\mathrm{\Theta }$ are approximated respectively by

Let us introduce the vector of the unknown coefficients ${\widehat{W}}_k$, ${\widehat{\mathrm{\Psi }}}_k$, ${\widehat{\mathrm{\Theta }}}_k$ by

With this notation the equation (18) becomes

where the matrix $𝐋$ is defined by

which means a homogeneous linear algebraic system of $3(N+1)$ equations and $3(N+1)$ unknowns.

The corresponding secular equation has the form

All the scalar products from (22) are based on the properties of Chebyshev polynomials [19], i.e.,

with $c_i=\{\begin{array}{c}2,\text{ if }i=0,\hfill \\ 1,\text{ if }i\ge 1.\hfill \end{array}$

For the Galerkin method the secular equation becomes

ii) Weighted residuals (Petrov-Galerkin) methods (WRPG). We considered two variants. At first, in the method (WRPG1), the set of trial functions is of the form ${\{{\mathrm{\Psi }}_k\times {\mathrm{\Psi }}_k\times {\mathrm{\Psi }}_k\}}_{k=0,1,\mathrm{\dots },N}\in X_M$ with the Heinrich’s basis ${\{{\mathrm{\Psi }}_k\}}_{k=0,1,\mathrm{\dots },N}$ ([14]) defined as

and the set of test functions is $F_1$ defined as before. For the second variant (WRPG2) we choose the same set of trial functions as for (WRPG1) and as a set of test functions we choose $F_2$.

The secular equation is obtained the same way as in the WRG method.

iii) Collocation methods. The Chebyshev collocation method, (CC) for short, reduces the eigenvalue problem (8)-(9) to an algebraic problem of the form (13) where the matrices $𝐀$ and $𝐁$ are now

with the same significance for matrices $\mathrm{𝐃𝟐},$ $𝐈$ and $𝐙$ but $\mathrm{𝐃𝐱}$ is the diagonal matrix $\mathrm{𝐃𝐱}=\mathrm{𝐝𝐢𝐚𝐠}\left(\left(x+1\right)/2\right).$ The numerical results obtained for various values of $\epsilon$ and $a^2$ with $N=256$ are reported in the last column of the Table 3. They were obtained using the same MATLAB code eig.

With respect to the sensitivity defined by (15) we remark that this problem is less sensitive than the singularly perturbed Viola’s problem. For the first indicator we obtained $S^{\left(1\right)}=5.789210$ which increases to $S^{\left(2\right)}=15.846905$ and $S{1}^{\left(i\right)}$ slightly oscillates around $0.59$. It is important that they remains practically independent of $N.$

Numerical evaluations of the Rayleigh number for the particular case in which the gravity variation is $\epsilon h(x)=-\epsilon (2x-1)$ are given. For comparison, we present in Table 2, our previous evaluations of the Rayleigh number from [6], [9]. These evaluations are obtained by spectral methods based on shifted Chebyshev polynomials (SCP), shifted Legendre polynomials (SLP) and trigonometric series (TS). The numerical results obtained here by the weighted residual Galerkin methods (WRG)(presented in Table 3) (for the first Shen basis (WRG(SB1)) and for the second Shen basis (WRG(SB2)) are similar to the ones in Table 2 and underline the same physical conclusion that a linearly decreasing variable gravity field enlarges the domain of stability. The evaluations obtained by the Chebyshev collocation method (CC) are presented in the last column of the same table.

Then numerical evaluations of $Ra$, obtained for the two variants of the weighted residual Petrov-Galerkin method (WRPG1 and WRPG2), are displayed in the Table 4. These numerical results are similar. However, the use of the (WRPG) method simplifies the numerical evaluations and reduces the time of computation due to the fact that a larger number of scalar products in the secular equation vanish.

All in all, comparing the Tables 2, 3 and 4 it is fairly clear that the Chebyshev collocation method and the Galerkin method based on trigonometric polynomials provide the closest results.

The numerical scheme for the generation of the spectrum in the case of weighted residual methods based on trigonometric and polynomial expansions is a Maple code which implements the LU decomposition. In Table 5 the Rayleigh number for the classical case [4], i.e. $\epsilon =0$, $a^2=4.92$, obtained by (WRG) with the first Shen basis ${\{S_{1k}\}}_{k=0,1,\mathrm{\dots },N}$ are shown. The results are presented, for various spectral orders, in comparison with the necessary time for each evaluation. It is noticed that a larger spectral order is time consuming without having a distinctive improvement of the numerical results.

We can conclude that for a relatively small $N$, i.e., $N=12,$ fairly accurate numerical results are obtained using these methods. However, for larger $N$ the computational time highly increases, without leading to a significant improvement of the numerical values (less than ${10}^{-5}$).

A better situation is observed with respect to the collocation method. The spectral order $N$ used is much more large than that used in the weighted residual methods, i.e., $N$ is of order $2^m$ with $5\le m\le 10$. In spite of this, our numerical experiments showed that the collocation methods are much more rapid (the CPU time increases moderately with $N$) and more flexible with respect to the implementation process. However, both classes of methods provide reliable approximations for the lower part of the spectra. More than that, they do not require assumptions on the scale resolution, i.e., that $Ra/N^2$ is small, as it is the case with the methods considered by Melenk, Kirkner and Schwab in [20].

In Fig. 4 graphical representations of the neutral curves for different values of the scale parameter are given, pointing out its influence on the stability domain.

In this paper a numerical study of a linear hydrodynamic stability problem of a convective flow in the presence of a varying gravity field is carried out. The sixth order eigenvalue problem with Dirichlet and hinged boundary conditions is cast into a second order eigenvalue problem supplied only with Dirichlet boundary conditions. It is solved by some weighted residual methods (Galerkin type) methods as well as by a collocation method. The algorithm is fairly general and its Chebyshev collocation implementation very fast and flexible. Formally this implementation produces some spurious eigenvalues but this fact is explicable and expected. However, the lower part of the spectrum, which is highly important, is accurately approximated. More than that, no scale resolution restrictions are required.

The Chebyshev collocation version of our method was validated in solving the Viola eigenvalue problem. It is a singular eigenvalue problem with high sensitive eigenvalues. About this problem we conjecture that the first eigenvalue ${\lambda }_1\to 1$ as $\theta \to 1^{-}.$

Acknowledgements

The research of the second author was supported in part by the CEEX Grant 2-CEx06-11-96/19.09.2006. The authors are fairly grateful to the two anonymus referees for their careful reading of the paper and trenchant suggestions for improvement.