Nice Integer Eigenvectors - Mmmm - 03-26-2015 02:38 PM
This program gives you nice integer eigenvectors for 2x2 and 3x3 matrices rather than the normalised ones given by the EIGENVV program (or the seemingly random mess that the HP39gii gives out).
EIGEN() will output a list of a matrix of column eigenvectors and its diagonalisation if possible or a vector of eigenvalues if not possible.
EIGVALS() outputs a vector of eigenvalues
EIGVECT() outputs a single eigenvector matrix (vectors as columns).
it only works for matrices with real integer eigenvalues. If it can't solve the problem it falls back to EIGENVV (ie if it is non-integer, complex or not 2x2 or 3x3).
eg: Comparison of output: EIGENVV() Vs EIGEN()
[attachment=1833]
It will also solve some problems that EIGENVV fails on.
eg: EIGENVV() and EIGENVAL() Vs EIGEN() and EIGVALS()
[attachment=1832]
Code:
EIGVALS(matrix); //declare_EIGVALS_here_to_preserve
//function_order_in_user_menu
SIMPLIFY(vector)
BEGIN
LOCAL p,size,elements,factor;
size:=SIZE(vector);
elements:=size(1);
//find common factor of vector elements
factor:=vector(1);
FOR p FROM 2 TO elements DO
factor:=igcd(vector(p),factor);
END;
//simplify
IF factor≠0 THEN
vector:=vector/factor;
END;
//make first element non-zero
FOR p FROM 1 TO elements DO
IF vector(p)≠0 THEN
IF vector(p)≠ABS(vector(p)) THEN
vector:=-1*vector;
END;
BREAK;
END;
END;
RETURN(vector);
END;
EXPORT EIGEN(matrix)
BEGIN
LOCAL m,vals,n,vect,eigmat,row,dim;
LOCAL col,nonparallel,valsmat,size;
//find dimension of matrix
size:=SIZE(matrix);
dim:=size(1);
IF size(2)≠dim THEN
dim:=0;
END;
//initialise vars
eigmat:=IDENMAT(dim);
vect:=MAKEMAT(1,dim);
row:=0;
//populate eigenvalue vector
vals:=EIGVALS(matrix);
//check matrix is compatible with EIGEN
//if not run EIGENVV instead
IF dim<2 OR dim>3 THEN //use EIGENVV instead
//leave row==0
ELSE
FOR row FROM 1 TO dim DO
IF IM(vals(row))≠0 THEN
BREAK;
END;
IF FP(vals(row))≠0 THEN
BREAK;
END;
END;
END;
IF row<(dim+1) THEN
RETURN(EIGENVV(matrix));
END;
//main loop - find each eigenvector
FOR n FROM 1 TO dim DO
//find characteristic matrix
m:=(matrix-vals(n)*IDENMAT(dim));
//Simplify rows of m
FOR row FROM 1 TO dim DO
m(row):=SIMPLIFY(m(row));
END;
IF RANK(m)==0 THEN //matrix == kI so eigenvectors=I
eigmat:=IDENMAT(dim);
BREAK;
END;
//find first non-zero element
FOR row FROM 1 TO dim DO
FOR col FROM 1 TO dim DO
IF m(row,col)≠0 THEN
BREAK;
END;
END;
IF col<(dim+1) THEN
BREAK;
END;
END;
IF dim==2 THEN //2x2 eigenvector calculation
vect(1):=-m(row,2);
vect(2):=m(row,1);
ELSE //3x3 eigenvector calculations
IF RANK(m)==1 THEN //only 1 li row - find a vector perp to m
nonparallel:=m(row);
nonparallel(1+(col MOD dim)):=1+m(row,(1+(col MOD dim)));
vect:=CROSS(m(row),nonparallel);
ELSE //2 li rows - use cross product
IF m(row)≠m(row+1) AND m(row+1)≠[0,0,0] THEN
vect:=CROSS(m(row),m(row+1));
ELSE
vect:=CROSS(m(row),m(row+2));
END;
END;
END;
vect:=SIMPLIFY(vect);
// in case of repeated eigenvalue find vector perp
// to charmatrix and repeated eigenvector,
// but only if we have 3 non li row vects in m and dim>2
IF n>2 AND dim>2 THEN
IF vect==eigmat(n-2) AND RANK(m)==1 THEN //find a perp vector
vect:=CROSS(vect,m(row));
vect:=SIMPLIFY(vect);
END;
END;
IF n>1 AND dim>2 THEN
IF vect==eigmat(n-1) AND RANK(m)==1 THEN
vect:=CROSS(vect,m(row));
vect:=SIMPLIFY(vect);
END;
END;
//populate eigenvector matrix (rows)
eigmat(n):=vect;
END;
//zero repeated eigenvectors
FOR n FROM 1 TO dim DO
IF eigmat(n)==eigmat(1+(n MOD dim)) THEN
FOR col FROM 1 TO dim DO
eigmat(MAX(n,(1+(n MOD dim))),col):=0;
END;
END;
END;
//change to a column matrix
eigmat:=TRN(eigmat);
//display result
IF DET(eigmat)≠0 THEN
valsmat:=eigmat*matrix*eigmat;
//diagonalise and show as a list
CONCAT(eigmat,valsmat);
ELSE
CONCAT(eigmat,vals);
END;
END;
EXPORT EIGVECT(matrix)
BEGIN
//brief function to return the vector matrix only
LOCAL result;
result:=EIGEN(matrix);
result(1);
END;
EXPORT EIGVALS(matrix)
BEGIN
LOCAL a,b,c,d,eqn,size,dim;
//determine the dimensions of the matrix
size:=SIZE(matrix);
dim:=size(1);
IF size(2)≠dim THEN
dim:=0;
END;
//use built in eigenvalue command if matrix not compatible
IF dim<2 OR dim>3 THEN
RETURN(EIGENVAL(matrix));
END;
IF dim==2 THEN
//Generate coefficients for characteristic equation (quadratic eqn)
a:=1;
b:=-(matrix(1,1)+matrix(2,2));
c:=DET(matrix);
//form quadratic coefficient vector
eqn:=[0,0,0];
eqn(1):=ROUND(a,8);
eqn(2):=ROUND(b,8);
eqn(3):=ROUND(c,8);
ELSE //(if dim==3)
//Generate coefficients for characteristic equation (cubic eqn)
a:=-1;
b:=TRACE(matrix);
c:=-matrix(1,1)*matrix(3,3)-matrix(1,1)*matrix(2,2)-matrix(2,2)*matrix(3,3)+matrix(1,3)*matrix(3,1)+matrix(1,2)*matrix(2,1)+matrix(2,3)*matrix(3,2);
d:=matrix(1,1)*matrix(2,2)*matrix(3,3)+matrix(1,2)*matrix(2,3)*matrix(3,1)+matrix(1,3)*matrix(2,1)*matrix(3,2)-matrix(1,1)*matrix(2,3)*matrix(3,2)-matrix(3,3)*matrix(1,2)*matrix(2,1)-matrix(2,2)*matrix(1,3)*matrix(3,1);
//form cubic coefficient vector
eqn:=[0,0,0,0];
eqn(1):=ROUND(a,8);
eqn(2):=ROUND(b,8);
eqn(3):=ROUND(c,8);
eqn(4):=ROUND(d,8);
END;
//solve polynomial
ROUND(POLYROOT(eqn),8);
END;
RE: Nice Integer Eigenvectors - Han - 03-31-2015 12:55 PM
The command eigenvects() -- in lower case -- gives the eigenvectors in rational form and works for larger sized matrices. You could probably write a simpler function that then converts rational vectors to a parallel integer vector (and apply gcd for smaller integers).
RE: Nice Integer Eigenvectors - Mmmm - 03-31-2015 09:59 PM
(03-31-2015 12:55 PM)Han Wrote: The command eigenvects() -- in lower case -- gives the eigenvectors in rational form and works for larger sized matrices. You could probably write a simpler function that then converts rational vectors to a parallel integer vector (and apply gcd for smaller integers).
WOW! That is a great improvement over the ole EIGENVV, why didn't anyone tell me about this command?
But unfortunately it falls short - it seems to ignore repeated eigenvalues.
eg.
running it on
5|-4|0
1|0|2
0|2|5
gives only one eigenvector (-4,-5,2) which corresponds to the eigenvalue 0. There are also two more eigenvalues: 5 and 5 for which there is one more eigenvector.
My program gives (2,0,-1) as well as (4,5,-2).
I've tried it on some other matrices with repeated eigenvalues too and it always ignores them!!
Shame...
|