share/matrix/eigen.mac

   1 /* modified for doe macsyma with define_variable */
   2
   3 /*      this is the file eigen > dsk:share;.
   4         this is the source code for the new eigen package and it is
   5         macsyma batchable, i.e. batch(eigen,>,dsk,share);.  if you do not want
   6         to waste time (and/or paper...) the fastloadable version is on the file
   7         eigen fasl dsk:share;.  you can load the latter using macsyma's
   8         loadfile command, i.e. loadfile(eigen,fasl,dsk,share);. the functions
   9         are described in the file eigen usage dsk:share;, and the demo file in
  10         which the functions are demonstrated is eigen demo dsk:share;.     */
  11
  12 /* commented out of doe macsyma
  13 eval_when(translate_file,
  14         modedeclare([hermitianmatrix,nondiagonalizable,knowneigvals,
  15         knowneigvects],boolean,
  16         [index1,index2,index3,index4,dimnsn,count,%rnum],fixnum),
  17         declare([hermitianmatrix,nondiagonalizable,knowneigvals,
  18         knowneigvects,listeigvects,listeigvals,%rnum,listarith,programmode,
  19         algebraic,realonly,multiplicities,solveexplicit],special))$ */
  20
  21 eval_when([translate,batch,demo,load,loadfile],
  22 mi(var)::=buildq([var],mode_identity(fixnum,var)),
  23 dv(var)::=buildq([var],define_variable(var,false,boolean)))$
  24
  25 /* commented out of doe macsyma
  26 sstatus(feature,eigen)$
  27
  28 hermitianmatrix:false$
  29
  30 nondiagonalizable:false$
  31
  32 knowneigvals:false$
  33
  34 knowneigvects:false$
  35
  36 listeigvects:[]$
  37
  38 listeigvals:[]$
  39 */
  40
  41 dv(hermitianmatrix)$ dv(nondiagonalizable)$ dv(knowneigvals)$
  42 dv(knowneigvects)$
  43
  44 define_variable(listeigvects,[],list)$
  45 define_variable(listeigvals,[],list)$
  46 define_variable(rightmatrix,[],any)$
  47 define_variable(leftmatrix,[],any)$
  48
  49 /* conjugate(x):=sublis('[%i=-%i],x)$ */
  50
  51 innerproduct(x,y):=block([listarith],listarith:true,ratsimp(conjugate(x).y))$
  52 /*
  53 unitvector(x):=block([listarith,intrn],listarith:true,intrn:innerproduct(x,x),
  54 intrn:sqrt(intrn),x/intrn)$
  55 */
  56 unitvector(x):=block([listarith],listarith:true,x/sqrt(innerproduct(x,x)))$
  57
  58 columnvector(x):=transpose(matrix(x))$
  59
  60
  61 gramschmidt(x%,[myinnerproduct]):=
  62                 block([listarith,dimnsn,listall,intern,count,denom,unit,index1,
  63                 index2],
  64                 if myinnerproduct=[] then myinnerproduct:innerproduct
  65                 else myinnerproduct:first(myinnerproduct),
  66                 mode_declare([dimnsn,count,index1,index2],fixnum,
  67                               [listall],list,[intern,denom,unit],any),
  68                 listarith:true,dimnsn:mi(length(x%)),listall:[part(x%,1)],
  69                 count:1,if dimnsn=1 then return(x%)
  70                 else (for index1:2 thru dimnsn do
  71                 (unit:part(x%,index1),for index2 thru count do
  72                 (intern:part(listall,index2),denom : apply (myinnerproduct, [intern, intern]),
  73                 unit:factor(ratsimp(unit - apply (myinnerproduct, [intern, unit])*intern/denom
  74                 ))),
  75                 count:count+1,listall:endcons(unit,listall)),
  76                 return(listall)))$
  77
  78
  79 eigenvalues(mat):=
  80                 block([dimnsn,listall,solution,multiplicities,solveexplicit,
  81                 dummy:?gensym(),index2],
  82                 mode_declare([dimnsn,index2],fixnum,[listall,solution],list,
  83                              [dummy],any),
  84                 listall:[],
  85                 solveexplicit:true,
  86                 dimnsn:mi(length(mat)),
  87                 solution:block([programmode:true],
  88                 solve(charpoly(mat,dummy),dummy)),
  89                 if solution=[] then
  90                 (print(" "),print("eigenvalues: solve is unable to find the roots of the characteristic polynomial."),
  91                 return(listall))
  92                 else if apply("+",multiplicities) < dimnsn then
  93                 (print(" "),print("eigenvalues: solve is unable to find some of the roots of the characteristic polynomial."),
  94                 dimnsn:apply("+",multiplicities)),
  95                 for index2 thru dimnsn do
  96                 (dimnsn:mi(dimnsn-part(multiplicities,index2)+1),
  97                 listall:endcons(rhs(part(solution,index2)),listall)),
  98                 listall:endcons(multiplicities,[listall]),
  99                 return(listall))$
 100
 101 /* EIGENVECTORS - Compute eigenvectors of the matrix mat.
 102
 103   The return value is a list of two elements. The first is a list of the
 104   eigenvalues of mat and a list of the multiplicities of the eigenvalues,
 105   found by calling eigenvalues().  The second is a list of lists of
 106   eigenvectors. There is one list of eigenvectors for each eigenvalue.
 107   There may be one or more eigenvectors in each list.
 108
 109   The following flags are used:
 110
 111   nondiagonalizable is set to true when the matrix is nondiagonalizable
 112   after eigenvectors returns, or false otherwise.
 113
 114   hermitianmatrix when true, causes the degenerate eigenvectors of the
 115   Hermitian matrix to be orthogonalized using the Gram-Schmidt algorithm.
 116
 117   knowneigvals when true causes the eigen package to assume the eigenvalues
 118   of the matrix are known to the user and stored under the global name
 119   listeigvals. listeigvals should be set to a list similar to the output
 120   eigenvalues.
 121  */
 122 eigenvectors(mat):=
 123                 block([eigvecs,equations,unknowns,solution,eigvals,dimnsn,
 124                 count,vectr,index3,index4,index2,index1,matrx,mmatrx,notknwn,
 125                 unit,multiplicities,%rnum,realonly,algebraic,interm,intern],
 126                 local(solution, mmatrx),
 127                 mode_declare([equations,unknowns,solution,eigvals,
 128                               unit,interm],list,
 129                              [dimnsn,count,index3,index4,index2,index1],fixnum,
 130                              [vectr,matrx,mmatrx,intern,notknwn],any),
 131                 unknowns:[],dimnsn:mi(length(mat)),
 132                 count:mi(dimnsn),notknwn:?gensym(),
 133                 mat:rationalize(mat),
 134                 if knowneigvals then eigvals:listeigvals
 135                 else eigvals:eigenvalues(mat),
 136                 if eigvals=[] then (nondiagonalizable:true,return(eigvals))
 137                 else (multiplicities:part(eigvals,2),
 138                 count:length(multiplicities),
 139                 for index1 thru dimnsn do
 140                 unknowns:endcons(concat(notknwn,index1),unknowns),
 141                 vectr:columnvector(unknowns),matrx:mat.vectr,
 142                 nondiagonalizable:false,
 143                 eigvecs:[],realonly:false,algebraic:true,
 144                 for index1 thru count do
 145                 (mmatrx:matrx-part(eigvals,1,index1)*vectr,
 146                 equations:[],
 147                 for index2 thru dimnsn do
 148                 equations:cons(mmatrx[index2,1],equations),%rnum:0,
 149                 solution:algsys(equations,unknowns),
 150                 interm:map('rhs,solution[1]),
 151                 unit:[],if %rnum#part(multiplicities,index1)
 152                 then nondiagonalizable:true,
 153                 for index3 thru %rnum do
 154                 (intern:substvectk(%rnum_list,index3,interm),
 155                 unit:append(unit,[intern])),
 156                 if unit=[] then
 157                 (print(" "),print("eigenvectors: the eigenvector(s) for the",
 158                 index1,"th eigenvalue will be missing.")),
 159                 if hermitianmatrix and %rnum>1 then unit:gramschmidt(unit),
 160                 eigvecs:append(eigvecs,[unit])),
 161                 return([eigvals, eigvecs])))$
 162
 163
 164 /* the first arg is of the form [r1,r2,r3].
 165    we want to construct [r1=0,r2=1,r3=0] for example. */
 166
 167 substvectk(l,n,exp):=(mode_declare(l,list,n,fixnum,exp,any),
 168                 block([sub_list:[],j:0],mode_declare(sub_list,list,j,fixnum),
 169                 for var in l do (mode_declare(var,any),
 170                  j:j+1,sub_list:cons(var = if j=n then 1 else 0,sub_list)),
 171                 sublis(sub_list,exp)))$
 172
 173
 174 uniteigenvectors (M):= block ([uv], local(uv),
 175     mode_declare (uv, list, [i, j], fixnum),
 176
 177     if knowneigvects
 178         then uv : copy (listeigvects)
 179         else uv : copy (eigenvectors (M)),
 180
 181     if uv = []
 182         then []
 183         else
 184            (for i thru length (uv[2])
 185                 do for j thru length (uv[2][i])
 186                     /* modification in-place (OK due to copy above) */
 187                     do uv[2][i][j] : ratsimp (unitvector (uv[2][i][j])),
 188             uv));
 189
 190
 191 similaritytransform(mat):=
 192                 block([listvec,listuevec], local(listuevec),
 193                 mode_declare([listvec,listuevec],list),
 194                 listuevec:uniteigenvectors(mat),
 195                 if nondiagonalizable then return(listuevec)
 196                 else (listvec: apply (append, listuevec[2]),
 197                 rightmatrix:transpose(apply('matrix,listvec)),
 198                 if hermitianmatrix then
 199                 leftmatrix:conjugate(transpose(rightmatrix))
 200                 else leftmatrix:rightmatrix^^-1,
 201                 return(listuevec)))$
 202
 203
 204 conj(x):=conjugate(x)$
 205
 206 inprod(x,y):=innerproduct(x,y)$
 207
 208 uvect(x):=unitvector(x)$
 209
 210 covect(x):=columnvector(x)$
 211
 212 gschmit(x%, [f]):= if f = [] then gramschmidt(x%) else gramschmidt (x%, first(f))$
 213
 214 eivals(mat):=eigenvalues(mat)$
 215
 216 eivects(mat):=eigenvectors(mat)$
 217
 218 ueivects(mat):=uniteigenvectors(mat)$
 219
 220 simtran(mat):=similaritytransform(mat)$
 221