# What is mathematica command for finding eigenvalues of a generalized eigenvalue problem?

The provided code below is designed to compute the eigenvalues for the problem I inquired about earlier in this post: How to find eigenvalues?. I believe my code is mostly accurate, with the exception of the final line. Could someone please assist in determining the correct command for solving the generalized eigenvalue problem, [Lambda]UX = WX, where W is a non-singular, invertible matrix? Matlab includes the QZ algorithm for addressing such problems. What is the Mathematica counterpart?

ClearAll["Global*"];
ClearAll[recurrenceTables];
recurrenceTables[x_?VectorQ, n_] := Module[{T, Tprime}, T = ConstantArray[0., {n + 1}];
Tprime = ConstantArray[0., {n + 1}];
T[[1]] = ConstantArray[1., Length[x]];
T[[2]] = x;
Tprime[[1]] = ConstantArray[0., Length[x]];
Tprime[[2]] = ConstantArray[1., Length[x]];
Do[T[[i]] = 2. x T[[i - 1]] - T[[i - 2]];
Tprime[[i]] = 2. T[[i - 1]] + 2. x Tprime[[i - 1]] - Tprime[[i - 2]];, {i, 3, n + 1}];
{Transpose[T], Transpose[Tprime]}];

n = 5;(* Collocation points*)
B = 0.11; ita = 1.0; F = 0.28; L = 10;(*parameters*)

(*Definig Collocation points*)
y = (L/2) + (L/2)*Cos[N[Pi/n] Range[0., n]];
{T, Tprime} = recurrenceTables[y, n];

(*Defining all the functions we need*)
Z[x_] := (L*(x + 1.0))/2.0;
v[Z_] := Exp[(F/4.0 ita) Z];
a[Z_] := Exp[-(F/4.0 ita) Z];
h[Z_] := Exp[-(F/4.0 ita) Z] Sqrt[((4.0 B)/F) + (1 - ((4.0 B)/F)) Exp[(F/4.0 ita) Z]];

Zy = Map[Z, y];
Vbar = DiagonalMatrix[SparseArray[Map[v, Zy]]];
Vbarprime = DiagonalMatrix[SparseArray[Map[v', Zy]]];

Abar = DiagonalMatrix[SparseArray[Map[a, Zy]]];
Abarprime = DiagonalMatrix[SparseArray[Map[a', Zy]]];

Hbar = DiagonalMatrix[SparseArray[Map[h, Zy]]];
Hbarprime = DiagonalMatrix[SparseArray[Map[h', Zy]]];

(*creating entries of the final U and W matrix below*)
A1 = 2.0 Vbar . Tprime + Vbarprime . T;
V1 = 2.0 Abar . Tprime + Abarprime . T;
T2 = Abar . T;
V2 = Hbar . Abarprime . T - 2.0 Abar . Hbar . Tprime;
A2 = 2.0 Vbar . Hbar . Tprime - Hbar . Vbarprime . T;
H2 = Vbar . Abarprime . T - Abar . Vbarprime . T + (F/(2.0*ita)) T;
T31 = Hbar . Hbar . T;
T32 = 2.0 Abar . Hbar . T;
F2 = -(1/(2.0 *ita)) T;
H3 = 2.0  Vbar . Abarprime . Hbar . T - 2.0 Abar . Vbar . Hbarprime . T - 4.0 Abar . Vbar . Hbar . Tprime;
V3 = Hbar . Hbar . Abarprime . T - 2.0 Abar . Hbar . Hbarprime . T;
A3 = 2.0 Vbar . Hbar . Hbar . Tprime - 2.0 Vbar . Hbar . Hbarprime . T - (2.0/ita) B Abar . T;

(*final generalized eigenvalue problem \[Lambda]UX=WX *)
U = ArrayFlatten[{{T, 0. T, 0. T}, {0.0, 0.0, T2}, {T31, 0.0, -T32}}];
W = ArrayFlatten[{{-A1, -V1, 0.0}, {-A2, -V2, -H2}, {-A3, -V3, -H3}}];
Vals = Eigenvalues[{W, U}]

• Help says Eigenvalues[{m,a}] gives the generalized eigenvalues of m with respect to a. I do not know if this applies to what you have. Commented Aug 23, 2023 at 5:43
• Thank you so much but I could not understand how to apply this document of WOLFRAM for generalized eigenvalues @Nasser:reference.wolfram.com/applications/insydes/ReferenceManual/… Commented Aug 23, 2023 at 6:30
• You are looking at some old documentation of some specific package/application it seems. I meant the main documentation Eigenvalues.html see the second entry at the top. Commented Aug 23, 2023 at 6:34
• Okay, I see. Thanks for the clarification @Nasser Commented Aug 23, 2023 at 6:40
• @Fun123 You solution is valid for L=1 only. To extend it to L=10 use Z[x_] :=L ((x + 1.0))/2.0. Unfortunately, matrix W is not invertible. Therefore we can try to solve this problem with NMInimize as in my answer on mathematica.stackexchange.com/questions/289044/… Commented Aug 23, 2023 at 10:11

There are several typos in the code with matrix U,W definition. After cleaning, the code looks like this

ClearAll["Global*"];
ClearAll[recurrenceTables];
recurrenceTables[x_?VectorQ, n_] :=
Module[{T, Tprime}, T = ConstantArray[0., {n + 1}];
Tprime = ConstantArray[0., {n + 1}];
T[[1]] = ConstantArray[1., Length[x]];
T[[2]] = x;
Tprime[[1]] = ConstantArray[0., Length[x]];
Tprime[[2]] = ConstantArray[1., Length[x]];
Do[T[[i]] = 2. x T[[i - 1]] - T[[i - 2]];
Tprime[[i]] =
2. T[[i - 1]] + 2. x Tprime[[i - 1]] - Tprime[[i - 2]];, {i, 3,
n + 1}];
{Transpose[T], Transpose[Tprime]}];

n = 9;(*Collocation points*)
B = 0.11; ita = 1.0; F = 0.28; L = 10;(*parameters*)

(*Definig Collocation points*)
y = -Cos[N[Pi/n] Range[0., n]];
{T, Tprime} = recurrenceTables[y, n];
Tprime = 2 Tprime/L;
(*Defining all the functions we need*)
Z[x_] := L (x + 1)/2;
v[Z_] := Exp[(F/4.0 ita) Z];
a[Z_] := Exp[-(F/4.0 ita) Z];
h[Z_] :=
Exp[-(F/4.0 ita) Z] Sqrt[((4.0 B)/
F) + (1 - ((4.0 B)/F)) Exp[(F/4.0 ita) Z]];

Zy = Map[Z, y];
Vbar = DiagonalMatrix[SparseArray[Map[v, Zy]]];
Vbarprime = DiagonalMatrix[SparseArray[Map[v', Zy]]];

Abar = DiagonalMatrix[SparseArray[Map[a, Zy]]];
Abarprime = DiagonalMatrix[SparseArray[Map[a', Zy]]];

Hbar = DiagonalMatrix[SparseArray[Map[h, Zy]]];
Hbarprime = DiagonalMatrix[SparseArray[Map[h', Zy]]];

(*creating entries of the final U and W matrix below, \
solution={a1,h1,v1}*)
eq1 = \[Lambda] a1[z] + a[z] v1'[z] + a'[z] v1[z] + v[z] a1'[z] +
v'[z] a1[z];
eq2 = \[Lambda] a[z] h1[z] - a[z] h[z] v1'[z] - a[z] v'[z] h1[z] -
h[z] v'[z] a1[z] + v[z] h[z] a1'[z] + v[z] a'[z] h1[z] +
h[z] a'[z] v1[z] + (F/(2.0*ita)) h1[z];
eq3 = \[Lambda] (h[z]^2 a1[z] - 2.0 a[z] h[z] h1[z]) -
2.0 a[z] v[z] h[z] h1'[z] - 2.0 a[z] v[z] h'[z] h1[z] -
2.0 a[z] h[z] h'[z] v1[z] - 2.0 v[z] h[z] h'[z] a1[z] +
v[z] h[z]^2 a1'[z] + 2.0 v[z] a'[z] h[z] h1[z] +
h[z]^2 a'[z] v1[z] - ((2.0*B)/ita) a[z] a1[z];

T2 = Abar . T; T31 = Hbar . Hbar . T;
T32 = -2.0 Abar . Hbar .
T; U = -ArrayFlatten[{{T, 0. T, 0. T}, {0. T, T2, 0. T}, {T31, T32,
0. T}}];
A1 = Vbar . Tprime + Vbarprime . T; H1 = 0. T;
V1 = Abar . Tprime + Abarprime . T;
V2 = Hbar . Abarprime . T - Abar . Hbar . Tprime;
A2 = Vbar . Hbar . Tprime - Hbar . Vbarprime . T;
H2 = Vbar . Abarprime . T - Abar . Vbarprime . T + (F/(2.0*ita)) T;
H3 = -2. Abar . Vbar . Hbar . Tprime -
2. Abar . Vbar . Hbarprime . T + 2. Vbar . Abarprime . Hbar . T;
V3 = Hbar . Hbar . Abarprime . T - 2.0 Abar . Hbar . Hbarprime . T;
A3 = Vbar . Hbar . Hbar . Tprime -
2.0 Vbar . Hbar . Hbarprime . T - (2.0 B/ita) Abar . T; W =
ArrayFlatten[{{A1, H1, V1}, {A2, H2, V2}, {A3, H3, V3}}];


With this code we compute eigenvalues as follows

(*final generalized eigenvalue problem \[Lambda]UX=WX*)Vals =
Union[Eigenvalues[{W, U}]]

(*Out[]= {-2.24273 + 0. I, -1.68633 + 0. I, -1.28402 +
0. I, -0.970793 - 0.514233 I, -0.970793 + 0.514233 I, -0.888949 +
0. I, -0.665724 + 0. I, -0.449759 + 0. I, -0.269871 -
0.459168 I, -0.269871 + 0.459168 I, -0.245564 + 0. I, -0.119623 +
0. I, -0.0396647 + 0. I, -0.00916759 + 0. I, -0.00134319 + 0. I,
1.94998*10^-7 + 0. I, 0.0615641 - 0.304778 I, 0.0615641 + 0.304778 I,
0.203549 - 0.105815 I, 0.203549 + 0.105815 I, ComplexInfinity}*)


To be sure of the results we get, we have compiled second code

ClearAll["Global*"];
n = 9;(*Collocation points*)
B = 0.11; ita = 1.0; F = 0.28; L = 10;(*parameters*)

(*Definig Collocation points*)
y = -Cos[N[Pi/n] Range[0., n]];

vec[x_] := Table[ChebyshevT[n, x], {n, 0, 9}]; Map[vec, y];
vec1[x_] := Table[i ChebyshevU[-1 + i, x], {i, 0, 9}]; Map[vec1, y];

(*Defining all the functions we need*)
Z[x_] := L (x + 1)/2;
v[Z_] := Exp[(F/4.0 ita) Z];
a[Z_] := Exp[-(F/4.0 ita) Z];
h[Z_] :=
Exp[-(F/4.0 ita) Z] Sqrt[((4.0 B)/
F) + (1 - ((4.0 B)/F)) Exp[(F/4.0 ita) Z]];

Zy = Map[Z, y];
V0 = DiagonalMatrix[SparseArray[Map[v, Zy]]];
V0p = DiagonalMatrix[SparseArray[Map[v', Zy]]];

A0 = DiagonalMatrix[SparseArray[Map[a, Zy]]];
A0p = DiagonalMatrix[SparseArray[Map[a', Zy]]];

H0 = DiagonalMatrix[SparseArray[Map[h, Zy]]];
H0p = DiagonalMatrix[SparseArray[Map[h', Zy]]];

AA1 = Array[aa, n + 1]; HH1 = Array[hh, n + 1]; VV1 =
Array[vv, n + 1]; a1[y_] := AA1 . vec[y]; h1[y_] := HH1 . vec[y];
v1[y_] := VV1 . vec[y]; a1p[y_] := AA1 . vec1[y];
h1p[y_] := HH1 . vec1[y]; v1p[y_] := VV1 . vec1[y];
A1 = Map[a1, y]; H1 = Map[h1, y]; V1 = Map[v1, y]; A1p =
2/L Map[a1p, y]; H1p = 2/L Map[h1p, y]; V1p = 2/L Map[v1p, y];

eq1 = \[Lambda] a1[z] + a[z] v1'[z] + a'[z] v1[z] + v[z] a1'[z] +
v'[z] a1[z];
eq2 = \[Lambda] a[z] h1[z] - a[z] h[z] v1'[z] - a[z] v'[z] h1[z] -
h[z] v'[z] a1[z] + v[z] h[z] a1'[z] + v[z] a'[z] h1[z] +
h[z] a'[z] v1[z] + (F/(2.0*ita)) h1[z];
eq3 = \[Lambda] (h[z]^2 a1[z] - 2.0 a[z] h[z] h1[z]) -
2.0 a[z] v[z] h[z] h1'[z] - 2.0 a[z] v[z] h'[z] h1[z] -
2.0 a[z] h[z] h'[z] v1[z] - 2.0 v[z] h[z] h'[z] a1[z] +
v[z] h[z]^2 a1'[z] + 2.0 v[z] a'[z] h[z] h1[z] +
h[z]^2 a'[z] v1[z] - ((2.0*B)/ita) a[z] a1[z];

eqs = Join[\[Lambda] A1 + A0 . V1p + A0p . V1 + V0 . A1p +
V0p . A1, \[Lambda] A0 . H1 - A0 . H0 . V1p - A0 . V0p . H1 -
H0 . V0p . A1 + V0 . H0 . A1p + V0 . A0p . H1 +
H0 . A0p .
V1 + (F/(2.0*ita)) H1 , \[Lambda] (H0 . H0 . A1 -
2.0 A0 . H0 . H1 ) - 2.0 A0 . V0 . H0 . H1p -
2.0 A0 . V0 . H0p . H1 - 2.0 A0 . H0 . H0p . V1 -
2.0 V0 . H0 . H0p . A1 + V0 . H0 . H0 . A1p +
2.0 V0 . A0p . H0 . H1 +
H0 . H0 . A0p . V1 - ((2.0*B)/ita) A0 . A1 ];

(*creating entries of the final U and W matrix below*)var =
Join[AA1, HH1, VV1]; {v, mat} = CoefficientArrays[eqs, var]

W = Normal[mat] /. \[Lambda] -> 0; U =
W - Normal[mat] /. \[Lambda] -> 1;


Solution

{val, fun} = Eigensystem[{W, U}];

val // Union

(*Out[]= {-2.24273 + 0. I, -1.68633 + 0. I, -1.28402 +
0. I, -0.970793 + 0.514233 I, -0.970793 - 0.514233 I, -0.888949 +
0. I, -0.665724 + 0. I, -0.449759 + 0. I, -0.269871 -
0.459168 I, -0.269871 + 0.459168 I, -0.245564 + 0. I, -0.119624 +
0. I, -0.0396647 + 0. I, -0.00916761 + 0. I, -0.00134317 + 0. I,
1.94951*10^-7 + 0. I, 0.0615641 - 0.304778 I, 0.0615641 + 0.304778 I,
0.203549 - 0.105815 I, 0.203549 + 0.105815 I, ComplexInfinity}*)

• Many many thanks @Alex Trounev Commented Aug 29, 2023 at 3:39
• I also am looking for eigenvalue(dominant eigenvalue) which will remain almost unchanged with the increase of collocation point n. How to get that? @Alex Trounev Commented Aug 29, 2023 at 3:46
• @Fun123 In your case it should be something around zero. But you need to formulate criteria, for example, $Min|\lambda|$, or $Max Re(\lambda), Re(\lambda)<0$ and so on. Commented Aug 29, 2023 at 7:06
• How to code the criteria in Mathematica? I am sorry, I think I did not quite get it.@Alex Trounev Commented Aug 30, 2023 at 2:26
• @Fun123 It depends on what you want to explore. I have code to compute eigenvalue with minimal Abs and Abs[Re[]]` part. Commented Aug 30, 2023 at 2:49