Nowadays the industry and ecology need more and more results from hydrodynamic stability theory for sophisticated fluid motions occurring in complicated circumstances. In certain of these situations, a direct application of numerical methods can lead one to false results due to the bifurcation problems of the stationary solutions set of the Navier-Stokes equations (or of some more general models) and to the dependence of the eigenvalue on physical parameters. That is why, an analytic along with a numerical study of the eigenvalue problems from hydrodynamic stability theory is highly requested. This is in fact one of the most difficult topic in hydrodynamic stability. The occurrence of false secular points in the linear stability of continua was first pointed out by Collatz in 1981 in his paper [4].

A considerable number of theoretical and numerical studies have been devoted to the interaction of electromagnetic fields with fluids. Rosenswieg [24] pointed out that there are three main categories on this subject, i.e. electrohydrodynamics (EHD), magnetohydrodynamics (MHD) and ferrohydrodynamics (FHD). Here we are interested in an eigenvalue problem in EHD which implies the presence of electric forces. The problem at hand consists into an eighth order differential equation, containing only even order derivatives, supplied with homogeneous boundary conditions for the even order derivatives up to sixth order, i.e. the so called hinged boundary conditions.

The main aim of this paper is to solve analytically as well as numerically this problem in order to get confidence in the later. In fact we are mainly interested in the low part of the spectrum of this problem.

With respect to the numerical methods we have to observe two important facts.

First, a straightforward application of the tau method based on Chebyshev polynomials leads to extremely ill conditioned matrices which are also fully populated (the asymptotic order of the entries varies between zero and $O\left(N^{2\times 8}\right)$ where $N$ is the spectral parameter). The eighth order Chebyshev differentiation matrices are also fairly bad conditioned. For $N=2^{9}$ the condition number of these matrices attains something of order $O\left(10^{30}\right)$ (see our paper [15]). Consequently, as we have shown, a direct application of the Chebyshev collocation method to this eighth order problem leads to huge instabilities.

Second, due to the fact that the boundary conditions imply a fairly challenging lacunar interpolation problem, the Galerkin and the collocation (pseudospectral) methods become directly inapplicable. It is worth noting in connection with this second point that Huang and Sloan, in their papers [20,21], consider a fairly general non lacunar interpolation problem with multiple nodes and then solve some fourth and sixth order eigenvalue problems. We also succeeded (see [13]) in solving a non-standard eigenvalue problem, i.e. an Orr-Sommerfeld equation supplied with boundary conditions containing the spectral parameter, concerning flows driven by surface tension gradients. As it is not the case with our problem, we follow our " $D^{2"}$ strategy from [15], i.e. we rewrite the problem as a second order differential system supplied only with Dirichlet boundary conditions (also called clamped boundary conditions). This way, all the three spectral methods, namely tau, Galerkin and collocation can be efficiently applicable. They provide fairly accurate approximations for the low part of the spectrum at a modest computational cost.

The rest of the paper is organized as follows. In Section 2 we introduce the problem and show that whenever the physical parameters satisfy an "ellipticity" condition the smallest eigenvalue is real and positive. In Section 3, a direct analytical method based on the characteristic equation is used in order to solve the eigenvalue problem. In Section 4, we carry out the above transformation of the problem into a second order system supplied with Dirichlet boundary conditions and solve this problem by analytical methods which use Fourier expansions. In Section 5 we solve this eigenvalue problem, first by standard Chebyshev tau method and try to explain its lack of accuracy. We observe that the matrices involved in this method are highly non-normal. In the last subsection we solve the eigenvalue problem using our " $D^{2}$ " strategy along with Galerkin methods. For various values of the physical parameters, the computed values of the critical Rayleigh number, are displayed in a table. For the considered numerical values of them, the numerical results confirm the analytical ones.

The linear stability of the stationary solution in an electrohydrodynamic convection model in a layer situated between the walls $z=\pm 0.5$, against normal mode perturbations, is governed by the following eigenvalue problem from [10]

Here $F$, which represents the amplitude of the temperature field perturbation, stands for the eigenfunction in (1). The physical parameter $a$ represents the wavenumber, $L$ is a parameter effectively measuring the potential difference between the planes and $R$ stands for the Rayleigh number.

Electrohydrodynamic systems have important industrial application in the construction of devices using the electroviscous effect or charge entrainment, for instance EHD clutch development and EHD high voltage generators.

The linear stability of their steady states typically leads to high order differential eigenvalue problems.
Thus, Roberts in [23] investigated two electrohydrodynamic convection models, based on the Gross’ experiments [17]. They were concerned with a layer of insulating oil confined between two horizontal conducting planes, heated from above and cooled from below. In the first model, the dielectric constant is allowed to vary with the temperature. The homogeneous insulating fluid is assumed to be situated in a layer of depth $d$ (the fluid occupies the region between the planes $z=\pm 0.5d$, which are maintained at uniform but different temperatures), with vertical, parallel applied gradients of temperature and electrostatic potential. The uniform electric field is applied in the $z$ direction. The second one, also investigated by Turnbull in [27,28], is characterized by a variation of the dielectric constant which is not important but the fluid is weakly conducting and its conductivity varies with temperature. The corresponding eigenvalue equation is simpler than (1) and has the form (see [26, p. 207])

with $M$ the dimensionless parameter measuring the variation of the electrical conductivity with temperature. The boundary conditions, written for the case of rigid boundaries at constant temperatures, read (cf. [26])

Roberts [23] solved the problem (2), (3) numerically in order to obtain the

He found that when the parameter $M$ is increasing in the range $0-1000$, the minimum value of $R,R_{\text{min }}$ is increasing from 1707.062 to 2065.034 . Unfortunately, a straightforward application of our " $D^{2}$ " strategy to problem (2), (3) is impossible due to the odd order of differentiation contained in both differential equation and boundary conditions.

Straughan in his monograph [26] also investigated these EHD convection problems, developing a fully nonlinear energy stability analysis for non-isothermal convection problems in a dielectric fluid. A bifurcation analysis of the problem (2), (3) was performed by us in [7] revealing the false secular points on certain parameters surfaces. A numerical study completed the investigation of the eigenvalue problem.

A first important remark is in order at this moment. By standard variational arguments (integration by parts and the imposing of boundary conditions) a weak formulation of the problem (1) can be obtained. It reads:
find $F\in H_{0}^{4}(I)$ and $R\in\mathbb{C}$ such that

with $H_{0}^{4}(I)$ representing the closure in the Sobolev space $H^{4}(I)$ of the set of infinitely differentiable compactly supported functions on I (cf. [2]).

Moreover, whenever the parameters $a$ and $L$ satisfy the "ellipticity" condition, i.e.

the problem reduces to a minimization one, i.e.

where

and

It is clear that the smallest (algebraic) eigenvalue is real and positive, i.e. there exists a first $R>0$ which satisfies (1).
In spite of the fact that in both formulations (4) and (6) the order of differentiation is halved, it remains too high in the perspective of numerical computations.

Most of the problems of hydrodynamic, electrohydrodynamic or hydromagnetic stability as well as bifurcation of solutions can be reduced to eigenvalue problems defined by systems of ordinary differential equations including a large set physical parameters. In the presence of more than one physical parameter and with a high order of differentiation in the system, it is, however, difficult to analyze how the most relevant eigenvalue of the system depends on parameters, since, in general, these systems are not selfadjoint and may have variable coefficients. That is why, a numerical study is usually performed in order to obtain the critical eigenvalues defining the neutral manifold.

The direct method based on the characteristic equation was first systematically applied to hydrodynamic stability problems by A . Georgescu and then extensively used by her group e.g. $[11,12,7]$. The method is one of the most simple methods to treat two-point problems for linear ordinary differential equations with constant coefficients. The characteristic equation associated with (1) reads (see [10])

or, by denoting $\mu:=\lambda^{2}-a^{2}$, it becomes

By means of the direct method, we write the general form of the solution of the two-point boundary value problem (1) in terms of the roots of the characteristic equation. The general solution of the equation depends on the multiplicity of the roots $\lambda_{i}$ of the characteristic equation associated with the eigenvalue problem. Whence the importance of discussing the multiplicity of $\lambda_{i}$ (see also [11]).

If the $\lambda_{i}$ ’s are distinct, then the general solution is a linear combination of hyperbolic sine and cosine functions. If $\lambda_{i}$ are multiple, of $m_{\lambda_{i}}$ order respectively, then the general solution is a product of a polynomial in $z$ by $\cosh\left(\lambda_{i}z\right)$ and $\sinh\left(\lambda_{i}z\right)$, of degree $m_{\lambda_{i}}-1$. Further, the introduction of the general solution into the boundary conditions leads to the secular equation.

The neutral manifolds, in particular the neutral curves, separate the domain of stability from the domain of instability. The bifurcation manifolds of the characteristic equation, or part of it, may be false secular manifolds.

In [10] an analytical investigation of the multiplicities of the Eq. (7) was performed. Although vanishing parameters ( $a=0$, $R=0,L=0$ ) are physically meaningless, for bifurcation reasons, these cases were also considered. When $a=0$, writing the general solution, it was proven in [10] that no secular points exists. For $a>0,L=0,R>0$ we point out that only specific points of the corresponding bifurcation manifolds of the characteristic equation are secular. In this case, the roots of the characteristic equation are

where

The secular equation has the form

The equation $\cosh\left(\lambda_{i}/2\right)=0$, the only one that can lead to secular points, has solutions for $i=3$, i.e $\cosh\left(\lambda_{3}/2\right)=0$ if and only if $\cos\left(\lambda_{3}/2\right)=0$. Thus $\lambda_{3}=-(2n-1)^{2}\pi^{2}$, and the only secular points are those situated on the manifolds

The critical value of the Rayleigh number belongs to $NC_{1}$ identical to the classical one from Chandrasekhar [3]. The secular manifold just found is confirmed in the following section by analytical methods based on Fourier series expansions. In order to obtain possible bifurcation sets of the characteristic equation, some other particular cases were investigated in [10]. The number of physical parameters, greater than two, as well as the high order of differentiation in the governing equation, hinder a systematic analytical investigation of the bifurcations of the involved manifolds. That is why, the problem required numerical methods, specific to bifurcation theory, in order to separate numerical solution of the characteristic equation and of the secular equation.

In order to apply these methods we rewrite the two-point boundary value problem (1) in a different form. First, we introduce the new vector function

with

This substitution turns the eighth order equation from (1) into the following system of second order ordinary differential equations

Consequently, the boundary conditions simplify considerably and become

It is important to underline that, in this way, the high order (hinged) boundary conditions from the problem (1) transform into simple Dirichlet boundary conditions. This was also the successful strategy in our previous work [15].

Most of the eigenvalue problems occurring in electrohydrodynamic stability theory consists of higher order ordinary differential equations with constant coefficients depending on physical parameters. Herein, the problem of linear stability in the electrohydrodynamic convection is an eighth order eigenvalue problem with constant coefficients depending on three parameters, namely $a,R$ and $L$. Consequently, the discussion of multiplicity of the roots of the characteristic equation when all the parameters are different from zero, becomes fairly laborious. In this way, the application of the analytic direct method becomes quite obscure and alternative methods must also be used. Some of these methods are based on Fourier series. In the
present context we consider two of the most known such methods: the Chandrasekhar-Galerkin method and the Budian-sky-DiPrima method.

With respect to the first one, the unknown functions are written in the form of expansions in Fourier series upon a complete set of trigonometric functions (in $L^{2}[I]$ ) which satisfy all boundary conditions. The characteristic equation is represented by the equality to zero of an infinite determinant.

The second one, is a direct expansion approach too, but the trigonometric functions are chosen such that they satisfy only part of the boundary conditions of the problem. In this case, the characteristic equation is the equality to zero of an infinite series.

Let us apply the Chandrasekhar-Galerkin method to (11), (12). To this aim, we expand the unknown functions $U_{1},U_{2},U_{3}$, $U_{4}$ upon sets of orthonormal functions that satisfy all boundary conditions. Further we write each of the unknown functions in problem (11), (12) as a sum of an odd function and an even function. In this way, taking into account that an even function is equal to an odd one only if they are both null, the problem splits into the following two-point boundary value problems

for odd part, and

for even part.
Assume that $U_{1},U_{2},U_{3}$ and $U_{4}$ are even functions of $z$ and for the sake of simplicity we will omit the indices in (15), (16). Anyway an analog analysis can be performed for the odd case.

However, better results are obtained for the even case (see [6] for such an example). Taking into account the boundary conditions (12), the unknown functions are expanded upon the orthonormal set in $L^{2}(I),\left\{E_{2n-1}\right\}_{n\geqslant 1}$, where

i.e.

The series expansions of the derivatives occurring in (11) are obtained by the backward integration technique (see the monograph [9]). We substitute these expressions in (11) and impose the condition that the residuals be orthogonal to $E_{2m-1}$, $m=1,2,\ldots$. In this way, for the unknown coefficients $U_{j_{n}}$, we get the system

where $A_{n}=(2n-1)^{2}\pi^{2}+a^{2}$. The following linear algebraic equation for the constant coefficients $U_{4_{n}}$ is obtained by eliminating $U_{j_{n}},j=1,2,3$ between the equations of (17)

This means the secular equation $R=\frac{A_{n}^{4}-La^{4}}{a^{2}A_{n}}$. So, for $n=1$, the surface

gives us the critical points for $L\neq 0$. For $L=0$, we get

which implies that the critical secular points ( $R,0,a$ ) belong to $NC_{1}$. This validates the conclusion obtained in the previous section with the direct method.

The Budianski-DiPrima method is based on the Fourier expansion of all unknown functions upon total sets of functions which do not satisfy all boundary conditions of the problem. The unfulfilled boundary conditions lead to some constraints for the Fourier coefficients. In our case, we consider the unknown functions expanded upon the total set $\left\{F_{2n-1}\right\}_{n\in\mathbb{N}}$ in $L^{2}(I)$, where

Then we substitute these expressions in (11), and impose the condition that the residuals be orthogonal onto $F_{2m-1}$, $m=1,2,\ldots$. We obtain the system

where $\alpha=DU_{4}(0.5),\beta=DU_{1}(0.5),\gamma=DU_{2}(0.5),\delta=DU_{3}(0.5)$ are arbitrary constants generated by the unfulfilled boundary conditions in the backwards integration technique. The determinant of the system has the form

Let us assume that $\Delta\neq 0$. Solving the system (21) and replacing the solution $\mathbf{U}$ into the constraints resulting from the boundary conditions, i.e.

the following linear algebraic system in the unknown $\alpha,\beta,\gamma,\delta$ is obtained

However, taking into account that the constant $\alpha,\beta,\gamma,\delta$ are not all null, i.e. the determinant of the algebraic system (22) vanishes, we find that $\Delta=0$. Consequently, the same conclusion arises: for $L\neq 0$ the critical secular points belong to $NC_{L}$ and for $L=0$ the critical point belongs to the manifold $NC_{1}$ given by (10).

The problem represents a specific class of eigenvalue problems for which the parity of the derivatives in the ordinary differential equations and in the boundary conditions allows a splitting of the problems in two parts - an even part and an odd one, both leading to the same neutral curve. This choice simplified the numerical analysis, since simplified trigonometric orthogonal functions were used. They confirmed the bifurcation analysis in the case of a vanishing physical parameter.

The linear stability problem of electrohydrodynamic convection between two parallel walls was accurately solved by two classes of analytical methods as well as by spectral methods. The eighth order eigenvalue problem was transformed into a second order system of differential equations supplied with (only!) Dirichlet boundary conditions. This new problem was first solved analytically by two methods based on Fourier series. The Fourier approximation was split into an even and an odd part. The analytical computations carried out on these two problems led to the same neutral curve. Then the Chebyshev Galerkin and Chebyshev tau along with our " $D^{2}$ " collocation strategy, introduced in our previous paper [15], were used in
order to find the smallest eigenvalue of the original problem. From the physical point of view the computed numerical values of the Rayleigh number lead to the following conclusion: when the parameter which effectively measures the potential difference between the walls is increasing the stability domain is decreasing.

For some values of the physical parameters the bifurcation results are accurately confirmed by our numerical computations.

The accuracy of the numerical results, as well as, the efficiency and the modest cost of the implementation of the " $D^{2}$ " strategy, shows once again its superiority over other spectral method to solve high order eigenvalue problems.

The first author would like to thank Prof. Adelina Georgescu for her encouragement, advice and support during the years of the author’s Ph.D. study. The work of the second author was partly supported by the Grant 2-CEx06-11-96.

[1] D. Bourne, Hydrodynamic stability, the Chebyshev tau method and spurious eigenvalues, Continuum Mech. Thermodyn. 15 (2003) 571-579.
[2] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics, SpringerVerlag, 2007.
[3] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, 1961.
[4] L. Collatz, Numerical methods for free boundary problems, in: Proceedings of Free Boundary Problems: Theory and Applications, Montecatini, Italy, 1981.
[5] J.J. Dongarra, B. Straughan, D.W. Walker, Chebyshev tau-QZ algorithm methods for calculating spectra of hydrodynamic stability problem, Appl. Numer. Math. 22 (1996) 399-434.
[6] F.I. Dragomirescu, Rayleigh number in a stability problem for a micropolar fluid, Turk. J. Math. 31 (2) (2007) 123-137.
[7] F.I. Dragomirescu, Bifurcation and numerical study in an EHD convection problem, An. St. Univ. "Ovidius" Constanta 16 (2) (2008) 47-57.
[8] D.R. Gardner, S.A. Trogdon, R.W. Douglas, A modified tau spectral method that eliminates spurious eigenvalues, J. Comput. Phys. 80 (1989) 137-167.
[9] A. Georgescu, Hydrodynamic Stability Theory, Kluwer, Dordrecht, 1985.
[10] A. Georgescu, D. Pasca, S. Gradinaru, M. Gavrilescu, Bifurcation manifolds in multiparametric linear stability of continua, ZAMM 73 (1993). 7/8 T767T768.
[11] A. Georgescu, L. Palese, On a method in linear stability problems. Application to a natural convection in a porous medium, Rapp. Int. Dept. Math. Univ. Bari 9 (1996).
[12] A. Georgescu, M. Gavrilescu, L. Palese, Neutral thermal hydrodynamic and hydromagnetic stability hypersurface for a micropolar fluid layer, Indian J. Pure Appl. Math. 296 (1998) 575-582.
[13] C.I. Gheorghiu, I.S. Pop, A Modified Chebyshev-Tau Method for a Hydrodynamic Stability Problem, in: Proceedings of ICAOR 1996, vol. II, 1996, pp. 119-126.
[14] C.I. Gheorghiu, Spectral Methods for Differential Problems, Casa Cartii de Stiinta Publishing House, Cluj-Napoca, 2007.
[15] C.I. Gheorghiu, F.I. Dragomirescu, Spectral methods in linear stability. Applications to thermal convection with variable gravity field, Appl. Numer. Math. 59 (2009) 1290-1302.
[16] D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods, SIAM, Philadelphia, PA., 1977.
[17] M.J. Gross, Mantles of the Earth and Terrestrial Planets, Wiley, 1967.
[18] P. Henrici, Bounds for iterates, inverses, spectral variation and fields of values of non-normal matrices, Numer. Math. 4 (1962) 24-40.
[19] A.A. Hill, B. Straughan, A legendre spectral element method for eigenvalues in hydromagnetic stability, J. Comput. Appl. Math. 193 (2003) 363-381.
[20] W. Huang, D.M. Sloan, The pseudospectral method for third-order differential equations, SIAM J. Numer. Anal. 29 (6) (1992) 1626-1647.
[21] W. Huang, D.M. Sloan, The pseudospectral method for solving differential eigenvalue problems, J. Comput. Phys. 111 (1994) 399-409.
[22] S.A. Orszag, Accurate solution of the Orr-Sommerfeld stability equation, J. Fluid Mech. 50 (1971) 689-703.
[23] P.H. Roberts, Electrohydrodynamic convection, Q.J. Mech. Appl. Math. 22 (1969) 211-220.
[24] R. Rosensweig, Ferrohydrodynamics, Cambridge Univ. Press, 1985.
[25] P.J. Schmid, D.S. Henningson, Stability and Transition in Shear Flows, Springer-Verlag, 2001.
[26] B. Straughan, The Energy Method, Stability, and Nonlinear Convection, second ed., Springer, Berlin, 2003.
[27] R. Tunbull, Electroconvective instability with a stabilizing temperature gradient. I. Theory, Phys. Fluids 11 (1968) 2588-2596.
[28] R. Turnbull, Electroconvective instability with a stabilizing temperature gradient. II. Experimental results, Phys. Fluids 11 (1968) 2597-2603.
[29] L.N. Trefethen, Computation of Pseudospectra, Acta Numer. (1999) 247-295.
[30] J.A.C. Weideman, S.C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. Math. Softw. 26 (2000) 465-519.