/* modified for DOE MACSYMA */ /* THIS IS THE FILE EIGEN DEMO DSK:SHARE;. (YOU CAN BATCH OR DEMO THIS FILE, I.E. BATCH(EIGEN,DEMO,DSK,SHARE);, OR DEMO(EIGEN,DEMO,DSK,SHARE);. NOTE THAT IN THE DEMO MODE YOU HAVE TO HIT THE SPACE KEY AFTER EACH STEP...) THE FUNCTIONS IN THE NEW EIGEN PACKAGE ARE DEMONSTRATED HERE. THE DESCRIPTION OF THE FUNCTIONS CAN BE FOUND IN THE FILE EIGEN USAGE DSK:SHARE;, THE SOURCE CODE IS ON THE FILE EIGEN > DSK:SHARE; AND THE FASTLOAD FILE IS EIGEN FASL DSK:SHARE;. ( YOU CAN LOAD THIS ONE USING MACSYMA'S LOADFILE COMMAND, I.E. LOADFILE(EIGEN,FASL,DSK,SHARE);.) WE START WITH LOADING THE EIGEN PACKAGE : */ /* IF NOT STATUS(FEATURE,EIGEN) THEN LOADFILE(EIGEN,FASL,DSK,SHARE); */ load("eigen.mc")$ /* LET US START WITH THE FIRST FUNCTION. (SEE THE DESCRIPTIONS...) FIRST LET'S DEFINE A COMPLEX VARIABLE... */ Z:A+%I*B; /* THE CONJUGATE FUNCTION SIMPLY RETURNS THE COMPLEX CONJUGATE OF ITS ARGUMENT... */ CONJUGATE(Z); /* NOTE THAT Z COULD BE A MATRIX, A LIST, ETC... */ Z:MATRIX([%I,0],[0,1+%I]); CONJUGATE(Z); /* THE NEXT FUNCTION CALCULATES THE INNER PRODUCT OF TWO LISTS...*/ LIST1:[A,B,C,D]; LIST2:[F,G,H,I]; INNERPRODUCT(LIST1,LIST2); /* THE ELEMENTS OF THE LISTS COULD BE COMPLEX ALSO... */ LIST1:[A+%I*B,C+%I*D]; INNERPRODUCT(LIST1,LIST1); /* THE NEXT FUNCTION TAKES A LIST AS ITS ARGUMENT AND RETURNS A UNIT LIST... */ LIST1:[1,1,1,1,1]; UNITVECTOR(LIST1); LIST2:[1,%I,1-%I,1+%I]; UNITVECTOR(LIST2); /* THE NEXT FUNCTION TAKES A LIST AS ITS ARGUMENT AND RETURNS A COLUMN VECTOR... */ LIST1:[A,B,C,D]; COLUMNVECTOR(LIST1); /* THE NEXT FUNCTION TAKES A LIST OF LISTS AS ITS ARGUMENT AND ORTHOGONALIZES THEM USING THE GRAM-SCHMIDT ALGORITHM...*/ LISTOFLISTS:[[1,2,3,4],[0,5,4,7],[4,5,6,7],[0,0,1,0]]; /* NOTE THAT THE LISTS IN THIS LIST OF LISTS ARE NOT ORTHOGONAL TO EACH OTHER... */ INNERPRODUCT([1,2,3,4],[0,5,4,7]); INNERPRODUCT([1,2,3,4],[4,5,6,7]); /* BUT AFTER APPLYING THE GRAMSCHMIDT FUNCTION... */ ORTHOGONALLISTS:GRAMSCHMIDT(LISTOFLISTS); INNERPRODUCT(PART(ORTHOGONALLISTS,1),PART(ORTHOGONALLISTS,2)); INNERPRODUCT(PART(ORTHOGONALLISTS,2),PART(ORTHOGONALLISTS,3)); /* NOTE THAT ORHTOGONALLISTS CONTAINS INTEGERS THAT ARE FACTORED. IF YOU DO NOT LIKE THIS FORM, YOU CAN SIMPLY RATSIMP THE RESULT : */ RATSIMP(ORTHOGONALLISTS); /* THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE EIGENVALUES OF THAT MATRIX... */ 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]); /* THIS IS THE MATRIX THAT CAUSED A LOT OF TROUBLE FOR THE OLD EIGEN PACKAGE... IT TOOK ~170 SECONDS TO FIND THE EIGEN VECTORS OF THIS MATRIX... YOU SHOULD BE ABLE TO DO IT IN YOUR HEAD IN ABOUT 20 SECONDS [note: PDP-10 timings] ... THE NEW EIGEN PACKAGE HANDLES IT IN ABOUT 10 SECONDS... ANYWAY, LET'S KEEP GOING... */ EIGENVALUES(MATRIX1); /* THE FIRST SUBLIST IN THE ANSWER IS THE EIGENVALUES, SECOND LIST IS THEIR MULTIPLICITIES IN THE CORRESPONDING ORDER... THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE EIGEN VALUES AND THE EIGEN VECTORS OF THAT MATRIX... */ EIGENVECTORS(MATRIX1); /* FIRST SUBLIST IN THE ANSWER IS THE OUTPUT OF THE EIGENVALUES COMMAND THE OTHERS ARE THE EIGEN VECTORS CORRESPONDING TO THOSE EIGEN VALUES... NOTICE THAT THIS COMMAND IS MORE POWERFUL THAN THE EIGENVALUES COMMAND BECAUSE IT DETERMINES BOTH THE EIGEN VALUES AND THE EIGEN VECTORS... IF YOU ALREADY KNOW THE EIGEN VALUES, YOU CAN SET THE KNOWNEIGVALS FLAG TO TRUE AND THE GLOBAL VARIABLE LISTEIGVALS TO THE LIST OF EIGEN VALUES... THIS WILL MAKE THE EXECUTION OF EIGENVECTORS COMMAND FASTER BECAUSE IT DOESN'T HAVE TO FIND THE EIGEN VALUES ITSELF... */ /* commented out here and placed in [share]eigen.dm1 because lack of GC in 8.11 version of DOE MACSYMA 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"$