Quadratic eigenvalue problem - biorthogonal system

Question

I need to solve the following quadratic eigenvalue problem: $$\left(-\lambda^2\lbrack \mathbf{X}\rbrack+\text{i}\lambda\lbrack \mathbf{Y}\rbrack+\lbrack \mathbf{Z}\rbrack\right)\left(\mathbf{u}\right)=\left(\mathbf{0}\right)$$ To solve this I used a linearization to transform this quadratic eigenvalue problem into the following generalized eigenvalue problem: $$\lbrack \mathbf{A}\rbrack(\mathbf{u})=\lambda\lbrack \mathbf{B}\rbrack(\mathbf{u})$$ with $\lbrack \mathbf{A}\rbrack$ and $\lbrack \mathbf{B}\rbrack$ respectively: $$\lbrack \mathbf{A}\rbrack=\begin{bmatrix}\text{i}\lbrack \mathbf{Y}\rbrack&\lbrack \mathbf{Z}\rbrack\\-\lbrack \mathbf{I}\rbrack&\lbrack \mathbf{0}\rbrack\end{bmatrix}~~~~~~~~~\lbrack \mathbf{B}\rbrack=\begin{bmatrix}-\lbrack \mathbf{X}\rbrack&\lbrack \mathbf{0}\rbrack\\\lbrack \mathbf{0}\rbrack&\lbrack \mathbf{I}\rbrack\end{bmatrix}$$ Due to the circumstance that $\lbrack \mathbf{A}\rbrack$ and $\lbrack \mathbf{B}\rbrack$ aren't hermitian in general in need to solve the adjoint eigenvalue problem $(\lbrack \mathbf{A}\rbrack^*)^\text{T}(\mathbf{v})=\lambda^*(\lbrack \mathbf{B}\rbrack^*)^\text{T}(\mathbf{v})$ aswell to receive a biorthogonal system.

• $\lbrack \mathbf{X}\rbrack, \lbrack \mathbf{Y}\rbrack, \lbrack \mathbf{Z}\rbrack$ - random $n\times n$ matrices
• $\lbrack \mathbf{A}\rbrack, \lbrack \mathbf{B}\rbrack$ - corresponding matrices based on the linearization ($2n\times 2n$)
• $\lbrack \mathbf{I}\rbrack$ - Identity Matrix ($n\times n$)
• $\lambda$ - eigenvalue
• $(\mathbf{u}), (\mathbf{v})$ - right- and left-eigenvectors

This is my implementation in Mathematica:

SeedRandom[123];

n = 15;
XMat = RandomReal[{1*10^-7, 1*10^-6}, {n, n}];
YMat = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];
ZMat = RandomReal[{1*10^-12, 1*10^-11}, {n, n}];

AMat = Join[Join[I*YMat, ZMat, 2], Join[-IdentityMatrix[n], ConstantArray[0, {n, n}], 2]];
BMat = Join[Join[-XMat, ConstantArray[0, {n, n}], 2], Join[ConstantArray[0, {n, n}], IdentityMatrix[n], 2]];

{EWR, EVR} = Eigensystem[{AMat, BMat}];
{EWL, EVL} = Eigensystem[{ConjugateTranspose[AMat], ConjugateTranspose[BMat]}];

perm = Flatten[Nearest[EWL -> "Index", EWR]];
EWR = EWR[[perm]];
EVR = EVR[[perm]];

My main problem is that the corresponding eigenvalues of the mentioned eigenvalue problems are the same and they should be complex conjugated. Did I make a mistake in my Mathematica code or is there a mistake in my theoretical approach?

Update:

I corrected the Mathematica code and changed Transpose to ConjugateTranspose but this didn't solve the problem for all eigenvalues. As shown below most of the corresponding eigenvalues are still identical.

0.00456861 + 3.28732*10^-6 I    0.00456861 - 3.2873*10^-6 I
-0.00456861 + 3.28731*10^-6 I   -0.00456861 - 3.28733*10^-6 I
-0.00257144 + 0.00316279 I      -0.00257543 + 0.00316029 I
0.00257144 + 0.00316279 I       0.00257543 + 0.00316029 I
0.00257543 - 0.00316029 I       0.00257144 - 0.00316279 I
-0.00257543 - 0.00316029 I      -0.00257144 - 0.00316279 I
-0.00372817 + 0.00101919 I      -0.00374614 + 0.00102208 I
0.00372817 + 0.00101919 I       0.00374614 + 0.00102208 I
0.00374614 - 0.00102208 I       0.00372817 - 0.00101919 I
-0.00374614 - 0.00102208 I      -0.00372817 - 0.00101919 I
0.002509 + 0.00236923 I         0.00250875 + 0.00242257 I
-0.002509 + 0.00236923 I        -0.00250875 + 0.00242257 I
0.00250875 - 0.00242257 I       0.002509 - 0.00236923 I
-0.00250875 - 0.00242257 I      -0.002509 - 0.00236923 I
-0.00323767 - 4.43767*10^-6 I   -0.00323767 + 4.43769*10^-6 I
0.00323767 - 4.43767*10^-6 I    0.00323767 + 4.43767*10^-6 I
-0.00237108 - 0.00201776 I      -0.00236452 - 0.00204887 I
0.00237108 - 0.00201776 I       0.00236452 - 0.00204887 I
0.00236452 + 0.00204887 I       0.00237108 + 0.00201776 I
-0.00236452 + 0.00204887 I      -0.00237108 + 0.00201776 I
1.82496*10^-11 - 0.00262931 I   -5.2503*10^-12 - 0.00262943 I
-2.53022*10^-14 + 0.00262943 I  3.85978*10^-17 + 0.00262931 I
0.00196643 + 0.00104558 I       0.00196388 + 0.00105333 I
-0.00196643 + 0.00104558 I      -0.00196388 + 0.00105333 I
0.00196388 - 0.00105333 I       0.00196643 - 0.00104558 I
-0.00196388 - 0.00105333 I      -0.00196643 - 0.00104558 I
-0.000751694 - 0.00133752 I     -0.000759215 - 0.00134487 I
0.000751694 - 0.00133752 I      0.000759215 - 0.00134487 I
0.000759215 + 0.00134487 I      0.000751694 + 0.00133752 I
-0.000759215 + 0.00134487 I     -0.000751694 + 0.00133752 I

Answer 1

This answer is similar to the solution you have found but has some advantages. Your equation often occurs with second order differential equations such as the equations for a vibrating system

$$ X u''+Y u'+ Z u=0 $$

Here the matrices X, Y and Z are the mass, damping and stiffness matrices. In many dynamic cases they are symmetric matrices which yours may or may not be. A good way of organising the eigenvalue problem is to write as a system, as you have done, but in an alternative manner as follows.

$$ \left( \begin{array}{cc} -Z & 0 \\ 0 & X \\ \end{array} \right) \left( \begin{array}{c} U' \\ U'' \\ \end{array} \right)+\left( \begin{array}{cc} 0 & Z \\ Z & Y \\ \end{array} \right) \left( \begin{array}{c} U \\ U' \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) $$

This form has the advantage of keeping the symmetry if the individual matrices are symmetric. This is important since it avoids left and right eigenvectors.

Here is an example. First make the three individual matrices, here symmetric.

n = 3; (* Number of rows*)
a = RandomReal[{-1, 1}, {n, n}];
x = a + Transpose[a];
a = RandomReal[{-1, 1}, {n, n}];
y = a + Transpose[a];
a = RandomReal[{-1, 1}, {n, n}];
z = a + Transpose[a];
zero = 0 x;

Then we form the symmetric matrices and solve the eigen problem

a = ArrayFlatten[{{-z, zero}, {zero, x}}];
b = ArrayFlatten[{{zero, z}, {z, y}}];
{vals, vecs} = Eigensystem[{b, -a}];

As we should we get diagonal matrices

vecs . a . Transpose[vecs] // Chop // MatrixForm
vecs . b . Transpose[vecs] // Chop // MatrixForm

Hope that helps.