share/matrix/rtest_eigen.mac

   1 (kill (all), load (eigen), 0);
   2 0;
   3
   4 (x : matrix ([1, 2, 3], [9, 18, 30], [12, 48, 60]),
   5  y : gramschmidt (x),
   6  expand (y, 0, 0));
   7 ''([[1,2,3],[-3^2/(2*7),-3^2/7,3*5/(2*7)],[-2^4*3/5,2^3*3/5,0]]);
   8
   9 map (innerproduct, [y[1], y[2], y[3]], [y[2], y[3], y[1]]);
  10 [0, 0, 0];
  11
  12 (kill (x),
  13  ip (f, g) := integrate (f * g, 'x, a, b),
  14  ev (y : gramschmidt ([1, sin(x), cos(x)], ip), a= -%pi/2, b=%pi/2));
  15 [1, sin(x), (%pi*cos(x) - 2)/%pi];
  16
  17 map (ip, [y[1], y[2], y[3]], [y[2], y[3], y[1]]), a= -%pi/2, b=%pi/2;
  18 [0, 0, 0];
  19
  20 (load("nchrpl"),0);
  21 0$
  22
  23 mattrace(matrix());
  24 0$
  25
  26 mattrace(matrix([]));
  27 0$
  28
  29 mattrace(matrix([a,b],[c,d]));
  30 a + d$
  31
  32 mattrace(matrix([a,b],[c,d],[p,q]));
  33 a + d$
  34
  35 /* Regression test for SF bug 2721670 */
  36 (a[1,1] : 42, mattrace(matrix([1,2],[5,6])));
  37 7$
  38
  39 (remarray(a),0);
  40 0$
  41
  42 /* from mailing list 2009-08-17 "bugs on wxMaxima 5.19.1" */
  43
  44 block ([D1, uv], local (dmatrix),
  45  /* diag_matrix and diag involve loading big blobs of stuff, and diagmatrix isn't suitable ... sigh. */
  46  dmatrix (L) := block ([M, n : length (L)], M : zeromatrix (n, n), for i thru n do M[i, i] : L[i], M),
  47  D1 : matrix ([3, 2, 1, -5], [-1, 2, 1, -3], [5, -1, 2, 4], [3, 1, -1, 5]),
  48  uv : block([ratmx : true], similaritytransform (D1)),
  49  /* ratsimp or radcan could help here, but ratsimp simplifies leftmatrix . D1 . rightmatrix
  50   * incorrectly ... sigh. Also mat_norm could help but it's embedded in a big blob ... sigh.
  51   */
  52  float (dmatrix (uv[1][1]) - leftmatrix . D1 . rightmatrix),
  53  expand (%%),
  54  matrixmap (cabs, %%),
  55  maplist (lmax, %%),
  56  lmax (%%),
  57  is (%% < 1e-10));
  58 true;
  59
  60 /* from mailing list 2019-05-11: "Problems with gramschmidt" */
  61
  62 (kill(x),
  63  basis: makelist (x^k, k, 0, 5),
  64  ip(u,v) := integrate (u*v, x, -1, 1),
  65  foo: ratsimp (gramschmidt (basis, ip)));
  66 [1,x,(3*x^2-1)/3,(5*x^3-3*x)/5,(35*x^4-30*x^2+3)/35,(63*x^5-70*x^3+15*x)/63]$
  67
  68 bar: map (lambda ([e], ip(e, e)), foo);
  69 [2,2/3,8/45,8/175,128/11025,128/43659]$
  70
  71 (baz: map (lambda ([a, b], a/sqrt(b)), foo, bar),
  72  genmatrix (lambda ([i, j], ip (baz[i], baz[j])), 5, 5));
  73 matrix ([1, 0, 0, 0, 0],
  74         [0, 1, 0, 0, 0],
  75         [0, 0, 1, 0, 0],
  76         [0, 0, 0, 1, 0],
  77         [0, 0, 0, 0, 1]);
  78
  79 /* SF bug #3654: "uniteigenvectors fails if uv[1] is used in the user main program." */
  80
  81 (kill(m, mm, r, uv, foo, theta),
  82  assume (m > 0, r > 0),
  83  mm:matrix(
  84           [-(r-2*m)/r,    sqrt(2*m/r),    0,  0],
  85           [sqrt(2*m/r),   1,  0,  0],
  86           [0, 0,  r^2,    0],
  87           [0, 0,  0,  r^2*sin(theta)^2]
  88       ),
  89  foo: uniteigenvectors(mm),
  90  0);
  91 0;
  92
  93 /* go through some gyrations here to verify result of uniteigenvectors */
  94
  95 (kill(S),
  96  [e1, e2, e3, e4]: [foo[2][1][1], foo[2][2][1], foo[2][3][1], foo[2][4][1]],
  97  S: transpose (matrix (e1, e2, e3, e4)),
  98  ratsimp (S . transpose (S)));
  99 matrix ([1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]);
 100
 101 (kill(Lambda),
 102  Lambda: ident(4),
 103  for i thru 4 do Lambda[i, i]: foo[1][1][i],
 104  ratsimp (S . Lambda . transpose (S) - mm));
 105 matrix ([0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]);
 106
 107 /* now test in presence of uv array */
 108
 109 (kill(uv), uv[1]:0);
 110 0;
 111
 112 member ('uv, arrays);
 113 true;
 114
 115 (bar: uniteigenvectors(mm),
 116  is (bar = foo));
 117 true;
 118
 119 kill(uv);
 120 done;