I am trying to solve the Orr-Sommerfeld equation numerically, using the techniques given in this article. This leads to solving a generalized eigenvalue problem, that is, given two matrices $\mathbf A,\mathbf B$ with possibly complex entries, we would like to find a pair $(\lambda,\boldsymbol v)$ where both $\lambda $ and $\boldsymbol v$ are possibly complex, satisfying $$\mathbf A\boldsymbol v=\lambda \mathbf B\boldsymbol v\tag 1$$
Note that $\mathbf B$ is not necessarily invertible, i.e, we cannot (in general) simply multiply both sides by $\mathbf B^{-1}$ to reduce it to an ordinary eigenvalue problem. And, even if we could, these matrices $\mathbf A,\mathbf B$ are fairly large, say about $100\times 100$ (and $\mathbf B$ is, in practice, probably quite sparse), so this would probably be a waste of computational resources anyway. (I would also like to add that in my application I am more interested in the eigenvalue $\lambda$ than I am in the eigenvector $\boldsymbol v$.)
A common way of solving this kind of generalized eigenvalue problem is with generalized Schur decomposition, also known as "QZ decomposition" in which we write
$$\mathbf A=\mathbf {QS}\mathbf Z^* \\ \mathbf B=\mathbf {QT}\mathbf Z^*\tag {2}$$ Where $\mathbf Q,\mathbf Z$ are unitary and $\mathbf S,\mathbf T$ are upper triangular. (The asterisk denotes the conjugate transpose.) The eigenvalues of $(1)$ then can be given as the ratio of the diagonal entries of these two upper triangular matrices, $$\lambda_i=\frac{T_{ii}}{S_{ii}}$$
Solving the problem numerically in python.
Python's scipy
package has the module linalg.eig
, which, according to the documentation, is able to solve generalized eigenvalue problems of the form of $(1)$. The right-hand-side matrix can be given in the optional argument b
. The default value of this argument b
is None
, in which the module will solve the standard eigenvector problem, i.e $\mathbf B=\mathbf I$ in $(1)$.
However, the scipy
package also includes the linalg.qz
module, which, according to the documentation is able to perform a decomposition of the form of $(2)$.
My question, is, would I be better off simply inputting my matrices $\mathbf A,\mathbf B$ as the a
and b
arguments into the scipy.linalg.eig
module, or would I be better off first finding their Schur decompositions with scipy.linalg.qz
and then inputting the matrices $\mathbf S,\mathbf T$ into scipy.linalg.eig
separately, and setting the argument b=None
? Which one is faster?
Thanks.