This is an update to make it handle geometric multiplicity>1.
When an eigenvalue has geometric multiplicity=1, then this will use the basic method of solving $(A-\lambda I) v_2 = v_1$ to find the generalized eigenvector $v_2$ given $v_1$, and repeats this as many times as needed per each defective $\lambda$.
To find the first generalized $v_2$, in the above $v_1$ will be the normal eigenvector, which is found by Eigensystem. If more generalized eigenvectors are needed (because the eigenvalue could have large defect), then the next step will use the now found $v_2$ in order to find $v_3$ and so on. This is done for each defective $\lambda$.
When an eigenvalue has geometric multiplicity>1 however, the generalized eigenvector is found from the null space of the transpose of the matrix made up of the current set of known eigenvectors. We just pick one random row from the nullspace and use it (after transposing) as new generalized eigenvector and add it to the set of known eigenvectors for that $\lambda$. If we need more generalized eigenvectors, we repeat this process. The function NullSpace is used in Mathematica for this.
Example 1
mat = {{4, 1, 1}, {-2, 1, -2}, {1, 1, 4}};
findV[mat, V];
Gives
Example 2
mat = {{3, 1}, {0, 3}}
findV[mat, V];
Example 3
mat = {{2, 1}, {-1, 4}};
findV[mat, V];
Example 4
From https://en.wikipedia.org/wiki/Generalized_eigenvector
mat = {{1, 0, 0, 0, 0}, {3, 1, 0, 0, 0}, {6, 3, 2, 0, 0}, {10, 6, 3,
2, 0}, {15, 10, 6, 3, 2}};
findV[mat, V];
Example 5
mat = {{5, 1, -4}, {4, 3, -5}, {3, 1, -2}};
findV[mat, V];
Example 6
mat = {{5, 1, 3, 2}, {0, 5, 0, -3}, {0, 0, 5, 1}, {0, 0, 0, 5}};
findV[mat, V];
Example 7
mat = {{39, 8, -16}, {-36, -5, 16}, {72, 16, -29}};
findV[mat, V];
Example 8
mat = {{1, 0, 0}, {1, 3, 1}, {-2, -4, -1}};
findV[mat, V];
Example 9
mat = {{3, 1, 0, 0}, {0, 3, 1, 0}, {0, 0, 3, 0}, {0, 0, 0, 3}};
findV[mat, V];
Note about Jordan form
Jordan form is unique up to permutation of its Jordan blocks.
When there is one block, like in example 9 above, since there is only one $\lambda$, the arrangement of columns of $P$ matrix (the matrix whose columns are the eigenvectors of this $\lambda$) is arbitrary. Only condition is that they are all linearly independent of each others. So you might not get same Jordan form as Mathematica gives if the order of eigenvectors in $P$ happened to be different than what Mathematica used. Here is for the above.
(A = {{3, 1, 0, 0}, {0, 3, 1, 0}, {0, 0, 3, 0}, {0, 0, 0, 3}}) // MatrixForm
JordanDecomposition[A][[2]] // MatrixForm
Here is the $P$ matrix of the eigenvectors found for $\lambda$
(p = {{0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {1, 0, 0, 0}}) // MatrixForm
Now
(Inverse[p] . A . p) // MatrixForm
This does not have same arrangement what Mathematica gave.
But if we simply swap the 3rd column with the 4th, now it will give same Joprdan matrix as Mathematica
(p = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}}) // MatrixForm
(Inverse[p] . A . p) // MatrixForm
Code
Bug reports are welcome. In case code gets messed up when pasting here as I use code cell, here is link to the notebook also. Used V 13.2.1
(*Change history *)
(*may 10, 2023. Initial Version *)
(*may 10, 2023. Added check for geometric multiplicity >1 for now until supported*)
(*may 11, 2023. Added support for geometric multiplicity >1, improved formatting*)
findV[mat_,V_]:=Module[{data,keys},
{data,keys}=generateInitialTable[mat];
processTable[mat,data,keys,V]
]
(*****************************)
printTable[data_]:=Module[{header,new},
header={"\[Lambda]","algebraic multiplicity","geometric multiplicity","normal eigenvectors","Defective \[Lambda]?"};
Print@Text[Style[
"We first use Mathematica to find the eigenvalues and the norma\n"<>
"eigenvectors and examine the result to determine which eigenvalues\n"<>
"are defective and which are not. This results in this table ",TextAlignment->Left]];
Print[Grid[Join[{header},
MapAt[MatrixForm[Transpose[#]]&,data,{All,4}]
],
Frame->All]]
]
(******************************)
generateInitialTable[mat_]:=Module[{lambda,v,p,keys,k,n,data,idx,theV,header,x},
{lambda,v}=Eigensystem[mat];
Print["The Matrix is ",MatrixForm@mat];
Print["its eigenvalues are ",MatrixForm@lambda];
Print["its normal Eigenvectors are ",MatrixForm@Transpose[v]];
p=PositionIndex[lambda];
keys=Keys[p];
(*this table will hold summary of Eigensystem finding. One row per each lambda*)
(*1=lambda,2=algebraic multiplicity,3=geometric multiplicity,4=normal V,5=defective?*)
data=Table[{0,0,0,0,0},Length@keys];
(*Now build the table *)
Do[
k=keys[[n]];
data[[n,1]]=k;
data[[n,2]]=Length@p[k];(*algebraic multiplicity*)
idx=p[k];
theV=v[[idx[[1]];;idx[[-1]]]];
theV=Cases[theV,x_/;Not@AllTrue[x,#==0&]];
data[[n,3]]=Length@theV;(*geometric multiplicity*)
data[[n,4]]=theV;
data[[n,5]]=If[data[[n,2]]>data[[n,3]],True,False]
,
{n,Length@keys}
];
printTable[data];
{data,keys}
]
(******************************)
(*geometric multiplicity =1, easy case *)
handleCaseOne[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,setOfV},
howMuchDefective=data[[n,2]]-data[[n,3]];
Print["Since eigenvalue ",data[[n,1]]," is defective, and the degree of defect is ",
howMuchDefective," then we need to find ",howMuchDefective,
" additional generalized eigenvectors for this eigenvalue."];
currentV=data[[n,4]]; (*normal V found*)
setOfV=currentV;
Do[
g=Table[V[i],{i,1,size}];
eq=mat.g==currentV;
Print["Now we solve for generalized eigevector ",MatrixForm[g],
" Using (A-\[Lambda]I).g=v1. This gives"];
Print[MatrixForm[mat],"-(",data[[n,1]],")",
MatrixForm[IdentityMatrix[size]],".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];
Print["The above equation becomes\n",MatrixForm[mat-data[[n,1]]*IdentityMatrix[size]],
".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];
sol=LinearSolve[(mat-data[[n,1]]*IdentityMatrix[size]),currentV[[1]]];
Print["Using LinearSolve gives ",MatrixForm[g],"=",Style[MatrixForm@sol,Bold]];
currentV={sol};
setOfV=Join[setOfV, currentV]
,
{m,1,howMuchDefective}
];
setOfV
]
(******************************)
(* when geometric multiplicity >1, hard case *)
handleCaseTwo[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,newV,i},
Print["This case of geometric multiplicity >1."];
howMuchDefective=data[[n,2]]-data[[n,3]];
Print["Since eigenvalue ",data[[n,1]]," is defective, and the\n",
"degree of defect is ",howMuchDefective," then we need to find ",howMuchDefective,"\n",
"additional generalized eigenvectors for this eigenvalue."];
currentV=data[[n,4]]; (*normal V's found. More than one*)
g=Table[V[i],{i,1,size}];
Do[
Print["Now we solve need to find a generlized eigenvector ",MatrixForm@g,"\n",
"which is linearly independent to the set of eigenvectors already found ",
MatrixForm[Transpose[currentV]],"\n",
"Which we do by looking at the NullSpace of the transpose of the above as a matrix\n",
"Which gives ",MatrixForm[NullSpace[currentV]],"\n",
"And we need to just pick one of these rows as our new eigenvector.\n"
];
newV=NullSpace[currentV];
newV=newV[[1,All]];
currentV=Join[currentV,{newV}];
Print["Taking the first row, and transposing it gives the generalized eigenvector as ",MatrixForm@g," = ",
Style[MatrixForm[Transpose[{newV}]],Bold],"\nwhich we now add to the set of known eigenvectors found giving\n",
MatrixForm[Transpose@currentV]
];
,
{m,1,howMuchDefective}
];
currentV
]
(******************************)
processTable[mat_,data_,keys_,V_]:=Module[{n,setOfV},
Do[
If[data[[n,5]]==True,(*i.e. defective eigenvalue*)
(*check if geometric multiplicity=1 or not. These cases are handled differently*)
If[data[[n,3]]==1,(*standard case, where one normal vector is found*)
(*geometric multiplicity =1 case*)
setOfV=handleCaseOne[data,mat,V,n];
Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
Style[MatrixForm[Transpose@setOfV],Bold]
]
,
(*geometric multiplicity >1 case*)
setOfV=handleCaseTwo[data,mat,V,n];
Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
Style[MatrixForm[Transpose@setOfV],Bold]
]
]
,
Print["Since eigenvalue ",data[[n,1]]," is not defective, then we already have its normal "];
Print["eigenvectors from the table above, which was found by Eigensystem. Nothing to do."]
]
,
{n,Length@keys}
]
]
Eigensystemdoes not return the generalized eigenvectors. $\endgroup$