share/matrix/eigen.dem

   1 /* modified for doe macsyma */
   2 /* this is the file eigen demo dsk:share;.
   3  (you can batch or demo this file, i.e. batch(eigen,demo,dsk,share);, or
   4  demo(eigen,demo,dsk,share);.  note that in the demo mode you have to hit the
   5  space key after each step...)
   6  the functions in the new eigen package are demonstrated here. the  description
   7  of the functions can be found in the file eigen usage dsk:share;, the
   8  source code is on the file eigen > dsk:share; and the fastload file is
   9  eigen fasl dsk:share;.  ( you can load this one  using macsyma's loadfile
  10  command, i.e. loadfile(eigen,fasl,dsk,share);.)
  11
  12
  13         we start with loading the eigen package :      */
  14 /*
  15 if not status(feature,eigen) then loadfile(eigen,fasl,dsk,share);
  16 */
  17 load("eigen")$
  18 /* let us start with the first function.  (see the descriptions...)
  19         first let's define a complex variable... */
  20
  21
  22 z:a+%i*b;
  23
  24 /* the conjugate function simply returns the complex conjugate of its
  25         argument...    */
  26
  27
  28 conjugate(z);
  29
  30 /* note that z could be a matrix, a list, etc...   */
  31
  32 z:matrix([%i,0],[0,1+%i]);
  33 conjugate(z);
  34
  35 /* the next function calculates the inner  product of two lists...*/
  36
  37 list1:[a,b,c,d];
  38 list2:[f,g,h,i];
  39 innerproduct(list1,list2);
  40
  41 /* the elements of the lists could be complex also...  */
  42
  43 list1:[a+%i*b,c+%i*d];
  44 innerproduct(list1,list1);
  45
  46 /* the next function takes a list as its argument and returns a unit
  47         list...   */
  48
  49 list1:[1,1,1,1,1];
  50 unitvector(list1);
  51 list2:[1,%i,1-%i,1+%i];
  52 unitvector(list2);
  53
  54 /* the next function takes a list as its argument and returns a column
  55         vector...      */
  56
  57 list1:[a,b,c,d];
  58 columnvector(list1);
  59
  60 /* the next function takes a list of lists as its argument and
  61         orthogonalizes them using the gram-schmidt algorithm...*/
  62
  63 listoflists:[[1,2,3,4],[0,5,4,7],[4,5,6,7],[0,0,1,0]];
  64
  65 /* note that the lists in this list of lists are not orthogonal to each
  66         other...   */
  67
  68 innerproduct([1,2,3,4],[0,5,4,7]);
  69 innerproduct([1,2,3,4],[4,5,6,7]);
  70
  71 /* but after applying the gramschmidt function...  */
  72
  73 orthogonallists:gramschmidt(listoflists);
  74 innerproduct(part(orthogonallists,1),part(orthogonallists,2));
  75 innerproduct(part(orthogonallists,2),part(orthogonallists,3));
  76
  77 /* note that orhtogonallists contains integers that are factored.
  78         if you do not like this form, you can simply ratsimp the result :  */
  79
  80 ratsimp(orthogonallists);
  81
  82 /* the next function takes a matrix as its argument and returns the
  83         eigenvalues of that matrix...    */
  84
  85 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]);
  86
  87 /* this is the matrix that caused a lot of trouble for the old eigen
  88         package...  it took ~170 seconds to find the eigen vectors of this
  89         matrix...  you should be able to do it in your head in about 20 seconds
  90         [note: pdp-10 timings]
  91         ...  the new eigen package handles it in about 10 seconds...  anyway,
  92         let's keep going... */
  93
  94
  95 eigenvalues(matrix1);
  96
  97 /* the first sublist in the answer is the eigenvalues, second list is
  98         their multiplicities in the corresponding order...
  99         the next function takes a matrix as its argument and returns the
 100         eigen values and the eigen vectors of that matrix...  */
 101
 102 eigenvectors(matrix1);
 103
 104 /* first sublist in the answer is the output of the eigenvalues command
 105         the others are the eigen vectors corresponding to those eigen values...
 106         notice that this command is more powerful than the eigenvalues command
 107         because it determines both the eigen values and the eigen vectors...
 108         if you already know the eigen values, you can set the knowneigvals flag
 109         to true and the global variable listeigvals to the list of eigen
 110         values...  this will make the execution of eigenvectors command faster
 111         because it doesn't have to find the eigen values itself... */
 112 /* commented out here and placed in [share]eigen.dm1 because lack of gc
 113 in 8.11 version of doe macsyma
 114 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]);
 115
 116 matrix2:matrix([1,2,3,4],[0,3,4,5],[0,0,5,6],[0,0,0,9]);
 117
 118 /* the next function takes a matrix as its argument and returns the
 119         eigenvalues and the unit eigen vectors of that matrix... */
 120
 121 uniteigenvectors(matrix2);
 122
 123 /* if you already know the eigenvectors you can set the flag
 124         knowneigvects to true and the global variable listeigvects to the
 125         list of the eigen vectors...
 126         the next function takes a matrix as its argument and returns the eigen
 127         values and the unit eigen vectors of that matrix.  in addition if
 128         the flag nondiagonalizable is false,two global matrices leftmatrix and
 129         rightmatrix will be generated.  these matrices have the property that
 130         leftmatrix.(matrix).rightmatrix is a diagonal matrix with the eigen
 131         values of the (matrix) on the diagonal...  */
 132
 133 similaritytransform(matrix1)$
 134 nondiagonalizable;
 135 ratsimp(leftmatrix.matrix1.rightmatrix);
 136
 137 /* now that you know how to use the eigen package, here are some
 138         examples about how not to use it.
 139         consider the following matrix :   */
 140
 141 matrix3:matrix([1,0],[0,1]);
 142
 143 /* as you've undoubtedly noticed, this is the 2*2 identity matrix.
 144         let's find the eigen values and the eigen vectors of this matrix...
 145 */
 146
 147 eigenvectors(matrix3);
 148
 149 /* "nothing special happened", you say.  everyone knows what the eigen
 150         values and the eigen vectors of the identity matrix are, right?
 151         right.  now consider the following matrix :  */
 152
 153 matrix4:matrix([1,e],[e,1]);
 154
 155 /* let e>0, but as small as you can imagine.  say 10^(-100).
 156         let's find the eigen values and the eigen vectors of this matrix :
 157 */
 158
 159 eigenvectors(matrix4);
 160
 161 /* since e~10^(-100), the eigen values of matrix4 are equal to the
 162         eigen values of matrix3 to a very good accuracy.  but, look
 163         at the eigen vectors!!!  eigen vectors of matrix4 are nowhere
 164         near the eigen vectors of matrix3.  there is angle of %pi/4
 165         between the corresponding eigen vectors.  so, one learns
 166         another fact of life :
 167
 168         matrices which have approximately the same eigen values do not
 169         have approximately the same eigen vectors in general.
 170
 171         this example may seem artificial to you, but it is not if you think
 172         a little bit more about it.  so, please be careful when you
 173         approximate the entries of whatever matrix you have.  you may
 174         get good approximations to its eigen values, however the eigen
 175         vectors you get may be entirely spurious( or some may be correct,
 176         but some others may be totally wrong ).
 177
 178         now, here is another sad story :
 179         let's take a look at the following matrix :  */
 180
 181 matrix5:matrix([5/2,50-25*%i],[50+25*%i,2505/2]);
 182
 183 /* nice looking matrix, isn't it?  as usual, we will find the eigen
 184         values and the eigen vectors of it :  */
 185
 186 eigenvectors(matrix5);
 187
 188 /* well, here they are.  suppose that this was not what you wanted.
 189         instead of those sqrt(70)'s, you want the numerical values of
 190         everything.  one way of doing this is to set the flag "numer"
 191         to true and use the eigenvectors command again :  */
 192
 193 numer:true;
 194 eigenvectors(matrix5);
 195
 196 /* ooops!!!  what happened??  we got the eigen values, but there are
 197         no eigenvectors.  nonsense, there must be a bug in eigen, right?
 198         wrong.  there is no bug in eigen.  we have done something which
 199         we should not have done.  let me explain :
 200         when one is solving for the eigen vectors, one has to find the
 201         solution to homogeneous equations like :  */
 202
 203 equation1:a*x+b*y=0;
 204 equation2:c*x+d*y=0;
 205
 206 /* in order for this set of equations to have a solution other than
 207         the trivial solution ( the one in which x=0 and y=0 ), the
 208         determinant of the coefficients ( in this case a*d-b*c ) should
 209         vanish.  exactly.   if the determinant does not vanish the only
 210         solution will be the trivial solution and we will get no eigen
 211         vectors.  during this demo, i did not set a,b,c,d to any
 212         particular values.  let's see what happens when we try to solve
 213         the set above :  */
 214
 215 algsys([equation1,equation2],[x,y]);
 216
 217 /* you see?  the infamous trivial solution.  now let me set a,b,c,d
 218         to some numerical values :  */
 219
 220 a:4;
 221 b:6;
 222 c:2;
 223 d:3;
 224 a*d-b*c;
 225 equation1:ev(equation1);
 226 equation2:ev(equation2);
 227 algsys([equation1,equation2],[x,y]);
 228
 229 /* now we have a nontrivial solution with one arbitrary constant.
 230         ( %r(something) ).  what happened in the previous case is that
 231         the numerical errors caused the determinant not to vanish, hence
 232         algsys gave the trivial solution and we got no eigen vectors.
 233         if you want a numerical answer, first calculate it exactly,
 234         then set "numer" to true and evaluate the answer.  */
 235
 236 numer:false;
 237 notnumerical:eigenvectors(matrix5);
 238 numer:true;
 239 ev(notnumerical);
 240
 241 /* you see, it works now.  actually, if you have a matrix with
 242         numerical entries and you can live with reasonably accurate
 243         answers, there are much better (faster) programs.  ask somebody
 244         about the imsl routines on the share directory...
 245         this is all...  if you think that the names of the functions are too
 246         long, there are shorter names for them and they are given in the file
 247         eigen usage dsk:share;.  good luck!!!!!!!!!!!!!...... yekta  */
 248 */
 249
 250 "end of demo"$