load("eigen")$

MATRIX1:MATRIX([M1,0,0,0,M5],[0,M2,0,0,M5],[0,0,M3,0,M5],[0,0,0,M4,M5],[0,0,0,0,0]);

MATRIX2:MATRIX([1,2,3,4],[0,3,4,5],[0,0,5,6],[0,0,0,9]);

/* THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE
	EIGENVALUES AND THE UNIT EIGEN VECTORS OF THAT MATRIX... */

UNITEIGENVECTORS(MATRIX2);

/* IF YOU ALREADY KNOW THE EIGENVECTORS YOU CAN SET THE FLAG 
	KNOWNEIGVECTS TO TRUE AND THE GLOBAL VARIABLE LISTEIGVECTS TO THE
	LIST OF THE EIGEN VECTORS... 
	THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE EIGEN
	VALUES AND THE UNIT EIGEN VECTORS OF THAT MATRIX.  IN ADDITION IF
	THE FLAG NONDIAGONALIZABLE IS FALSE,TWO GLOBAL MATRICES LEFTMATRIX AND
	RIGHTMATRIX WILL BE GENERATED.  THESE MATRICES HAVE THE PROPERTY THAT
	LEFTMATRIX.(MATRIX).RIGHTMATRIX IS A DIAGONAL MATRIX WITH THE EIGEN 
	VALUES OF THE (MATRIX) ON THE DIAGONAL...  */

SIMILARITYTRANSFORM(MATRIX1)$
NONDIAGONALIZABLE;
RATSIMP(LEFTMATRIX.MATRIX1.RIGHTMATRIX);

/* NOW THAT YOU KNOW HOW TO USE THE EIGEN PACKAGE, HERE ARE SOME
	EXAMPLES ABOUT HOW NOT TO USE IT.
	CONSIDER THE FOLLOWING MATRIX :   */

MATRIX3:MATRIX([1,0],[0,1]);

/* AS YOU'VE UNDOUBTEDLY NOTICED, THIS IS THE 2*2 IDENTITY MATRIX.
	LET'S FIND THE EIGEN VALUES AND THE EIGEN VECTORS OF THIS MATRIX...
*/

EIGENVECTORS(MATRIX3);

/* "NOTHING SPECIAL HAPPENED", YOU SAY.  EVERYONE KNOWS WHAT THE EIGEN
	VALUES AND THE EIGEN VECTORS OF THE IDENTITY MATRIX ARE, RIGHT?
	RIGHT.  NOW CONSIDER THE FOLLOWING MATRIX :  */

MATRIX4:MATRIX([1,E],[E,1]);

/* LET E>0, BUT AS SMALL AS YOU CAN IMAGINE.  SAY 10^(-100).
	LET'S FIND THE EIGEN VALUES AND THE EIGEN VECTORS OF THIS MATRIX :
*/

EIGENVECTORS(MATRIX4);

/* SINCE E~10^(-100), THE EIGEN VALUES OF MATRIX4 ARE EQUAL TO THE
	EIGEN VALUES OF MATRIX3 TO A VERY GOOD ACCURACY.  BUT, LOOK
	AT THE EIGEN VECTORS!!!  EIGEN VECTORS OF MATRIX4 ARE NOWHERE
	NEAR THE EIGEN VECTORS OF MATRIX3.  THERE IS ANGLE OF %PI/4
	BETWEEN THE CORRESPONDING EIGEN VECTORS.  SO, ONE LEARNS
	ANOTHER FACT OF LIFE :

	MATRICES WHICH HAVE APPROXIMATELY THE SAME EIGEN VALUES DO NOT
	HAVE APPROXIMATELY THE SAME EIGEN VECTORS IN GENERAL.

	THIS EXAMPLE MAY SEEM ARTIFICIAL TO YOU, BUT IT IS NOT IF YOU THINK
	A LITTLE BIT MORE ABOUT IT.  SO, PLEASE BE CAREFUL WHEN YOU
	APPROXIMATE THE ENTRIES OF WHATEVER MATRIX YOU HAVE.  YOU MAY
	GET GOOD APPROXIMATIONS TO ITS EIGEN VALUES, HOWEVER THE EIGEN
	VECTORS YOU GET MAY BE ENTIRELY SPURIOUS( OR SOME MAY BE CORRECT,
	BUT SOME OTHERS MAY BE TOTALLY WRONG ).

	NOW, HERE IS ANOTHER SAD STORY :
	LET'S TAKE A LOOK AT THE FOLLOWING MATRIX :  */

MATRIX5:MATRIX([5/2,50-25*%I],[50+25*%I,2505/2]);

/* NICE LOOKING MATRIX, ISN'T IT?  AS USUAL, WE WILL FIND THE EIGEN
	VALUES AND THE EIGEN VECTORS OF IT :  */

EIGENVECTORS(MATRIX5);

/* WELL, HERE THEY ARE.  SUPPOSE THAT THIS WAS NOT WHAT YOU WANTED.
	INSTEAD OF THOSE SQRT(70)'S, YOU WANT THE NUMERICAL VALUES OF
	EVERYTHING.  ONE WAY OF DOING THIS IS TO SET THE FLAG "NUMER"
	TO TRUE AND USE THE EIGENVECTORS COMMAND AGAIN :  */

NUMER:TRUE;
EIGENVECTORS(MATRIX5);

/* OOOPS!!!  WHAT HAPPENED??  WE GOT THE EIGEN VALUES, BUT THERE ARE
	NO EIGENVECTORS.  NONSENSE, THERE MUST BE A BUG IN EIGEN, RIGHT?
	WRONG.  THERE IS NO BUG IN EIGEN.  WE HAVE DONE SOMETHING WHICH
	WE SHOULD NOT HAVE DONE.  LET ME EXPLAIN : 
	WHEN ONE IS SOLVING FOR THE EIGEN VECTORS, ONE HAS TO FIND THE
	SOLUTION TO HOMOGENEOUS EQUATIONS LIKE :  */

EQUATION1:A*X+B*Y=0;
EQUATION2:C*X+D*Y=0;

/* IN ORDER FOR THIS SET OF EQUATIONS TO HAVE A SOLUTION OTHER THAN
	THE TRIVIAL SOLUTION ( THE ONE IN WHICH X=0 AND Y=0 ), THE 
	DETERMINANT OF THE COEFFICIENTS ( IN THIS CASE A*D-B*C ) SHOULD
	VANISH.  EXACTLY.   IF THE DETERMINANT DOES NOT VANISH THE ONLY
	SOLUTION WILL BE THE TRIVIAL SOLUTION AND WE WILL GET NO EIGEN
	VECTORS.  DURING THIS DEMO, I DID NOT SET A,B,C,D TO ANY
	PARTICULAR VALUES.  LET'S SEE WHAT HAPPENS WHEN WE TRY TO SOLVE
	THE SET ABOVE :  */

ALGSYS([EQUATION1,EQUATION2],[X,Y]);

/* YOU SEE?  THE INFAMOUS TRIVIAL SOLUTION.  NOW LET ME SET A,B,C,D
	TO SOME NUMERICAL VALUES :  */

A:4;
B:6;
C:2;
D:3;
A*D-B*C;
EQUATION1:EV(EQUATION1);
EQUATION2:EV(EQUATION2);
ALGSYS([EQUATION1,EQUATION2],[X,Y]);

/* NOW WE HAVE A NONTRIVIAL SOLUTION WITH ONE ARBITRARY CONSTANT.
	( %R(SOMETHING) ).  WHAT HAPPENED IN THE PREVIOUS CASE IS THAT
	THE NUMERICAL ERRORS CAUSED THE DETERMINANT NOT TO VANISH, HENCE
	ALGSYS GAVE THE TRIVIAL SOLUTION AND WE GOT NO EIGEN VECTORS.
	IF YOU WANT A NUMERICAL ANSWER, FIRST CALCULATE IT EXACTLY,
	THEN SET "NUMER" TO TRUE AND EVALUATE THE ANSWER.  */

NUMER:FALSE;
NOTNUMERICAL:EIGENVECTORS(MATRIX5);
NUMER:TRUE;
EV(NOTNUMERICAL);

/* YOU SEE, IT WORKS NOW.  ACTUALLY, IF YOU HAVE A MATRIX WITH
	NUMERICAL ENTRIES AND YOU CAN LIVE WITH REASONABLY ACCURATE
	ANSWERS, THERE ARE MUCH BETTER (FASTER) PROGRAMS.  ASK SOMEBODY
	ABOUT THE IMSL ROUTINES ON THE SHARE DIRECTORY... 
	THIS IS ALL...  IF YOU THINK THAT THE NAMES OF THE FUNCTIONS ARE TOO 
	LONG, THERE ARE SHORTER NAMES FOR THEM AND THEY ARE GIVEN IN THE FILE
	EIGEN USAGE DSK:SHARE;.  GOOD LUCK!!!!!!!!!!!!!...... YEKTA  */


"end of demo"$