share/linearalgebra/rtest_matrixexp.mac

   1 (load("matrixexp.lisp"),0);
   2 0$
   3
   4 (remvalue(a,b,c,x),0);
   5 0$
   6
   7 (xequal(a,b) := block([gcd : 'spmod, algebraic : true],
   8    ratvars(),
   9    is(fullratsimp(a) = fullratsimp(b)) or
  10       is(fullratsimp(rectform(a)) = fullratsimp(rectform(b)))),0);
  11 0$
  12
  13 scalarmatrixp : false;
  14 false$
  15
  16 (check(mat,n) := block([s, okay, z : gensym(), i : 0, m],
  17
  18   s : spectral_rep(mat),
  19   okay : xequal(mat - s[1] . s[2] - s[3] , zeromatrix(n,n)),
  20   print(i : i + 1, okay),
  21
  22   okay : okay and xequal(mat . s[3], s[3] . mat),
  23   print(i : i + 1, okay),
  24
  25   okay : okay and xequal(s[3]^^n, zeromatrix(n,n)),
  26   print(i : i + 1, okay),
  27
  28   okay : okay and xequal(matrixfun(lambda([x],1), mat), ident(n)),
  29   print(i : i + 1, okay),
  30
  31   okay : okay and xequal(matrixfun(lambda([x], 1+x+x^2), mat), ident(n) + mat + mat . mat),
  32   print(i : i + 1, okay),
  33
  34   okay : okay and xequal(matrixfun(lambda([x], x^2/5), mat), mat . mat / 5),
  35   print(i : i + 1, okay),
  36
  37   okay : okay and xequal(matrixfun(lambda([x],x^8),mat), mat^^8),
  38   print(i : i + 1, okay),
  39
  40   m : errcatch(matrixfun(lambda([z],sqrt(z)), mat)),
  41   okay : okay and if m = [ ] then member(0, s[1]) else xequal(m[1] . m[1], mat),
  42   print(i : i + 1, okay),
  43
  44   m : errcatch(matrixfun(lambda([mat], 1/mat^2), mat)),
  45   okay : okay and if m = [ ] then member(0, s[1]) else xequal(m[1], mat^^-2),
  46   print(i : i + 1, okay),
  47
  48   m : fullratsimp(matrixexp(mat,z) . matrixexp(mat, -z)),
  49   okay : okay and xequal(m, ident(n)),
  50   print(i : i + 1, okay),
  51
  52   m : matrixexp(mat,z),
  53   okay : okay and xequal(diff(m,z) - m . mat, zeromatrix(n,n)),
  54   print(i : i + 1, okay),
  55
  56   okay), 0);
  57 0$
  58
  59 check(matrix([1]),1);
  60 true$
  61
  62 check(matrix([x]),1);
  63 true$
  64
  65 check(matrix([1,0],[0,1]),2);
  66 true$
  67
  68 check(matrix([1,0],[0,0]),2);
  69 true$
  70
  71 check(matrix([0,0],[0,1]),2);
  72 true$
  73
  74 check(matrix([0,1],[1,0]),2);
  75 true$
  76
  77 check(matrix([0,1],[-1,0]),2);
  78 true$
  79
  80 check(matrix([0,%i],[-%i,0]),2);
  81 true$
  82
  83 check(matrix([0,a],[0,b]),2);
  84 true$
  85
  86 check(matrix([1,2],[2,1]),2);
  87 true$
  88
  89 check(matrix([1,%i],[-%i,1]),2);
  90 true$
  91
  92 check(matrix([1,2],[3,1]),2);
  93 true$
  94
  95 check(matrix([1,6],[0,1]),2);
  96 true$
  97
  98 check(matrix([x,1],[1,1]), 2);
  99 true$
 100
 101 check(matrix([sqrt(3), sqrt(2)],[sqrt(2), 9]),2);
 102 true$
 103
 104 check(matrix([sqrt(6),1],[1, -9]),2);
 105 true$
 106
 107 check(matrix([sqrt(3), sqrt(2)],[sqrt(2), 9]),2);
 108 true$
 109
 110 check(matrix([a,b],[c,d]),2);
 111 true$
 112
 113 check(matrix([1,0,0],[0,1,0],[0,0,1]),3);
 114 true$
 115
 116 check(matrix([1,0,0],[0,1,0],[0,0,-1]),3);
 117 true$
 118
 119 check(matrix([1,1,0],[1,1,1],[0,1,1]),3);
 120 true$
 121
 122 check(matrix([1,%i,0],[-%i,1,%i],[0,-%i,1]),3);
 123 true$
 124
 125 check(matrix([1,2,3],[4,5,6],[7,8,9]),3);
 126 true$
 127
 128 check(matrix([1,a,b],[0,1,c],[0,0,2]),3);
 129 true$
 130
 131 check(matrix([1,a,b],[0,1,c],[0,0,1]),3);
 132 true$
 133
 134 check(matrix([1,a,b],[0,2,c],[0,0,1]),3);
 135 true$
 136
 137 check(matrix([2,a,b],[0,2,c],[0,0,1]),3);
 138 true$
 139
 140 check(matrix([1,a,b],[0,1,c],[0,0,2]),3);
 141 true$
 142
 143 check(matrix([1,0,0],[a,1,0],[b,c,1]),3);
 144 true$
 145
 146 check(matrix([1,0,0],[a,%i,0],[b,c,%i]),3);
 147 true$
 148
 149 check(matrix([1,0,0],[a,%i,0],[b,c,1]),3);
 150 true$
 151
 152 check(matrix([1,0,0],[a,%i,0],[b,c,sqrt(7)]),3);
 153 true$
 154
 155 check(matrix([1,0,0],[a,%i,0],[b,c,%i]),3);
 156 true$
 157
 158 /* spectral_rep checks its results--so this really does testing.*/
 159
 160 (spectral_rep(matrix([x,1,0],[1,1,1],[0,1,1])),0);
 161 0$
 162
 163
 164 check(matrix([1,1,0],[1,1,1],[0,1,1]),3);
 165 true$
 166
 167 (s : matrix([1,2],[2,3]), 0);
 168 0$
 169
 170 (m : s . matrix([1,0],[0,-1]) . s^^-1, 0);
 171 0$
 172
 173 xequal(s . matrix([exp(1),0],[0,exp(-1)]) . s^^-1, matrixexp(m));
 174 true$
 175
 176
 177 /* Matrices from Gosei Furuya's Jordan form code. */
 178
 179 (m : matrix([2,0,0,0,0,0,0,0],
 180          [1,2,0,0,0,0,0,0],
 181          [-4,1,2,0,0,0,0,0],
 182          [2,0,0,2,0,0,0,0],
 183          [-7,2,0,0,2,0,0,0],
 184          [9,0,-2,0,1,2,0,0],
 185          [-34,7,1,-2,-1,1,2,0],
 186          [145,-17,-16,3,9,-2,0,3]), 0)$
 187 0$
 188
 189 xequal(matrixfun(lambda([x],x^2), m), m.m);
 190 true$
 191
 192
 193 (m : matrix([2,1,0,0,0,0],
 194           [-1,4,0,0,0,0],
 195           [-1,1,2,1,0,0],
 196           [-1,1,-1,4,0,0],
 197           [-1,1,-1,1,3,0],
 198           [-1,1,-1,1,1,2]),0)$
 199 0$
 200
 201 xequal(matrixfun(lambda([x],x^2), m), m.m);
 202 true$
 203
 204 (m : matrix([0,1,0],
 205           [0,0,1],
 206           [1,-3,3]),0);
 207 0$
 208
 209 xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5);
 210 true$
 211
 212 (m : matrix([0,1,2],
 213           [0,0,0],
 214           [0,0,0]),0);
 215 0$
 216
 217 xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5);
 218 true$
 219
 220 (m : matrix([3,1,0,0],
 221           [-4,-1,0,0],
 222           [7,1,2,1],
 223           [-17,-6,-1,0]),0);
 224 0$
 225
 226 xequal(matrixfun(lambda([x],1+x^2), m), ident(4)+m^^2);
 227 true$
 228
 229 (m : matrix([2,1,2,0],
 230           [-2,2,1,2],
 231           [-2,-1,-1,1],
 232           [3,1,2,-1]),0);
 233 0$
 234
 235 xequal(matrixfun(lambda([x],1+x^2), m), ident(4)+m^^2);
 236 true$
 237
 238 (m : matrix([0,0,1,1,1],
 239           [0,0,0,1,1],
 240           [0,0,0,0,1],
 241           [0,0,0,0,0],
 242           [0,0,0,0,0]),0);
 243 0$
 244
 245 xequal(matrixfun(lambda([x],1+x^2), m), ident(5)+m^^2);
 246 true$
 247
 248 (m : matrix([0,1,0],
 249           [0,0,1],
 250           [-1,-3,-3]),0);
 251 0$
 252
 253 xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5);
 254 true$
 255
 256
 257
 258