geneig¶

Computes the generalized eigenexpansion of a system \(Ax = \lambda Bx\), with A and B real.

Synopsis¶

geneig (a, b, alpha, beta)

Required Arguments¶

float a[[]] (Input): Array of size n × n containing the coefficient matrix A.
float b[[]] (Input): Array of size n × n containing the coefficient matrix B.
complex alpha (Output): Vector of size n containing scalars \(\alpha_i\). If \(\beta_i \neq 0\), \(\lambda_i = \alpha_i / \beta_i\) for i = 0, …, n − 1 are the eigenvalues of the system.
float beta (Output): Vector of size n .

Optional Arguments¶

vectors (Output): An array of size n × n containing eigenvectors of the problem. Each vector is normalized to have Euclidean length equal to the value one.

Description¶

The function geneig uses the QZ algorithm to compute the eigenvalues and eigenvectors of the generalized eigensystem \(Ax = \lambda Bx\), where A and B are real matrices of order n. The eigenvalues for this problem can be infinite, so α and β are returned instead of λ. If β is nonzero, \(\lambda = \alpha / \beta\).

The first step of the QZ algorithm is to simultaneously reduce A to upper-Hessenberg form and B to upper-triangular form. Then, orthogonal transformations are used to reduce A to quasi-upper-triangular form while keeping B upper triangular. The generalized eigenvalues and eigenvectors for the reduced problem are then computed.

The function geneig is based on the QZ algorithm due to Moler and Stewart (1973), as implemented by the EISPACK routines QZHES, QZIT and QZVAL; see Garbow et al. (1977).

Examples¶

Example 1¶

In this example, the eigenvalue, λ, of system \(Ax = \lambda Bx\) is computed, where

\[\begin{split}A = \begin{bmatrix} 1.0 & 0.5 & 0.0 \\ -10.0 & 2.0 & 0.0 \\ 5.0 & 1.0 & 0.5 \\ \end{bmatrix} \text{ and } B = \begin{bmatrix} 0.5 & 0.0 & 0.0 \\ 3.0 & 3.0 & 0.0 \\ 4.0 & 0.5 & 1.0 \\ \end{bmatrix}\end{split}\]

from __future__ import print_function
from numpy import *
from pyimsl.math.geneig import geneig
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

a = [[1., 0.5, 0.],
     [-10., 2., 0.],
     [5., 1., 0.5]]
b = [[0.5, 0., 0.],
     [3., 3., 0.],
     [4., 0.5, 1.]]
alpha = []
beta = []

# Compute eigenvalues
eval = geneig(a, b, alpha, beta)
eval = zeros(len(alpha), dtype='complex64')
for i in range(0, len(alpha)):
    if (beta[i] != 0.0):
        eval[i] = alpha[i] / beta[i]
    else:
        print("Infinite eigenvalue at: ", i)

# Print eigenvalues
writeMatrixComplex("Eigenvalues", eval)

Output¶

 
                     Eigenvalues
                        1                          2
(      0.833,      1.993)  (      0.833,     -1.993)
 
                        3
(      0.500,      0.000)

Example 2¶

This example finds the eigenvalues and eigenvectors of the same eigensystem given in the last example.

from __future__ import print_function
from numpy import *
from pyimsl.math.geneig import geneig
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

a = [[1., 0.5, 0.],
     [-10., 2., 0.],
     [5., 1., 0.5]]
b = [[0.5, 0., 0.],
     [3., 3., 0.],
     [4., 0.5, 1.]]
alpha = []
beta = []
evec = []

# Compute eigenvalues and eigenvectors
eval = geneig(a, b, alpha, beta, vectors=evec)
eval = zeros(len(alpha), dtype='complex64')
for i in range(0, len(alpha)):
    if (beta[i] != 0.0):
        eval[i] = alpha[i] / beta[i]
    else:
        print("Infinite eigenvalue at: ", i)

# Print eigenvalues and eigenvectors
writeMatrixComplex("Eigenvalues", eval)
writeMatrixComplex("Egenvectors", evec)

Output¶

 
                     Eigenvalues
                        1                          2
(      0.833,      1.993)  (      0.833,     -1.993)
 
                        3
(      0.500,     -0.000)
 
                      Egenvectors
                           1                          2
1  (     -0.197,      0.150)  (     -0.197,     -0.150)
2  (     -0.069,     -0.568)  (     -0.069,      0.568)
3  (      0.782,      0.000)  (      0.782,      0.000)
 
                           3
1  (     -0.000,      0.000)
2  (     -0.000,      0.000)
3  (      1.000,      0.000)