Computing eigenvalues occurring in continuation methods with 
the Jacobi­Davidson QZ method 
Jos L.M. van Dorsselaer \Lambda 
February 1997 
Abstract. Continuation methods are a well­known technique for computing several sta­ 
tionary solutions of problems involving one or more physical parameters. In order to 
determine whether a stationary solution is stable, and to detect the bifurcation points 
of the problem, one has to compute the rightmost eigenvalues of a related, generalized 
eigenvalue problem. The recently developed Jacobi­Davidson QZ method can be very 
effective for computing several eigenvalues of a given generalized eigenvalue problem. In 
this paper we will explain how the Jacobi­Davidson QZ method can be used to compute 
the eigenvalues needed in the application of continuation methods. As an illustration, the 
two­dimensional Rayleigh­B'enard problem has been studied, with the Rayleigh number 
as a physical parameter. We investigated the stability of stationary solutions, and several 
bifurcation points have been detected. The Jacobi­Davidson QZ method turns out to be 
very efficient for this problem. 
1 Introduction 
In physical applications one is often interested in stationary solutions of partial differential 
equations, and how their behaviour depends on (some) physical parameter(s) in the model. 
For instance, one would like to know whether stationary solutions (if they exist) are stable. 
At some critical values of the physical parameter, the so­called bifurcation points, a stable 
solution may become unstable and vice versa; further the number of stationary solutions can 
change at bifurcation points. Clearly, the computation of stationary solutions and bifurcation 
points is important for analyzing the physical problem under consideration. 
In practice, continuation methods [10] are often used to compute stationary solutions 
for different values of the physical parameter(s). With this approach one may also find 
unstable stationary solutions; these unstable stationary solutions have no physical relevance, 
but they might change into stable ones, when the physical parameter passes a bifurcation 
point (an example of this will be given in Section 4.2). For the investigation of stability, and 
the determination of bifurcation points, one has to compute some eigenvalues of a certain 
generalized eigenvalue problem 
Aq = –Bq; 
(1.1) 
a stationary solution is stable if all eigenvalues – of (1.1) have negative real parts; if at least 
one of the eigenvalues has a positive real part, the stationary solution is unstable. When one 
\Lambda Mathematical Institute, Utrecht University, P.O. Box 80.010, NL­3508 TA Utrecht, The Netherlands. This 
research has been supported by the Netherlands Organization for Scientific Research (N.W.O.). 
1 

2 Jos L.M. van Dorsselaer 
of the eigenvalues of (1.1) equals zero (i.e. the matrix A is singular), the physical parameter 
is a bifurcation point, and the number of stationary solutions may change [10]. 
In one continuation step one has to solve a system of non­linear algebraic equations (in 
order to obtain a stationary solution) and to compute some eigenvalues (the rightmost ones 
and those closest to zero) of (1.1). The determination of these eigenvalues is the most expen­ 
sive part of the computation, both in CPU­time and memory requirements. In this paper we 
will focus on the computation of these eigenvalues. 
For small problems (1.1), one can compute all eigenvalues with the QZ method (see e.g. 
[6]). However, this is not feasible for larger problems (1.1), e.g. those obtained from partial 
differential equations in 2­ or 3D (the size of the matrices A and B is equal to the number of 
unknowns obtained after discretizing the partial differential equations). For these problems 
one should use other methods; a well­known technique is the power method (see e.g. [6]), and 
block versions of this method (like SIT, see [13, 3]) have been used to obtain more eigenvalues 
of a problem equivalent to (1.1) (see, e.g. [3]). However, these methods can be very slow in 
practice. The fact that (1.1) is a generalized eigenvalue problem (i.e. B is not the identity 
matrix I ; in fact, B is often singular) may cause some extra complications when applying 
these methods. In the last decade some promising eigenvalue methods have been developed 
(see [7] for an overview of and references to such methods). One of these methods is the 
so­called Jacobi­Davidson QZ method developed by Fokkema, Sleijpen and Van der Vorst [5]. 
The main purpose of this paper is to show that this Jacobi­Davidson QZ method can be very 
efficient for computing the rightmost eigenvalues of (1.1), and those closest to zero. 
As an application of the Jacobi­Davidson QZ method, we consider the two­dimensional 
Rayleigh­B'enard problem, with the Rayleigh number as a physical parameter. The bifurcation 
behaviour of this problem has been studied extensively in [3]; in this paper the SIT method has 
been used to compute some of the rightmost eigenvalues to (1.1). Our experiments show that 
the Jacobi­Davidson QZ method is more efficient than the SIT method for this example. The 
successful application of the Jacobi­Davidson QZ method to the 2D Rayleigh­B'enard problem 
suggests that this method might be suitable for investigating the bifurcation behaviour of 3D 
flow problems in the near future. 
This paper is organized as follows. In Section 2.1 we describe how stationary solutions 
can be found by using continuation methods. The relation between stability of stationary 
solutions and the eigenvalue problem (1.1), as well as the concept of bifurcation points, will 
be discussed in Section 2.2. Section 3 deals with the Jacobi­Davidson QZ method and its 
application to continuation methods. The Jacobi­Davidson method, an essential ingredient 
of the Jacobi­Davidson QZ method, will be described in Section 3.1, and the Jacobi­Davidson 
QZ method is presented in Section 3.2. The application of the Jacobi­Davidson QZ method 
in combination with continuation methods will be discussed in Section 3.3. In Section 4 we 
present an illustration of the Jacobi­Davidson QZ method. The two­dimensional Rayleigh­ 
B'enard problem, and its discretization will be described in Section 4.1. Our numerical ex­ 
periments are given in Section 4.2, and the main conclusions of the paper are summarized in 
Section 4.3. 

Computing eigenvalues in continuation methods with the JDQZ method 3 
2 Continuation methods and the related eigenvalue problem 
2.1 Computing stationary solutions with continuation methods 
Consider the systems of differential equations 
BY 0 (t) = F (Y (t); ¯) for t – 0; 
(2.1) 
with B 2 R n;n , ¯ 2 R, Y (t) 2 R n for t – 0, and F : R n+1 ! R n is a smooth function. For 
y 2 R n the n \Theta n matrix 
F 0 (y; ¯) 
stands for the Jacobian matrix of the function F with respect to y, and the vector F ¯ (y; ¯) 2 
R n contains the partial derivatives (@F j =@¯)(y; ¯). 
In (2.1), ¯ stands for a `physical' parameter. Although the solutions to (2.1) depend on ¯, 
we will not express this dependence in order to keep the notation transparent. Also (systems 
of) time­dependent partial differential equations lead to (2.1), after discretizing the spatial 
derivative(s) with e.g. a finite difference or a finite element method. 
In many applications one is interested in finding stationary solutions to (2.1), for different 
values of ¯. Let y = y(¯) be such a stationary solution, i.e. 
F (y(¯); ¯) = 0: 
(2.2) 
Continuation methods [10] are often used for computing stationary solutions. A typical 
example of a continuation method is given in Algorithm 1. Suppose a stationary solution 
y 0 = y(¯ 0 ) and \Delta¯ 6= 0 are given. 
¯ 1 = ¯ 0 + \Delta¯ 
~ 
y 1 = y 0 + \Delta¯ \Delta y 0 (¯ 0 ) 
solve F (y 1 ; ¯ 1 ) = 0 with Newton's method; use ~ 
y 1 as starting vector 
for j = 2; 3; : : : ; M do 
¯ j = ¯ j 
(y j ; ¯ j ) = 0 with Newton's method; use y j as starting vector 
end for 
Algorithm 1. A continuation method. 
Algorithm 1 generates M stationary solutions y j = y(¯ j ). The vector y 0 (¯ 0 ) is the solution 
to F 0 (y 0 ; ¯ 0 )y 0 (¯ 0 ) = 
are O(\Delta¯ 2 ) approximations to y j . The determination of y j requires less computational costs 
(no linear system has to be solved), and therefore y j is used as a starting vector in Algorithm 
1. 
In practice one often uses (pseudo­)arclength continuation (cf., e.g., [8, 10]). An advantage 
of this technique is that it can handle turning points [10], while Algorithm 1 cannot. A major 

4 Jos L.M. van Dorsselaer 
drawback of arclength continuation is that two linear systems of order n have to be solved in 
each Newton step (see [3]) --- instead of one in Algorithm 1. 
More discussion about continuation methods can be found in [10, Chapter 4]. 
The linear systems which occur in Algorithm 1 have the matrix F 0 (y; ¯) as a coefficient 
matrix. When (2.1) is obtained from a set of partial differential equations, this matrix is 
usually large and sparse. For some of these problems it is possible to solve the linear systems 
with a direct solver, using important properties of the Jacobian matrix (like sparsity, small 
bandwidth etc.). But, for large general problems one has to solve these linear systems iter­ 
atively, using e.g. a Krylov subspace method (such as GMRES [9] or BiCGstab(`) [12]) in 
combination with a suitable preconditioner. Finally we note that there exist methods for solv­ 
ing nonlinear equations which combine the ideas of Newton iteration with Krylov subspace 
techniques (see, e.g., [2, 4]). These methods may be more efficient for solving F (y j ; ¯ j ) = 0 
in Algorithm 1 than Newton's method. 
2.2 Stability and bifurcation points 
We call a stationary solution y to (2.1) stable if lim 
t!1 
Y (t) = y for all Y (0) close to y. In order 
to investigate whether y is stable, equation (2.1) is linearized around y (note that F (y; ¯) = 0): 
BY 0 (t) = F 0 (y; ¯) \Delta (Y (t) 
alue of (2.3) has a 
positive real part, the equilibrium y is unstable [10, p. 20]. Hence investigating the stability 
of equilibria amounts to computing the eigenvalues (with largest real part) of (2.3). 
The number ¯ is called a bifurcation point if the matrix F 0 (y; ¯), with y = y(¯) a stationary 
solution, is singular (cf., e.g., [10]). This matrix is singular if and only if 0 is an eigenvalue of 
(2.3). The number of stationary solutions to (2.1) may change at a bifurcation point. 
In order to investigate the stability of stationary solutions, and the computation of bifur­ 
cation points, we have to determine in each continuation step the rightmost eigenvalues to 
(2.3), and those closest to 0. The Jacobi­Davidson QZ method, to be described in the next 
section, is well suited for this. 
3 The Jacobi­Davidson QZ (JDQZ) method and its applica­ 
tion in continuation methods 
The JDQZ method has been developed recently by Fokkema, Sleijpen and Van der Vorst [5]. 
In this section, the JDQZ method will be described briefly; for more details and discussion, 
see [5]. 
With the JDQZ method one can compute several eigenvalues (and eigenvectors) of the 
generalized eigenvalue problem 
fiAq = ffBq; 
(3.1) 

Computing eigenvalues in continuation methods with the JDQZ method 5 
here A, B are n \Theta n matrices with complex entries, and ff, fi 2 C. The pair hff; fii is called 
an eigenvalue, and q is the corresponding eigenvector. We write the eigenvalue problem in 
the form (3.1), instead of (1.1) (note that – = ff=fi), because fi = 0 is possible; when B is 
singular, fi = 0 for at least one eigenvalue. In the Rayleigh­B'enard problem (see Section 4), 
and in many other applications from fluid dynamics, the matrix B is indeed singular. 
In Section 3.1 we describe the Jacobi­Davidson method (see [11]), a method to compute 
one eigenvalue of (3.1); this method is an essential ingredient of the JDQZ method. The topic 
of Section 3.2 is the JDQZ method. The application of the JDQZ method in combination 
with continuation methods will be discussed in Section 3.3. 
3.1 The Jacobi­Davidson (JD) method 
With the JD method [11] one tries to compute an approximation h~ff; ~ 
fii ß hff; fii, close to a 
specified target ø (i.e. ~ ff= ~ 
fi should be close to ø ), and an approximate eigenvector ~ 
q ß q. In 
each step a search subspace spanfV g containing the vector ~ 
q and a test subspace spanfWg 
are constructed; V and W are complex n \Theta j matrices with j ø n and V \Lambda V = W \Lambda W = I . 
The vector ~ 
q and h~ff; ~ 
fii are obtained from the projected eigenvalue problem 
~ 
fiW \Lambda AV u = ~ 
ffW \Lambda BV u: 
(3.2) 
The eigenpair (h~ff; ~ 
fii; u) closest to the target ø is selected (i.e. ~ 
ff= ~ 
fi should be as close as 
possible to ø ); the vector ~ 
q = V u is an approximation to the eigenvector q. Throughout this 
paper, the approximate eigenvalue h~ff; ~ 
fii is scaled such that j ~ 
ffj 2 + j ~ 
fij 2 = 1. 
In order to improve the approximations, the spaces spanfV g and spanfWg will be ex­ 
panded in the next step; compute the residual r = ~ 
fiA~q 
i 
1 + jø j 2 
j j 
the QZ 
method (see, e.g., [6]) to compute all eigenvalues and eigenvectors to (3.2). 
The space spanfV g is expanded with the vector v which is orthogonal to ~ 
q and satisfies 
(I 
ect to the columns of the n \Theta j matrix V , and then added to the matrix V . In a 
similar way, the matrix W is enlarged with the vector w. This procedure is repeated until 
krk 2 is small enough. 
When the space spanfV g becomes too large, it is possible to restart the JD method; one 
might, e.g., replace spanfV g by the vector ~ 
q and spanfWg by ~ 
z, and repeat the procedure 
described above. A more efficient restarting procedure will be discussed in Section 3.2. 
The JD method converges quadratically, if (3.3) is solved exactly. Solving (3.3) is not 
trivial, because of the different projections involved; see Section 3.2 for details. 
Other choices for Ÿ 0 , Ÿ 1 , and the projections in front of and after ~ 
fiA 
6 Jos L.M. van Dorsselaer 
3.2 The JDQZ method 
The purpose of the JDQZ method is to determine a partial generalized Schur form 
AQ k = Z k S k ; BQ k = Z k T k ; 
(3.4) 
here Q k and Z k are n \Theta k matrices with Q \Lambda 
k Q k = Z \Lambda 
k Z k = I , and S k and T k are k \Theta k upper 
triangular matrices. From (3.4) one easily obtains k eigenvalues of (3.1) (and, optionally, the 
corresponding eigenvectors): note that fiS k x = ffT k x (x 2 C k ) implies fiA Q k x = ffB Q k x. 
The columns of the matrix Q k are called generalized Schur vectors. 
The first column of Q k (`the first Schur vector') is an eigenvector of (3.1), and we use the 
JD method to compute this Schur vector. Suppose a partial Schur form 
AQ k Z k 
a Schur vector) of (3.5). 
To apply the JD method we construct n \Theta j matrices V and W with V \Lambda V = W \Lambda W = I , 
and the extra condition V \Lambda Q k 
same problem as (3.2) (in exact arithmetic), because V \Lambda Q k 
B~q, and scale 
the vectors ~ q and ~ z such that k~qk 2 = k~zk 2 = 1. The search space spanfV g will be expanded 
with the vector v satisfying 
Q \Lambda 
k 
q~q \Lambda )v = 
acceptable 
approximation for a new Schur vector q has been detected. This vector will be added to 
the matrix Q k 
(3.7) can be solved, we will 
discuss how to determine a new matrix V after the detection of a Schur vector. 
When a Schur vector has been found, we have to restart the JDQZ proces with a different 
matrix V , because the relation V \Lambda Q k = 0 is violated. One might replace V by a vector 
v satisfying v \Lambda Q k = 0, but with this choice one might discard information in V regarding 
the new Schur vectors. Also when the space V becomes too large one would like to restart 
without losing valuable information. Both kind of restarts can be done efficiently, using the 
generalized Schur form related to the projected system (3.2) (which is equivalent to (3.6)): 
W \Lambda AV Q = ZS; W \Lambda BV Q = ZT ; 
(3.8) 

Computing eigenvalues in continuation methods with the JDQZ method 7 
here Q and Z are j \Theta j matrices with Q \Lambda Q = Z \Lambda Z = I and S and T are j \Theta j upper triangular 
matrices with diagonal elements s i and t i , respectively. The generalized Schur form (3.8) will 
be ordered such that 
js 1 =t 1 
y of V Q it follows that ~ 
q is perpendicular to the other columns of V Q. Therefore 
we restart with V := V Q 2;j (Q 2;j is the matrix consisting of the 2 nd ; 3 rd ; : : : ; jth column of 
V Q); this new matrix V satisfies V \Lambda V = I , V \Lambda Q k = 0, and contains as much information 
of the old V as possible. Further we set W := WZ 2;j . In a similar fashion we can restart 
when the matrix V becomes too large. From (3.9) we may argue that the first columns of 
V Q contain more important information about the Schur vectors to be detected than the 
last columns. One might replace V by V Q 1;j min , and W by WZ 1;j min , where j min ! j, and 
continue the process. 
It is not clear how to solve (3.7) in practice, because of the projections involved. Let K 
be a nonsingular n \Theta n matrix and denote e 
Q k = [Q k 
\Lambda 
k 
e 
Y k . It is shown in [5] that 
(3.7) is equivalent to 
e 
Q \Lambda 
k v = 0 and (I 
solve (3.10). The performance of Krylov subspace methods can often be improved by using 
some kind of preconditioning. The matrix K in (3.10) may be interpreted as a preconditioner, 
and in [5] it is proposed to take K ß A 
factorization 
of the matrix A 
the JDQZ method. 
A pseudo­code for the JDQZ method is given in Algorithm 2. In order to apply the JDQZ 
method, the user has to supply some parameters, viz. 
''; ø; k max ; j min ; j max ; 
(3.11) 
and a starting vector ~ 
q 2 C n . The parameter '' is a stopping tolerance for the JD iteration 
(the JD iteration will be stopped when krk 2 Ÿ ''), ø is a target (the JDQZ method is supposed 
to compute the eigenvalues closest to ø ), k max is the number of eigenvalues one would like to 
compute, and j min and j max are the minimal (after restart) and maximal dimension of the 
search space spanfV g, respectively. If there is no approximate eigenvector known, one can 
take a random starting vector ~ 
q. The notation Q k = [Q k 
= I , Q 1 = [Q 0 ; ~ 
q] = ~ 
q etc. In case 
k max = 1, Algorithm 2 reduces to the JD method (see Section 3.1). 
The JDQZ method may converge very fast, even for interior eigenvalues and double eigen­ 
values (see [5]). In Section 3.1 we mentioned that the JD method converges quadratically, 

8 Jos L.M. van Dorsselaer 
k=1; j=1; 
choose a starting vector ~ 
q 
Ÿ 0 = 
i 
1 + jø j 2 
j 
j 2 + j ~ 
fij 2 = 1) 
select h~ff; ~ 
fii and u; ~ q = V u 
~ z = Ÿ 0 A~q + Ÿ 1 B~q 
~ 
q = ~ q=k~qk 2 and ~ z = ~ z=k~zk 2 
r = ( ~ 
fiA 
min and W = WZ 1;j min 
j = j min 
else 
obtain v from solving (3.10) (or (3.7)) approximately 
w = Ÿ 0 Av + Ÿ 1 Bv 
expand V with v and W with w; make V and W orthogonal 
j = j + 1 
end if 
solve the eigenvalue problem ~ 
fiW \Lambda AV u = ~ 
ffW \Lambda BV u (j ~ 
ffj 2 + j ~ 
fij 2 = 1) 
select h~ff; ~ 
fii and u; ~ 
q = V u 
~ z = Ÿ 0 A~q + Ÿ 1 B~q 
~ q = ~ 
q=k~qk 2 and ~ 
z = ~ z=k~zk 2 
r = (I 
k = k + 1 
if k ! k max 
V = V Q 2;j and W = WZ 2;j 
j = j 
1 B~q 
~ q = ~ 
q=k~qk 2 and ~ 
z = ~ z=k~zk 2 
r = (I 
Algorithm 2. JDQZ: kmax eigenvalues close to a target ø are computed. '' is a stopping 
tolerance. j min and j max determine the smallest (after restart) and largest dimension of 
the search subspace spanfV g, respectively. 

Computing eigenvalues in continuation methods with the JDQZ method 9 
if the correction equation (3.3) (or (3.7), (3.10) in the JDQZ setting) is solved accurately. 
Experiments show (see, e.g., [5]) that it is not necessary to solve (3.10) accurately in the 
beginning of the JDQZ process for obtaining fast convergence. The correction equation (3.10) 
can be solved iteratively, which allows the JDQZ method to be applicable to very large 
eigenvalue problems. An important question is how to choose a good stopping criterion for 
the iterative solution of (3.10); solving (3.10) accurately may reduce the number of steps in 
the JDQZ method, but the execution time of a single step may become higher. It is not 
clear which strategy leads to the best overall performance of the JDQZ method. In [5] it is 
suggested to solve (3.10) with a Krylov subspace method, and to stop the iterative proces 
when 
kr i k 2 ! 2 
to [5] for more details, discussions, variants, and illustrations of the JDQZ 
method. 
Remark 3.1. For standard eigenvalue problems (i.e. B = I in (3.1)) one can simplify the 
method described above. In [5], the Jacobi­Davidson QR (JDQR) method is proposed for 
computing a partial Schur form AQ k = Q k R k (here R k is a k \Theta k upper triangular matrix, 
and Q k is as above). Roughly speaking, this JDQR method can be obtained from the JDQZ 
method by replacing W by V , Z k by Q k , and h~ff; ~ 
fii by ~ 
– = ~ 
ff= ~ 
fi. Hence, in the JDQR method 
one might save both computation time and memory storage in comparison with the JDQZ 
method. See [5] for more details. 
3.3 Using JDQZ in continuation methods 
In continuation methods we have to compute the rightmost eigenvalues to (2.3), and investi­ 
gate whether an eigenvalue equals zero or not (see Section 2.2). In many physical applications, 
most of the eigenvalues have negative real parts, and only a few of them (if any) have a non­ 
negative real part. In this paper we will consider this situation, which means that eigenvalues 
close and equal to zero belong to the rightmost ones. We now discuss how the parameters 
(3.11) in the JDQZ algorithm should be chosen; we will focus on the target ø (note that A is 
the Jacobian matrix F 0 (y; ¯)). 
Since we are interested in eigenvalues close to 0, one might choose ø = 0. However, it 
may be safer to have ø in the right half plane, in order to avoid missing eigenvalues with 
positive real parts. On the other hand, when jø j is too large, the JDQZ method may not 
be able to designate between different eigenvalues close to 0, which may lead to a slower 
convergence rate (or even no convergence at all). Hence the choice ø = 1 seems reasonable. 
Moreover, in case ø = 1, the approximated eigenvalues h~ff; ~ 
fii closest to ø correspond to the in 
modulus largest approximated eigenvalues of the matrix (A 
the JDQZ method in 
continuation methods. 
The choice of k max depends on the number of bifurcation points one expects to compute; 
it is advised to take k max slightly larger than this number. It is possible to change k max during 
the continuation process. Standard choices for the other parameters can be found in [5] and 
Section 4.2. 

10 Jos L.M. van Dorsselaer 
In general one uses a randomly chosen vector ~ 
q as a starting vector for the JDQZ method. 
But, when using JDQZ in combination with continuation, one has already computed Schur 
vectors in the previous continuation step(s). Using these Schur vectors may lead to faster 
convergence (cf. the selection of a starting vector for computing stationary solutions in 
Algorithm 1). One might e.g. take the Schur vector which was computed first in the previous 
continuation step (this Schur vector is also an eigenvector) as a starting vector ~ 
q. Instead of 
starting JDQZ with one single vector, one might also start with a search subspace spanfV g, 
e.g. the space spanned by the Schur vectors from the previous continuation step. Note that 
it is not necessary to compute eigenvectors (apart from the first one, which is also a Schur 
vector); in fact, the columns of V have to be orthogonal. (Also extrapolation of Schur vectors 
from different continuation steps is possible.) Although these kind of starting procedures 
look attractive, we observed in our experiments that there is not much difference (in CPU­ 
time) between starting with an arbitrary vector ~ 
q or the first Schur vector of the previous 
continuation step. Starting with the subspace spanned by the old Schur vectors even leads to 
higher computation times. See Section 4.2.2 for our experiments and more discussion. 
A possible disadvantage of using previously computed Schur vectors is that they may be 
close to Schur vectors in the current step which do not correspond to the rightmost eigen­ 
values. This may slow down the convergence of the JDQZ method (because `wrong' Schur 
vectors are selected first), or, the JDQZ method may converge to undesired eigenvalues (those 
corresponding to Schur vectors close to the Schur vectors used for starting JDQZ). To illus­ 
trate this, consider, e.g., the 2 \Theta 2 matrices A = diag( 
is a Schur vector for all ¯ 2 R, but corresponds only to the rightmost eigenvalue for 
¯ ! 
4 An application: Rayleigh­B'enard convection 
4.1 The Rayleigh­B'enard problem 
In order to illustrate the JDQZ method in combination with continuation methods, we con­ 
sider the 2D Rayleigh­B'enard problem, which has been studied extensively in the literature 
(see e.g. [3]). A liquid layer in a two­dimensional rectangular box, with length 10 and height 
1, is heated from below. The temperature on the top and bottom of the box is constant, the 
sidewalls are isolated and all velocities are zero on the boundaries (no­slip condition). The 
horizontal and vertical velocity are denoted by u and w respectively, p stands for the (scaled) 
pressure, and the temperature is denoted by T . This leads to the following system of partial 
differential equations: 
Pr 
z 
' 
= 
(4.3) 
@T 
@t 
+ u 
@T 
@x 
+ w 
@T 
@z 
= r 2 T ; 
(4.4) 

Computing eigenvalues in continuation methods with the JDQZ method 11 
with boundary conditions: 
u = w = @T=@x = 0 at x = 0; 10; 
u = w = 0 and T = 1 at z = 0; 
(4.5) 
u = w = T = 0 at z = 1: 
Here Pr is the Prandtl number and Ra is the Rayleigh number; in this paper we take Pr = 5:5, 
and Ra will be our continuation parameter ¯ [3]. Note that p is not uniquely determined; this 
leads to, after discretization of the spatial variables, a Jacobian which is always singular. It 
is clear from Section 2 that this is not attractive for continuation, and therefore we prescribe 
p at a certain point: we set p(5; 1 
2 ) = 0. 
The equations (4.1) -- (4.5) are discretized on a staggered grid with uniform mesh sizes, 
using finite difference approximations. For the nonlinear terms we use first order upwind, and 
the other terms are discretized by second order central differences. In the grid cell containing 
the point (5; 1 
2 ) we replace the discretization of (4.3) by p(5; 1 
2 ) = 0. The discretized system 
obtained in this way can be written in the form (2.1), where Y (t) contains the velocities u, 
w, the pressure p and the temperature T at certain gridpoints, and ¯ = Ra. The dimension 
n of the system (2.1) equals 4n x n z , where n x and n z stand for the number of grid cells in the 
horizontal and vertical direction, respectively. The matrix B is a diagonal matrix, which is 
singular due to the boundary conditions (4.5) and the absence of time derivatives in (4.3). The 
unknowns are numbered per grid cell, and a grid cell is only coupled to six of its neighbours. 
The grid cells are ordered by column, from bottom to top, beginning with the first column. 
This leads to a Jacobian matrix of which the bandwidth is equal to 4n z + 3; when n z is small, 
it is feasible to compute LU­factorizations of this Jacobian F 0 (y; Ra) (needed for the Newton 
process in Algorithm 1), and the preconditioner K = A 
'enard problem was formulated using the temperature 
T , the streamfunction / and the vorticity ! (with u = @/=@z, w = 
reason for using primitive variables is that this approach can easily be extended 
to 3D, while this is not possible for the streamfunction­vorticity formulation. 
4.2 Numerical experiments 
In this section we will describe our experiments. In Section 4.2.1 we will give the actual results 
of our continuation code, and the performance of the JDQZ method (for different choices of 
parameters) will be discussed in Section 4.2.2. 
4.2.1 The continuation code applied to the 2D Rayleigh­B'enard problem 
In Table 1 we have listed the first two bifurcation points obtained with our code (for different 
grid sizes), and we compared these to the corresponding values from [3]. 
From Table 1 we see that n x = 129 and n z = 33 gives the most accurate results. Due 
to memory limitations, we were not able to perform experiments with smaller grid sizes. In 
particular, the memory requirements for the LU­factorization of A 
12 Jos L.M. van Dorsselaer 
n z = 17 suggest that is not useful to take n x ? 129. We took n x = 129 and n z = 33 in the 
experiments described in this section. 
n x n z Ra 1 Ra 2 
129 17 1698.3 1701.7 
257 17 1698.8 1701.8 
129 33 1720.0 1724.0 
[3] 1731.3 1734.1 
Table 1. Values of the first two bifurcation points (Ra 1 and Ra 2 ) for different grid sizes, 
and the corresponding values obtained in [3]. 
The following bifurcation behaviour for the 2D Rayleigh­B'enard problem has been found in 
[3]. The trivial, motionless, solution u j w j 0 and T j 1 
is stable for Ra ? Ra 1 , and the motionless stationary solution becomes unstable for 
Ra ? Ra 1 . At Ra = Ra 2 an unstable 9­cell solution branches off, and this solution splits into 
a stable and an unstable stationary solution at another bifurcation point. The motionless 
solution splits again at a third bifurcation point, which results in an 11­cell solution. 
1700 1720 1740 1760 1780 1800 1820 1840 1860 1880 1900 
0 
0.05 
0.1 
0.15 
0.2 
0.25 
10-cell 
11-cell 
9-cell 
Figure 1. A bifurcation diagram. The Rayleigh number Ra is on the horizontal axis, and 
the vertical velocity w at (0:0362; 0:75) (the same point as in [3]) is on the vertical axis. 
Each curve corresponds to a stationary solution; stable stationary solutions are indicated 
with a solid line, and unstable ones with a dashed line. The bifurcation points are given 
by small circles (ffi). 
With our code we found the same bifurcation behaviour. The 9­cell solution splits at Ra = 
1864:6, and the 11­cell solution branches off at Ra = 1781:0. Our results are visualized in 

Computing eigenvalues in continuation methods with the JDQZ method 13 
Figure 1, and three nontrivial stationary solutions are displayed in Figure 2. The pictures in 
Figure 2 look similar to the corresponding ones in [3]. Hence we may conclude that our code 
produces meaningful results. 
0 1 2 3 4 5 6 7 8 9 10 
0 
0.5 
1 
0 1 2 3 4 5 6 7 8 9 10 
0 
0.5 
1 
0 1 2 3 4 5 6 7 8 9 10 
0 
0.5 
1 
Figure 2. Three solutions, viz. a 10­cell solution at Ra = 1721 (the first picture), a 
9­cell solution at Ra = 1725 (the second picture), and an 11­cell solution at Ra = 1782. 
In each picture, the current of the fluid is given; a dashed curve means that the fluid 
moves clockwise, and a solid curve indicates an anticlockwise direction of the fluid. 
The bifurcation points can be computed as follows. At Ra = Ra 1 , the sign of the largest eigen­ 
value of (2.3) changes (the rightmost eigenvalues of (2.3) turned out to be real in our case), 
and the secant method has been used to determine at which value Ra the largest eigenvalue 
of (2.3) equals 0. To compute the second bifurcation point Ra = Ra 2 , we determined, using 
the secant method, when the second largest eigenvalue equals 0, etc. In order to determine 
more bifurcation points at an unstable branch, we need to compute several eigenvalues at 
each continuation step, and not only the rightmost one. 
4.2.2 Performance of the JDQZ method 
In this section we will consider the performance of the JDQZ method. We will deal with 
the case n x = 129 and n z = 17, instead of n x = 129 and n z = 33, because of memory 

14 Jos L.M. van Dorsselaer 
requirements and slower computation times in the latter case (in particular wall clock time, 
due to swapping). 
In Algorithm 1 we set \Delta¯ = \DeltaRa = 10, and the Newton iteration is stopped when two con­ 
secutive approximations y (k) 
j and y (k+1) 
j of the stationary solution y j satisfy ky (k+1) 
j 
In Algorithm 2 (the JDQZ method) we set 
ø = 1 (cf. our discussion in Section 3.3), j min = 10, j max = 20, and we take different values for 
k max and '' (see Table 2). In order to solve the correction equation (3.10) we used two different 
Krylov subspace methods, viz. GMRESm (at most m steps with full GMRES, no restarts) 
[9] with m = 5, and BiCGstab(`) [12] with ` = 2 and a maximum of 100 matrix­vector multi­ 
plications for solving (3.10) per JD step. Further, the stopping criterion (3.12) has been used 
for both methods. It is likely that (3.10) is solved more accurately with BiCGstab(2) than 
with GMRES 5 (more matrix­vector multiplications are allowed for BiCGstab(2)). Therefore 
one might expect that less JD steps are needed for BiCGstab(2), but, a single JD step is more 
expensive. A priori it is not clear which Krylov subspace method leads to the most efficient 
variant of JDQZ. In order to obtain a subspace spanfV g as soon as possible, GMRES 1 is 
used to solve (3.10) when j Ÿ j min (j is the dimension of spanfV g). 
We consider three different strategies to start the JDQZ method, viz. starting with a 
random vector, starting with the first column of the matrix Q kmax computed in the previous 
continuation step, and starting with the subspace spanned by the columns of Q kmax ; we 
will call these methods JDQZ1, JDQZ2 and JDQZ3 respectively (the latter two have been 
discussed in Section 3.3). 
In Table 2 we listed the CPU­times (in seconds) for an `average' continuation step. (We 
performed 15 steps, and the average of the last 10 steps is listed in the table, so that the 
effect of starting the continuation is ruled out). In each step, 3 Newton iterations were needed 
to compute a stationary solution with Algorithm 1. The computations were done on a SUN 
SPARC 1000E with 4 processors. 
solving CPU(s) CPU(s) CPU(s) 
k max '' (3.12) JDQZ1 JDQZ2 JDQZ3 
289 
6 10 
able 2. The CPU­time (in seconds) for an `average' continuation step. In the first line 
of the table the CPU­time for one step of Algorithm 1 has been displayed (no eigenvalues 
have been computed in this case). 
From the results in Table 2 we may conclude that the JDQZ method is well suited for com­ 
puting several eigenvalues accurately. 

Computing eigenvalues in continuation methods with the JDQZ method 15 
When we compare the CPU­times for GMRES 5 and BiCGstab(2) we see that the first one 
is the most efficient for solving (3.10) when JDQZ1 or JDQZ2 is used. In the experiments 
with JDQZ1 and JDQZ2 we observed that K = A 
`average' step, and this 
might explain why GMRES 5 is more efficient for JDQZ1 and JDQZ2. On the other hand, 
when JDQZ3 is used, BiCGstab(2) turns out to be more efficient than GMRES 5 . In these 
experiments we observed that the average number of JD steps for GMRES 5 is significantly 
larger than for BiCGstab(2), and this might explain why solving (3.10) with BiCGstab(2) is 
more efficient when JDQZ3 is used. 
5 10 15 20 25 30 35 40 45 50 55 
-12 
-10 
-8 
-6 
-4 
-2 
0 
5 10 15 20 25 30 35 40 45 50 55 
-12 
-10 
-8 
-6 
-4 
-2 
0 
Figure 3. The convergence of the JDQZ method with different starting strategies for the 
last continuation step (k max = 6, '' = 10 
krk 2 ) (r is 
the residual in the JDQZ method) is on the vertical axis. In both pictures, the solid curve 
corresponds to JDQZ1; the dashed curve in the upper picture corresponds to JDQZ2, and 
the dashed curve in the lower picture corresponds to JDQZ3. 
A somewhat surprising result is that starting the JDQZ method with Schur vectors from the 

16 Jos L.M. van Dorsselaer 
previous continuation step does not improve the efficiency of the method much. Starting with 
the first Schur vector (JDQZ2) leads to (almost) the same CPU­times (in particular when 
GMRES 5 is used to solve (3.10)), while starting with all Schur vectors (JDQZ3) leads to a 
significantly slower method. In order to understand this, we have plotted the convergence 
behaviour of the different JDQZ methods at the last continuation step in Figure 3. When we 
compare JDQZ1 with JDQZ2 we observe that the first Schur vector is detected earlier with 
JDQZ2, but JDQZ2 needs more JD steps to compute the other Schur vectors. This might 
be explained as follows: the first Schur vector of the previous continuation step might have 
a larger component in the direction of the first Schur vector than a random vector (observe 
that the norm of the initial residual is smaller), so that this first Schur vector is found earlier. 
On the other hand, it is likely that a random vector has larger components in the direction 
of the other Schur vectors, so the subspace spanfV g (after the detection of the first Schur 
vector) might contain more information about the other Schur vectors than the corresponding 
subspace in JDQZ2 (compare the norms of the residuals just after the first restart, when the 
first Schur vector has been removed from spanfV g). In the upper picture in Figure 3 we see 
that both methods need about the same number of JD steps in order to find 6 Schur vectors 
(with corresponding eigenvalues). 
For the JDQZ3 method we would not expect such behaviour, because the start subspace 
contains all Schur vectors computed in the previous continuation step. Again the first Schur 
vector is found earlier (in comparison with JDQZ1). It is possible that some Schur vectors 
of the previous continuation step are removed in spanfV g when this space is reduced, after 
restarting, from dimension 20 (= j max ) to 10 (= j min ), but this does not explain why JDQZ3 
needs more JD steps than JDQZ1 (or JDQZ2) to discover the other Schur vectors. Perhaps 
it is better to start JDQZ with the first old Schur vector, and add the second Schur vector of 
the previous continuation step to spanfV g when the first new Schur vector has been detected 
and removed from spanfV g etc. We have not tried this approach. On the other hand, JDQZ1 
performs very well for this example, so it could be hard to construct a method which performs 
better in this case. 
With both '' = 10 
the imaginary axis has been found. The 4 eigenvalues computed with 
k max = 4 were also detected with k max = 6, and they turned out to be the four rightmost 
ones. This shows that the choice k max = 4 is a reasonable one. 
Unfortunately, it in general is not possible to determine a complete LU­factorization 
of the matrix A 
to 
investigate the behaviour of JDQZ when incomplete factorizations are used, we have repeated 
some of our experiments with an incomplete factorization of A 
equations 
arising from Newton's method in Algorithm 1 have been solved with BiCGstab(8). The CPU­ 
time of an average continuation step with k max = 0 is 388 seconds, while only 33 seconds were 
needed for a direct factorization of the Jacobian (see Table 2). It is not an easy task to compute 
eigenvalues with JDQZ, using this preconditioner; we did not find any eigenvalues with the 
methods and parameters chosen as in Table 2 (i.e. at most 100 matrix­vector multiplications 
with BiCGstab(2)). Using BiCGstab(8) instead, with a maximum of 1000 matrix­vector 

Computing eigenvalues in continuation methods with the JDQZ method 17 
multiplications for solving (3.10) (instead of 100), and k max = 4, '' = 10 
or the 2D Rayleigh­B'enard problem, this incomplete factorization 
is not well suited as a preconditioner for solving the correction equation (3.10). Without a 
good preconditioner it can be very hard to obtain eigenvalues with the JDQZ method within 
a reasonable computation time. 
Finally we compare the JDQZ method to the SIT method [13], which has been used 
in [3] to compute the rightmost eigenvalues. In [3] the SIT method (which is essentially a 
block version of the power method, see [13, 3] for the details) has been applied to the matrix 
C = (B 
as not able to compute the rightmost eigenvalue accurately; in our 
experiments we observed that the error in the rightmost eigenvalue is slightly more than 1% --- 
even for the experiments with much higher CPU­times for an `average' continuation step than 
those corresponding to the JDQZ method with k max = 4 and '' = 10 
eigenvalues occurring in continuation methods for problems of the type 
that we have considered. 
4.3 Conclusions 
In this section the main conclusions of the paper are summarized. 
The JDQZ method can be a very efficient tool for computing several rightmost eigenvalues 
in continuation methods, provided a good preconditioner for the correction equation (3.10) 
is available. Without such a preconditioner the JDQZ method may behave rather poorly (cf. 
also [5]). For some (small) problems a good preconditioner can be constructed by a direct 
factorization of A 
t than the SIT 
method [13], which is not surprising because the SIT method is based on a block version of 
the power method, which converges linearly, while the JDQZ method often shows a quadratic 
convergence behaviour. 
Using one or more Schur vectors from the previous continuation step does not necessarily 
lead to a faster convergence of the JDQZ method. Starting the JDQZ method with one 
Schur vector or a random vector gives about the same computation times, while starting 
with a subspace containing all previously computed Schur vectors led to a significantly slower 
method. 
Acknowledgements. The author wishes to thank D.R. Fokkema and M.B. van Gijzen 
(Utrecht University, Dept. of Math.) for providing the JDQZ code, G. Tiesinga (Groningen 
University, Dept. of Math.) for the Rayleigh­B'enard code with the incomplete factorization 
mentioned in Section 4.2.2, and M. Vellinga, H.A. Dijkstra, and M.J. Molemaker (Utrecht 
University, IMAU) for the continuation code as used in [3]. Further he is indebted to G.L.G. 
Sleijpen, H.A. van der Vorst (Utrecht University, Dept. of Math.), and F.W. Wubs (Groningen 

18 Jos L.M. van Dorsselaer 
University, Dept. of Math.) for stimulating discussions, and suggestions concerning the 
presentation of this paper. 
References 
[1] C.W. Brand: An incomplete­factorization preconditioning using repeated red­black ordering. 
Numer. Math. 61, 433­454 (1992). 
[2] P.N. Brown & Y. Saad: Hybrid Krylov methods for nonlinear systems of equations. SIAM J. 
Sci. Stat. Comput. 11, 450­481 (1990). 
[3] H.A. Dijkstra, M.J. Molemaker, A. van der Ploeg & E.F.F. Botta: An efficient code 
to compute non­parallel steady flows and their linear stability. Computers and Fluids 24, 415­434 
(1995). 
[4] D.R. Fokkema, G.L.G. Sleijpen & H.A. van der Vorst: Accelerated inexact Newton 
schemes for large systems of nonlinear equations. Preprint nr. 918, Dept. Math., Utrecht Univer­ 
sity, Utrecht, July 1995. To appear in SIAM J. Sci. Comput. 
[5] D.R. Fokkema, G.L.G. Sleijpen & H.A. van der Vorst: Jacobi­Davidson style QR and 
QZ algorithms for the partial reduction of matrix pencils. Preprint nr. 941, Dept. Math., Utrecht 
University, Utrecht, January 1996. To appear in SIAM J. Sci. Comput. 
[6] G.H. Golub & C.F. van Loan: Matrix Computations, 2nd Ed. Baltimore: The John Hopkins 
University Press, 1989. 
[7] G.H. Golub & H.A. van der Vorst: 150 Years old and still alive: eigenproblems. Preprint 
nr. 981, Dept. Math., Utrecht University, Utrecht, October 1996. To appear in The State of the 
Art in Numerical Analysis (I.S. Duff & G.A. Watson eds). Oxford: Oxford University Press. 
[8] H.B. Keller: Numerical solution of bifurcation and nonlinear eigenvalue problems, in Applica­ 
tions of Bifurcation Theory (P.H. Rabinowitz ed.), 359­384. Academic Press, 1977. 
[9] Y. Saad & M. Schultz: GMRES: A generalized minimal residual algorithm for solving non­ 
symmetric linear systems. SIAM J. Sci. Stat. Comput. 7, 856­869 (1986). 
[10] R. Seydel: Practical Bifurcation and Stability Analysis, 2nd Ed. New York: Springer, 1994. 
[11] G.L.G. Sleijpen, A.G.L. Booten, D.R. Fokkema & H.A. van der Vorst: Jacobi­ 
Davidson type methods for generalized eigenproblems and polynomial eigenproblems. BIT 36, 
595­633 (1996). 
[12] G.L.G. Sleijpen & D.R. Fokkema: BiCGstab(`) for linear equations involving unsymmetric 
matrices with complex spectrum. Electronic Transactions in Numerical Analysis (ETNA) 1, 11­32 
(1993). 
[13] W.J. Stewart & A. Jennings: A simultaneous iteration algorithm for real matrices. ACM 
Trans. Math. Softw. 7, 184­198 (1981). 
[14] G. Tiesinga: Block preconditioned BiCGstab(2) for solving the Navier­Stokes equations. Z. 
angew. Math. Mech. 76, Supplement 1, 563­564 (1996).