share/matrix/eigen_1.dem

   1 load("eigen")$
   2
   3 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]);
   4
   5 matrix2:matrix([1,2,3,4],[0,3,4,5],[0,0,5,6],[0,0,0,9]);
   6
   7 /* the next function takes a matrix as its argument and returns the
   8         eigenvalues and the unit eigen vectors of that matrix... */
   9
  10 uniteigenvectors(matrix2);
  11
  12 /* if you already know the eigenvectors you can set the flag
  13         knowneigvects to true and the global variable listeigvects to the
  14         list of the eigen vectors...
  15         the next function takes a matrix as its argument and returns the eigen
  16         values and the unit eigen vectors of that matrix.  in addition if
  17         the flag nondiagonalizable is false,two global matrices leftmatrix and
  18         rightmatrix will be generated.  these matrices have the property that
  19         leftmatrix.(matrix).rightmatrix is a diagonal matrix with the eigen
  20         values of the (matrix) on the diagonal...  */
  21
  22 similaritytransform(matrix1)$
  23 nondiagonalizable;
  24 ratsimp(leftmatrix.matrix1.rightmatrix);
  25
  26 /* now that you know how to use the eigen package, here are some
  27         examples about how not to use it.
  28         consider the following matrix :   */
  29
  30 matrix3:matrix([1,0],[0,1]);
  31
  32 /* as you've undoubtedly noticed, this is the 2*2 identity matrix.
  33         let's find the eigen values and the eigen vectors of this matrix...
  34 */
  35
  36 eigenvectors(matrix3);
  37
  38 /* "nothing special happened", you say.  everyone knows what the eigen
  39         values and the eigen vectors of the identity matrix are, right?
  40         right.  now consider the following matrix :  */
  41
  42 matrix4:matrix([1,e],[e,1]);
  43
  44 /* let e>0, but as small as you can imagine.  say 10^(-100).
  45         let's find the eigen values and the eigen vectors of this matrix :
  46 */
  47
  48 eigenvectors(matrix4);
  49
  50 /* since e~10^(-100), the eigen values of matrix4 are equal to the
  51         eigen values of matrix3 to a very good accuracy.  but, look
  52         at the eigen vectors!!!  eigen vectors of matrix4 are nowhere
  53         near the eigen vectors of matrix3.  there is angle of %pi/4
  54         between the corresponding eigen vectors.  so, one learns
  55         another fact of life :
  56
  57         matrices which have approximately the same eigen values do not
  58         have approximately the same eigen vectors in general.
  59
  60         this example may seem artificial to you, but it is not if you think
  61         a little bit more about it.  so, please be careful when you
  62         approximate the entries of whatever matrix you have.  you may
  63         get good approximations to its eigen values, however the eigen
  64         vectors you get may be entirely spurious( or some may be correct,
  65         but some others may be totally wrong ).
  66
  67         now, here is another sad story :
  68         let's take a look at the following matrix :  */
  69
  70 matrix5:matrix([5/2,50-25*%i],[50+25*%i,2505/2]);
  71
  72 /* nice looking matrix, isn't it?  as usual, we will find the eigen
  73         values and the eigen vectors of it :  */
  74
  75 eigenvectors(matrix5);
  76
  77 /* well, here they are.  suppose that this was not what you wanted.
  78         instead of those sqrt(70)'s, you want the numerical values of
  79         everything.  one way of doing this is to set the flag "numer"
  80         to true and use the eigenvectors command again :  */
  81
  82 numer:true;
  83 eigenvectors(matrix5);
  84
  85 /* ooops!!!  what happened??  we got the eigen values, but there are
  86         no eigenvectors.  nonsense, there must be a bug in eigen, right?
  87         wrong.  there is no bug in eigen.  we have done something which
  88         we should not have done.  let me explain :
  89         when one is solving for the eigen vectors, one has to find the
  90         solution to homogeneous equations like :  */
  91
  92 equation1:a*x+b*y=0;
  93 equation2:c*x+d*y=0;
  94
  95 /* in order for this set of equations to have a solution other than
  96         the trivial solution ( the one in which x=0 and y=0 ), the
  97         determinant of the coefficients ( in this case a*d-b*c ) should
  98         vanish.  exactly.   if the determinant does not vanish the only
  99         solution will be the trivial solution and we will get no eigen
 100         vectors.  during this demo, i did not set a,b,c,d to any
 101         particular values.  let's see what happens when we try to solve
 102         the set above :  */
 103
 104 algsys([equation1,equation2],[x,y]);
 105
 106 /* you see?  the infamous trivial solution.  now let me set a,b,c,d
 107         to some numerical values :  */
 108
 109 a:4;
 110 b:6;
 111 c:2;
 112 d:3;
 113 a*d-b*c;
 114 equation1:ev(equation1);
 115 equation2:ev(equation2);
 116 algsys([equation1,equation2],[x,y]);
 117
 118 /* now we have a nontrivial solution with one arbitrary constant.
 119         ( %r(something) ).  what happened in the previous case is that
 120         the numerical errors caused the determinant not to vanish, hence
 121         algsys gave the trivial solution and we got no eigen vectors.
 122         if you want a numerical answer, first calculate it exactly,
 123         then set "numer" to true and evaluate the answer.  */
 124
 125 numer:false;
 126 notnumerical:eigenvectors(matrix5);
 127 numer:true;
 128 ev(notnumerical);
 129
 130 /* you see, it works now.  actually, if you have a matrix with
 131         numerical entries and you can live with reasonably accurate
 132         answers, there are much better (faster) programs.  ask somebody
 133         about the imsl routines on the share directory...
 134         this is all...  if you think that the names of the functions are too
 135         long, there are shorter names for them and they are given in the file
 136         eigen usage dsk:share;.  good luck!!!!!!!!!!!!!...... yekta  */
 137
 138
 139 "end of demo"$