In a linear hydrodynamic stability analysis one first considers the linearization of the governing equations around a steady-state flow. The perturbation variables are described by a linear system of coupled differential equations supplied with a set of boundary conditions. The discretization of this system frequently leads to some non-Hermitian generalized eigenvalue problems with singular second order matrix in pencil. This fact was first observed in [1] in connection with stability of convective motion in a variable gravity field, in [2] where the authors analyze the stability of convective flows in porous media and in [3, 4, 5] where the authors consider various physical situations such as thermal convection in a box, a multi-component convection-diffusion and the Hadley flow. The same phenomenon is noticed in [6, 7, 8] for various incompressible flows. When the linear stability analysis of two-phase flows and free-surface flows (in MHD) are considered, the presence of eigenvalue at infinity is observed in [9] and [10].

The level of discretization needed to describe accurately the perturbed fields gives rise to large matrices, i.e., of order of thousands. This large dimension of problems rules out the calculation of the entire spectrum: usually, only the leading eigenvalues (those with the largest or the smallest real part) are computed. Up to our knowledge, the $QZ$ method was used in quasitotality of generalized eigenvalue problems encountered in linear hydrodynamic stability, see [11] for one of the first applications of $QZ$ in this context. Our papers [12] and [13] are not an exception to this situation. Unfortunately, the $QZ$ method has the main drawback that its complexity is of order $O(n^3)$ and due to large CPU times required it is hardly applicable to large problems. In order to eliminate non-physical eigenvalues in such non-Hermitian singular problems some ad-hoc methods were designed. For instance, in [8] the authors introduce a technique based essentially on the physics of the flow in order to eliminate these spurious eigenvalues. In [14] the authors use a transformation that maps the spurious eigenvalues to an arbitrary location in the complex plane. In [15] and [16] for spectral Galerkin and Chebyshev tau (Lanczos) respectively, the authors suggest practical tricks in order to remove the spurious modes. The authors of [6] obtain reliable results for the rightmost eigenvalue of Orr-Sommerfeld equation provided that a scale resolution assumption is fulfilled, i.e., the ration between the Reynolds number and the square of the cut-off parameter is ”small”.

Hence, there is need for a more specialized algorithm which computes only a few specific eigenvalues of the algebraic problem. In the case of hydrodynamic stability problems considered hereafter, one is typically interested in the eigenvalues closest to zero. Consequently, in this paper we propose to use (a real variant [17] of) the Jacobi-Davidson method [18] to solve the arising generalized eigenvalue problems. Advantages of the Jacobi-Davidson method are that it is applicable to large-scale eigenvalue problems and that it can deal with spurious eigenvalues at infinity [19]. Concerning the eigenvalues at infinity, in principle $JD$ can also be hampered by their presence. However, by clever selection strategies this can be avoided. Furthermore, in the paper [19] on purification it is described how shift-and-invert/Cayley transformation strategies can be employed to avoid convergence to eigenvalues at infinity.

However, we want to stress that, up to our knowledge, the $JD$ method has not been applied previously to these hydrodynamic stability problems and hence it is a novelty of our paper.

On the other hand, we are thoroughly aware of the existence of several methods designed to compute a specified subset of eigenspectrum, such as Krylov based methods (see the paper of Golub and van Loan [20] for the main research developments in the area). The reason we use $JD$ consists in the expertise we have got working with this method.

A partially similar study with ours is that of Valenttaro et al. [21]. The authors solve a two-dimensional eigenvalue problem by spectral collocation based on Legendre and Chebyshev polynomials. The generalized eigenvalue problem is solved by a Krylov subspace type method (incomplete Arnoldi-Chebyshev). For this method a detailed analysis of round-off error is provided. But, while our study refers to singular pencils Valenttaro et al. considers regular pencils.

The paper is organized as follows. In Section 2, we first introduce a linear hydrodynamic stability problem, comment on its weak and minimization formulations, and use the $\mathrm{`}\mathrm{`}{D}^2\mathrm{”}$ strategy from [12] to transform it into a second order system supplied with Dirichlet boundary conditions. By Chebyshev collocation the generalized algebraic eigenvalue problem is obtained. Then, a second and a third singular problems are discretized with the same strategy. In Section 3, the Jacobi-Davidson method is described. Numerical experiments are reported in Section 4 along with comments concerning the convergence and round-off errors involved. Conclusions are contained in Section 5.

In [12] the following even order eigenvalue problem was considered:

In these equations, as well as in the following two problems, $D$ stands for the derivative with respect to the variable $z$. $Ra:=R^2$ is the Rayleigh number which represents the eigenparameter of the problem (1)-(2), $a$ is the wavenumber, $W$ and $\mathrm{\Theta }$ are the amplitudes of the vertical velocity and temperature perturbation and $\epsilon h(z)$ signifies the gravity variation in the layer of fluid.

In order to reduce the order of differentiation in this problem we introduce the new variable

Consequently, the two point boundary value problem (1)-(2) can be rewritten as a second order system

supplied with the homogeneous Dirichlet boundary conditions

It is important to observe that the problem (1)-(2) can also be formulated as a sixth order differential equation

supplied with the hinged boundary conditions, namely

A weak formulation for the new problem (6)-(7) reads: find $w\in K\backslash \left\{0\right\}$ and $Ra\in ℂ$ such that

where the set $K$ is a closed subspace of $H^4(0,1)$ defined by

More than that, the critical Rayleigh number can be characterized by the following minimization problem

whenever the quantity $1+\epsilon h(z)$ remains positive. In fact, in all our experiments this condition is fulfilled.

The $\mathrm{`}\mathrm{`}{D}^2\mathrm{”}$ strategy implemented with Chebyshev collocation method to solve the problem (4)-(5) leads to the following generalized eigenvalue problem

where the matrices $𝐀$ and $𝐁$ are

The submatrices $\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}\},$ the unitary matrix of order $N-1$ and the matrix with all zero entries of order $N-1$, $N$ being as usual the spectral (cut-off or resolution) parameter. They have the same significance in the forthcoming two problems. The matrix $\mathrm{𝐃𝐱},$ signifies the diagonal matrix $\mathrm{𝐝𝐢𝐚𝐠}\left(h(𝐱)\right)$. 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. All the differentiation matrices we have used come from [22]. The matrix $𝐁$ is singular of rank $N-1.$

As a second test problem we consider another representative singular hydrodynamic stability problem. In this case we solve an eight order differential system which models multi-component convection-diffusion in a porous medium. Defining $a^2$ to be the wavenumber and $\lambda$ the eigenparameter, the non-dimensional linear perturbation equations read ([4]):

supplied with the boundary conditions:

In the system (12), $W$, $S$, ${\mathrm{\Psi }}^1$, ${\mathrm{\Psi }}^2$ are the perturbations of velocity, temperature and two solutes, $R$ and $R_{\beta }$ are thermal and solute Rayleigh numbers, respectively, and $P_{\beta }$ are salt Prandtl numbers, $\beta =1,2$ . Eventually $\zeta$ signifies a quantity connected with the temperature of the upper boundary. The $\mathrm{`}\mathrm{`}{D}^2\mathrm{”}$ strategy casts this problem into a generalized eigenvalue one with the matrices $𝐀$ and $𝐁$ defined as follows:

As a third test problem we consider the stability of the so called Hadley flow. It refers to convection in a layer of porous medium where the basic temperature field varies in the vertical (i.e. $z$ direction) as well as along one of the horizontal directions which is defined as $x$ direction. The nondimensional perturbation equations are ([4])

supplied with boundary conditions

In (15), $W\left(z\right)$ and $S\left(z\right)$ are respectively the third component of velocity and temperature field perturbations. The steady-state solutions have the form

where $R_H$ and $R_V$ are the horizontal and vertical Rayleigh numbers, respectively, $a^2:=k^2+m^2$ with $k$ and $m$ being the $x$ and $y$ wavenumbers and $\sigma$ is the eigenparameter. The same strategy casts this problem into a generalized eigenvalue one with the matrices $𝐀$ and $𝐁$ defined as follows:

where $𝐗:=diag(x^2/4)$, and $\mathrm{𝐃𝟏}$ is the first order differentiation matrix. We finish this section with the remark that the problem (15)-(16) can be reduced to a fourth order one in $W$ containing the leading term ${\left(D^2-a^2\right)}^2.$

Our strategy is strongly motivated by the conditioning of the matrices $𝐀$, $\left(4\mathrm{𝐃𝟐}-{𝐚}^2𝐈\right)$ and its second, third and fourth powers. In Fig. 1 the log of the condition numbers of these matrices is depicted. In logarithmic terms, this picture shows that conditioning of the fourth order differentiation matrix doubles, and that conditioning of the sixth order differentiation matrix triples, the conditioning of second order matrix. The conditioning of the eighth order differentiation matrix is (roughly) only $3.43$ times larger than that of second order differentiation matrix. It means powers two, three or $3.43$, of the conditioning of second order differentiation matrix, for the conditioning of higher order differentiation matrices, respectively. It is the main reason for which we do not applied directly collocation scheme to the problem (6)-(7) or to its weak formulation (8) or even to its minimization formulation (9). A completely analogous remark holds for the problems (12)-(13) and (15)-(16).

As mentioned in the introduction, the major drawback of full space methods such as the $QR$ method (for ordinary eigenvalue problems) and the $QZ$ method (for generalized eigenvalue problems) [23] is their computational complexity of $O(n^3)$, where $n$ is the dimension of the problem. Consequently, for $n\gg 5000$, direct computation of the complete spectrum is no longer feasible using the $QR$ or $QZ$ method (see the uppermost curve in Fig. 2).

On the other hand, for the application described in this paper (and for many other applications as well), one is not interested in the complete spectrum, but in a few specific eigenvalues: the eigenvalues closest to zero. In fact the eigenvalue closest to zero is the most important. Before describing the proposed algorithm, first the requirements for an eigensolution method for the application of this paper are summarized:

1. 1.

   Given the number $k$ of wanted eigenvalues and a target (e.g., closest to zero), the methods must produce $k$ eigenvalues closest to the target.

2. 2.

   The method must be scalable with respect to the dimension of the eigenvalue problem.

3. 3.

   The method must produce real approximate values for real eigenvalues and complex (conjugated pairs of) approximate values for complex eigenvalues.

4. 4.

   The method must not be hampered by eigenvalues at infinity of the eigenvalue problems (10)-(11), (14) and (17).

In the following we briefly describe the Jacobi-Davidson method [18], an iterative method for the computation of a few specific eigenvalues. We also show why this method satisfies the above requirements. For more details, see [24, 18, 17].

In this section we will report some representative numerical results obtained in solving the generalized algebraic eigenvalue problems defined by (10), (14) and (17). All numerical computations reported in [12] and [13] were carried out using the MATLAB code eig, which is an implementation of the QZ method. Alternatively, we use in this paper the $JDQZ$ implementation from

the $rJDQZ$ implementation described in [17], and compare the results with those obtained using the MATLAB code eig for the $QZ$ method (see also Table 3 in [12]), and the MATLAB code $eigs$ for the Arnoldi method. All computations were carried out using MATLAB 2010a on a HP xw8400 Workstation with clock speed of 3.2Ghz/s.

The computed values of the first three eigenvalues, i.e. $R_1,R_2,R_3$, ${\lambda }_1,{\lambda }_2,{\lambda }_3$ corresponding to problems (10) and (14) respectively are displayed in Table 1. The eigenvalues ${\sigma }_1,{\sigma }_2,{\sigma }_3$ corresponding to problem (17) are reported in Table 2. The CPU times reported in these tables, as well as in Fig. 2, measure the time required for the computation of the first six eigenvalues except QZ method which computes the whole spectrum.

They are obtained for the following sets of fixed physical parameters: $a^2=21.344$, $\zeta =0.14286$, $R=228.009$, $R_1=-291.066$, $R_2=261.066$, $P_1=4.5454$, $P_2=4.7619$, in case of problem (14), $\epsilon =0.03$ and $a=4.92$ in case of problem (10), and $k=0$, $m=10$, $R_H=114.2$, $R_V=100$, in case of problem (17).

It is fairly clear from both Tables that most numerical approximations of the first three leading eigenvalue computed by both $JDQZ$, $QZ$ and Arnoldi are almost indistinguishable. Our numerical experiments concerning problem (14), the largest one, extended successfully for all four methods for larger values of spectral parameter $N$, i.e., for $N$ in the range $512$, $2500$. This involved matrices of order $n=(N-2)\ast 4$, i.e., up to $n=10000$. As the values of the first two eigenvalues, up to sixth digit, remain unchanged for such range of $N$, it implies convergence of the numerical method.

The CPU times required by the all four methods are depicted in Fig. 2. It is important to observe that $rJDQZ$ computed only the first two eigenvalues for spectral parameter $N$ larger than $1000$.

In spite of its considerable robustness for large matrices, $JDQZ$ failed to compute the first eigenvalue corresponding to problem (17) (see second column in Table 2). The reason consists in the very bad conditioning of problem defined by (17).

The real variant of $JDQZ$, $rJDQZ$, does produce real eigenvalues because it uses real search spaces, as explained in Section 3.2.

Both Arnoldi and the $JD$ based methods did not show any difficulties with spurious eigenvalues.

We also have to emphasize that, as far as we know, we have used the largest spectral order $N$ in such eigenvalue computations. Comparatively, in [2] and [4] values of $N$ attain $50$, in [9] the authors work with $N$ of order $100$, in [8] the authors consider that $N=602$ is fine enough to predict correctly the leading eigenvalues of the problem and in [6] this order attains $1000$. In older papers the authors simply work with $N$ much less than $100$.

Some more comments on the convergence and round-off errors in $rJDQZ$ and Arnoldi method when singular generalized eigenvalue problems are solved are now in order. We define first, the error of a method as the absolute value of the difference between the computed eigenvalue and the converged one, i.e., obtained with the largest resolution. Than, we observe that the pseudospectra of our pencils $(A,B)$ can be unbounded, i.e., they extend to the whole complex plane for sufficiently large perturbations (see van Doersselaer [30]). Thus, to get some information we cut up from the whole spectrum the region surrounding the first two eigenvalues and observe that our computation is backward stable in the sense that the line of level which equals the error encloses the computed eigenvalue. Using a resolution of $256$ which implies a matrix of order $1.e+03$ the pseudospectrum is depicted in Fig. 3.

As the diameter of the enclosed area provides a hint of the largest possible relative error on the respective eigenvalue, we can infer a more sensitive second eigenvalue than the first one.

The errors of the first two computed eigenvalues, by $rJDQZ$ and respectively Arnoldi, vs. $1/N$, as is usual in spectral methods, are depicted in Fig. 4.

For resolutions larger than $1500$ this picture shows that the computation of the second eigenvalue becomes unstable in both methods. It is in accordance with the conclusion provided by the pseudospectrum. In order to observe the order of convergence of both methods we depicted Fig. 5.

The $rJDQZ$ method for resolutions inferior to $1700$ have order of convergence that could attains unity, i.e., $Err$ is of order . These resolutions are satisfactory in hydrodynamic stability problems. Arnoldi method behaves better. It has a comparable order of convergence but looks more accurate. It is fairly clear that $rJDQZ$ method is seriously affected by rounding off errors in Chebyshev differentiation matrices.

With respect to Arnoldi method we have to observe that this method confirms its performances established in [21] also for generalized singular eigenvalue problems.

In this paper we have shown that the Jacobi-Davidson and Arnoldi methods can be successfully used to compute the smallest eigenvalues of large-scale matrix pencils arising in hydrodynamic stability problems.

Both methods are not hampered by eigenvalues at infinity. Since these methods are applicable to large-scale problems, they are preferred over the $QZ$ method, which is only applicable to small eigenvalue problems. The variant of the Jacobi-Davidson method that uses only real search spaces, $rJDQZ$, takes advantage of using such search spaces resulting in real approximations of real eigenvalues. Consequently, it is the fastest method for solving such large and singular generalized eigenvalue problems. Moreover, this method is completely independent of the physical formulation of the problem at hand.

Eventually, we want to underline that the Chebyshev collocation (discretization) method proved again to be very implementable and feasible, for fairly large values of spectral parameter, provided that accurate differentiation matrices are used.