From 1d48b730160e6ac9066e1fe68d6fad0ca2041b46 Mon Sep 17 00:00:00 2001 From: Eric Swenson Date: Fri, 22 Jun 2018 12:48:55 -0700 Subject: [PATCH] Resolves #1012: Add SHARE;EIGEN DEMO and SHARE;EIGEN USAGE. --- doc/share/eigen.usage | 246 +++++++++++++++++++++++++ src/share/eigen.demo | 404 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 650 insertions(+) create mode 100755 doc/share/eigen.usage create mode 100755 src/share/eigen.demo diff --git a/doc/share/eigen.usage b/doc/share/eigen.usage new file mode 100755 index 00000000..af0bdc7a --- /dev/null +++ b/doc/share/eigen.usage @@ -0,0 +1,246 @@ + + This is the file EIGEN USAGE DSK:SHARE; . + + The EIGEN package is described in this file. The source code +is in 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); . If you access the FASL +file via the functions EIGENVALUES or EIGENVECTORS, the FASL file will +autoload.) The DEMO file is 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 should hit the space key at each step... + + The new EIGEN package is written by Yekta G"ursel (YEKTA@MIT-MC) and +is faster and more memory efficient than the old EIGEN package. It also is +able to handle multiple eigenvalues and the eigenvectors corresponding +to those eigenvalues. It will work with any square matrix ( not necessarily +symmetric or hermitian ) and will tell whether the matrix is diagonalizable. +The calculated eigenvectors and the unit eigenvectors of the matrix are the +RIGHT eigenvectors and the RIGHT unit eigenvectors respectively. +( You should be aware of the fact that this program uses the MACSYMA functions +SOLVE and ALGSYS and if SOLVE can not find the roots of the characteristic +polynomial of the matrix or if it generates a rather messy solution the +EIGEN package may not produce any useful results. More info on this will be +given in the description of the commands... This package is DESIGNED to try +to get the EXACT solutions to the eigenvalue and eigenvector problems. If the +matrices you have contain FLOATING point numbers, it may not be able to solve +your problem. You should use the IMSL eigenvalue and eigenvector package for +numerical matrices with floating point numbers. These excellent routines will +find the APPROXIMATE solutions for numerical matrices with floating point +numbers. See MACSYMA Manual Version 10, Volume 2, Page V2-4-78.) + + If MODULUS is not false, i.e. MODULUS:7, then EIGENVALUES will +automatically call, EIGENFINITEFIELD, to give the eigenvalues of your +matrix over the finite field. EIGENFINITEFIELD calls FACTOR. This capability +gives exact answers for matrices over finite field. This was written by +Nicholas C. Strauss (NCS@MIT-MC). + + Another command of interest is JORDANFORM. For algebraically closed +fields, for example the complex numbers, this will always return a matrix +which is the Jordan canonical form of the input matrix. For non-algebraically +closed fields, a matrix will be returned only in the event that all eigenvalues +are in the field. + + Otherwise, a list of lists similar to the EIGENVALUES list will be +returned. The first entry is a list of eigenvalues together with a list which +contains the information on the number and type of Jordan block for each +eigenvalue. the second entry is the statement "Field Not Algebraically Closed." +For example, [[1,[0,1,2]], ... ] states that 1 is an eigenvalue with one +2x2 Jordan block, and two 3x3 Jordan blocks. The list paired with the +eigenvalue lists the Jordan blocks of that eigenvalue in ascending order. +The example [[5,[2]],...], states that the eigenvalue 5 has two 1x1 Jordan +blocks. Thus [[1,[1]],[1,[1]]] denotes the 2x2 Identity matrix. + + The Jordan form for real symmetric matrices is diagonal. +In general, the Jordan form is the similar matrix closest to diagonal form. +This was written by ncs. + + Description of the functions : + +CONJUGATE(X) returns the complex conjugate of its argument. + ( Note that %I's in the expressions should be explicit, since there is + no complex variable declaration in MACSYMA at the present time. This + is true for all the functions in this package...) + + +INNERPRODUCT(X,Y) takes two LISTS of equal length as its arguments and + returns their inner ( scalar ) product defined by + (Complex Conjugate of X).Y ( The "dot" operation is the same + as the usual one defined for vectors. ) . + + +UNITVECTOR(X) takes a LIST as its argument and returns a unit list. + ( i.e. a list with unit magnitude. ) + + +COLUMNVECTOR(X) takes a LIST as its argument and returns a column vector the + components of which are the elements of the list. The first element is + the first component,...etc...( This is useful if you want to use parts + of the outputs of the functions in this package in matrix calculations. + ..) + + +GRAMSCHMIDT(X) takes a LIST of lists the sublists of which are of + equal length and not necessarily orthogonal ( with respect to the + innerproduct defined above ) as its argument and returns a similar + list each sublist of which is orthogonal to all others. + ( Returned results may contain integers that are factored. + This is due to the fact that the MACSYMA function FACTOR is + used to simplify each substage of the Gram-Schmidt algorithm. + This prevents the expressions from getting very messy and + helps to reduce the sizes of the numbers that are produced + along the way. ) + + +EIGENVALUES(MAT) takes a MATRIX as its argument and returns a list of + lists the first sublist of which is the list of eigenvalues of + the matrix and the other sublist of which is the list of the + multiplicities of the eigenvalues in the corresponding order. + { The MACSYMA function SOLVE is used here to find the roots of + the characteristic polynomial of the matrix. Sometimes SOLVE + may not be able to find the roots of the polynomial;in that + case nothing in this package except CONJUGATE, INNERPRODUCT, + UNITVECTOR, COLUMNVECTOR and GRAMSCHMIDT will work unless + you know the eigenvalues. + In some cases SOLVE may generate very messy eigenvalues. You may + want to simplify the answers yourself before you go on. There + are provisions for this and they will be explained below. + ( This usually happens when SOLVE returns a not-so-obviously + real expression for an eigenvalue which is supposed to be real...)} + + +EIGENVECTORS(MAT) takes a MATRIX as its argument and returns a list of + lists the first sublist of which is the output of the EIGENVALUES + command and the other sublists of which are the eigenvectors + of the matrix corresponding to those eigenvalues respectively. + The flags that affect this function are : + + NONDIAGONALIZABLE[FALSE] will be set to TRUE or FALSE depending + on whether the matrix is nondiagonalizable or diagonalizable after + an EIGENVECTORS command is executed. + + HERMITIANMATRIX[FALSE] If set to TRUE will cause the degenerate + eigenvectors of the hermitian matrix to be orthogonalized using + the Gram-Schmidt algorithm. + + KNOWNEIGVALS[FALSE] If set to TRUE the EIGEN package will assume + the eigenvalues of the matrix are known to the user and stored + under the global name LISTEIGVALS. LISTEIGVALS should be set to + a list similar to the output of the EIGENVALUES command. + ( The MACSYMA function ALGSYS is used here to solve for the + eigenvectors. Sometimes if the eigenvalues are messy, ALGSYS may + not be able to produce a solution. In that case you are advised + to try to simplify the eigenvalues by first finding them using + EIGENVALUES command and then using whatever marvelous tricks you + might have to reduce them to something simpler. You can then use + the KNOWNEIGVALS flag to proceed further. ) + + + +UNITEIGENVECTORS(MAT) takes a MATRIX as its argument and returns a list of + lists the first sublist of which is the output of the EIGENVALUES + command and the other sublists of which are the unit eigenvectors + of the matrix corresponding to those eigenvalues respectively. + The flags mentioned in the description of the EIGENVECTORS command + have the same effects in this one as well. In addition there is one + more flag which may be useful : + + KNOWNEIGVECTS[FALSE] If set to TRUE the EIGEN package will assume that + the eigenvectors of the matrix are known to the user and are stored + under the global name LISTEIGVECTS. LISTEIGVECTS should be set to + a list similar to the output of the EIGENVECTORS command. + ( If KNOWNEIGVECTS is set to TRUE and the list of eigenvectors is + given the setting of the flag NONDIAGONALIZABLE may not be correct. + If that is the case please set it to the correct value. The author + assumes that the user knows what he is doing and will not try to + diagonalize a matrix the eigenvectors of which do not span the + vector space of the appropriate dimension...) + + +SIMILARITYTRANSFORM(MAT) takes a MATRIX as its argument and returns a list + which is the output of the UNITEIGENVECTORS command. In addition if + the flag NONDIAGONALIZABLE is FALSE two global matrices LEFTMATRIX + and RIGHTMATRIX will be generated. These matrices have the property + that LEFTMATRIX.MAT.RIGHTMATRIX is a diagonal matrix with the + eigenvalues of MAT on the diagonal. If NONDIAGONALIZABLE + is TRUE these two matrices will not be generated. + If the flag HERMITIANMATRIX is TRUE then LEFTMATRIX is the + complex conjugate of the transpose of RIGHTMATRIX. Otherwise + LEFTMATRIX is the inverse of RIGHTMATRIX. RIGHTMATRIX + is the matrix the columns of which are the unit eigenvectors of MAT. + The other flags have the same effects since SIMILARITYTRANSFORM + calls the other functions in the package in order to be able to + form RIGHTMATRIX... + + + Finally, for some of you who may think that the names of the + functions are too long, I also made shorter names...( In the following + list " := " means "is equivalent to" ...). Note that using these may + make your code pretty unreadable, you'll save 50% in typing though. + +CONJ(X):= CONJUGATE(X) + +INPROD(X,Y):= INNERPRODUCT(X,Y) + +UVECT(X):= UNITVECTOR(X) + +COVECT(X):= COLUMNVECTOR(X) + +GSCHMIT(X):= GRAMSCHMIDT(X) + +EIVALS(MAT):= EIGENVALUES(MAT) + +EIVECTS(MAT):= EIGENVECTORS(MAT) + +UEIVECTS(MAT):= UNITEIGENVECTORS(MAT) + +SIMTRAN(MAT):= SIMILARITYTRANSFORM(MAT) + + Well, I guess this is the end... Have fun with the EIGEN + package. I hope it will give you useful results. If you run into a + bug please let me know... + ( I do not think there is any, but nevertheless... ) + + + Here is the beginning of Eigenfinite field and Jordanform. + I hope it will give you useful and fascinating results. + --ncs. + +EIGENFINITE(MAT) By brute force plugging in of every member of the reduced + modulus p, p prime, eigenfinite finds the roots of the characteristic + equation generated from CHARPOLY. This is here only in the event + of error checking. MODULUS flag must be set. + +EIGENFINITEFIELD(MAT) Using FACTOR, eigenfinitefield finds the roots of the + characteristic equation generated from CHARPOLY. If entries are + in CRE form, eigenfinitefield, TOTALDISREP's them into general form. + MODULUS flag must be set. + +JORDANBLOCK(X,MAT) X is an eigenvalue of matrix MAT. JORDANBLOCK will return + a list which contains all the Jordan block information for the + eigenvalue X. For example, [0,1,2,3] states that there is one 2x2 + Jordan block, two 3x3 Jordan blocks, three 4x4 Jordan blocks for + the given eigenvalue. JORDANBLOCK calls PRANK. + +PRETTYJORDAN(LIST,MAT) This constructs the Jordan matrix from the + information given by JORDANFORM. + +JORDANFORM(MAT) This first calls EIGENVALUES, then JORDANBLOCK for each + eigenvalue, finally returning a matrix constructed by PRETTYJORDAN. + +JRANK(MAT,LAMBDA) Since the system RANK has zero equivalence problems, + JRANK is RANK with a user given lambda simplifier. For example, + RANK gives the wrong answer for the matrix + [sin(x),1-cos(x)],[cos(x)+1,sin(x)]. But by giving as argument, + lambda([x],fullratsimp(trigreduce(x))) to JRANK the proper answer, + one will be given. A global variable R is left loose by JRANK. + So this program should not be called. Call PRANK instead. + +PRANK(MAT) This program wraps a block around JRANK, so that there are + no loose global variables. The flag RANKPRO will normally + be set false and FULLRATSIMP will be used to determine zero + equivalence. If RANKPRO is set to a lambda expression then + that lambda will be used to simplify for zero equivalence. + + +Any comments, criticisms, or bugs should be addressed to ncs @ mc. \ No newline at end of file diff --git a/src/share/eigen.demo b/src/share/eigen.demo new file mode 100755 index 00000000..98276061 --- /dev/null +++ b/src/share/eigen.demo @@ -0,0 +1,404 @@ +/* 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); + +/* 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 + ... 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... */ + +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 */ + + +/* Part II. Matrices, over a finite field, Zp. If you are a physicist, or + engineer, your exposure to matrices has been limited to matrices + with entries in the reals or the complex domain. You probably know + that spiralling sinks are associated with complex eigenvalues with + length (called modulus) less than one. Spiralling sources similarly + are associated with complex eigenvalues with length >1. + Over a finite field, the geometric intuition is not there, but + what you think should happen does. Eigenvalues and eigenvectors + obey the same definitions. The only problem is that not all the + eigenvalues might lie in the field. + Finite fields can be found useful as a quantum model-- + there are a finite number of energy states, and the energies + cycle back after an electromagnetic emission. This is only + one hypothetical model. I am sure that more physical models + could be found. + In any event, mathematically these are well-defined + quantities, which behave much better than floating point + realities. Either enjoy them as they are, or take it as + a challenge to find them a physical reality!!! + --ncs */ + +/* */ + +reset()$ + +/* A MATRIX THAT I WORKED WITH AND FOUND SOME NICE THEOREMS ABOUT IS + THE PASCAL MATRIX J. + HERE, I WORK WITH THE 3 X 3 MATRIX MODULO 3, AND THE 5 X 5 MATRIX + MODULO 5. */ + +m[i,j]:= binomial(i+j-2,j-1)$ + +pascal[k]:= genmatrix(m,k,k); + +PASCAL[3]; + +jordanform(pascal[3]); + +modulus:3$ + +jordanform(pascal[3]); + +/* By introducing the modulus 3, the matrix pascal[3], + has a totally different Jordan form. Notice the off-diagonal + terms. */ + +/* NOW LET'S TRY PASCAL[5] WITH MODULUS 5. */ + +modulus:5$ + +PASCAL[5]; + +jordanform(pascal[5]); + +/* THERE JUST AREN'T ENOUGH ROOTS IN MODULULO 5. IN A SUITABLE LARGER + FIELD ALL THE ROOTS WOULD LIE. + HERE WE COULD TAKE THAT FIELD TO BE THE REALS BY modulus:false$ + THEN jordanform(pascal[5]); */ + +/* Knowing that the roots are 1, w, and w^2, the cube roots of unity, with */ +/* multiplicities 1, 2, and 2, I use EIGENVECTOR and macsyma's */ +/* ALGEBRAIC INTEGER capability */ + +tellrat(w^2+w+1)$ + +knowneigvals:true$ +listeigvals:[[1,w,w^2],[1,2,2]]$ + +eigenvectors(PASCAL[5]); + +/* Suppose we are interested in Finite Fermat Transforms (FFT). + In modulus 17, the second Fermat number, a fourth root is 4. + So for convolutions of length 4, the following transformation + matrix is important. */ + +kill(m)$ + +m[i,j]:= v^((i-1)*(j-1))$ + +/* HERE IS THE STRUCTURE OF THE FFT. NOTICE IT IS THE SAME FFT + AS COOLEY-TUKEY'S FFT OVER THE COMPLEX FIELD. A relevant + reference is E. Landau "NUMBER THEORY", especially Schur's + proof of Gauss's sum. */ + +genmatrix(m,4,4); + +modulus:5$ + +kill(m)$ + +m[i,j]:=2^((i-1)*(j-1))$ + +a:matrixmap(fullratsimp,genmatrix(m,4,4)); + +jordanform(a); + +/* MUCH MORE INFORMATION CAN BE OBTAINED IN THIS EXAMPLE BY USE OF + THE DOT OPERATOR. THE FOLLOWING IS a^2 */ + +matrixmap(fullratsimp,a.a); + +/* another cute example. The Vandermonde matrix, evaluated at the + fourth root of unity mod 17, and dimensioned + as a 4 x 4 matrix */ + +kill(m)$ + +m[i,j]:= 4^(i-1)$ + +modulus:17$ + +genmatrix(m,4,4); + +jordanform(%); + +/* now isn't that unbelievable!!!! THERE ARE FOUR EIGENVALUES, ALL ZERO, + TWO JORDAN BLOCKS OF SIZE ONE, + AND ONE JORDAN BLOCK OF SIZE TWO. */ + +/* SO THIS IS ANOTHER INTERESTING EXAMPLE. */ + +kill(m)$ + +m[i,j]:= i+j$ + +modulus:5$ + +matrixmap(fullratsimp,genmatrix(m,4,4)); + +jordanform(%); + +/* THE OLD RANK WASN'T A VERY VERSATILE COMMAND. IT COULD TELL THE RANK OF + MANY MATRICES, BUT WHOEVER WROTE + IT EXPECTED THAT THE CONSUMER WOULD ONLY INPUT INTEGERS. THUS RANK WOULD + SOMETIMES GIVE THE WRONG ANSWER... + TO REMEDY THIS ILL, I WROTE THE COMMAND PRANK WHICH USES THE SAME ALGORITHM + AS THE OLD RANK, HOWEVER + IT WRAPS THE ZERO-EQUIVALENCE STEP IN A LAMBDA-WRAPPER. THIS ALLOWS THE + CONSUMER FREEDOM OF CHOICE.*/ + +FOOBAR:MATRIX([SIN(X),1-COS(X)],[COS(X)+1,SIN(X)]); + +/* THIS EXAMPLE DATES FROM 1982. RANK GIVES THE WRONG ANSWER AS CAN BE SEEN */ + +RANK(FOOBAR); + +/* NOW LETS TRY pRANK!! */ + +RANKPRO:LAMBDA([X],FULLRATSIMP(TRIGREDUCE(X)))$ + +PRANK(FOOBAR); + +/* SO BY TWIDDLING WITH rankpro YOU SHOULD HAVE COMPLETE FREEDOM OF CHOICE + AS TO YOUR SIMPLIFIERS. + HOWEVER, IT SHOULD BE MENTIONED THAT D. RICHARDSON PROVED IN 1968 THAT + ZERO-EQUIVALENCE IS + FORMALLY UNSOLVABLE. SEE "SOME UNSOLVABLE PROBLEMS INVOLVING ELEMENTARY + FUNCTIONS OF A REAL VARIABLE" + J. SYMBOLIC LOGIC,33,1968,PP.514-520 */ + +/* So that's all for now*/ \ No newline at end of file