1
$\begingroup$

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}]
$\endgroup$
12
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Aug 23, 2023 at 5:43
  • $\begingroup$ 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/… $\endgroup$
    – Fun123
    Commented Aug 23, 2023 at 6:30
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Aug 23, 2023 at 6:34
  • $\begingroup$ Okay, I see. Thanks for the clarification @Nasser $\endgroup$
    – Fun123
    Commented Aug 23, 2023 at 6:40
  • $\begingroup$ @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/… $\endgroup$ Commented Aug 23, 2023 at 10:11

1 Answer 1

1
$\begingroup$

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}*)
$\endgroup$
16
  • $\begingroup$ Many many thanks @Alex Trounev $\endgroup$
    – Fun123
    Commented Aug 29, 2023 at 3:39
  • $\begingroup$ 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 $\endgroup$
    – Fun123
    Commented Aug 29, 2023 at 3:46
  • $\begingroup$ @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. $\endgroup$ Commented Aug 29, 2023 at 7:06
  • $\begingroup$ How to code the criteria in Mathematica? I am sorry, I think I did not quite get it.@Alex Trounev $\endgroup$
    – Fun123
    Commented Aug 30, 2023 at 2:26
  • $\begingroup$ @Fun123 It depends on what you want to explore. I have code to compute eigenvalue with minimal Abs and Abs[Re[]] part. $\endgroup$ Commented Aug 30, 2023 at 2:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.