/* Original version of this file copyright 1999 by Michael Wester,
 * and retrieved from http://www.math.unm.edu/~wester/demos/MatrixTheory/problems.macsyma
 * circa 2006-10-23.
 *
 * Released under the terms of the GNU General Public License, version 2,
 * per message dated 2007-06-03 from Michael Wester to Robert Dodier
 * (contained in the file wester-gpl-permission-message.txt).
 *
 * See: "A Critique of the Mathematical Abilities of CA Systems"
 * by Michael Wester, pp 25--60 in
 * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
 * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
 */
/* ----------[ M a c s y m a ]---------- */
/* ---------- Initialization ---------- */
showtime: all$
prederror: false$
/* ---------- Matrix Theory ---------- */
ratmx: true$
keepfloat: true$
/* Extract the superdiagonal => [2, 6] */
mat_diag(matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 1);
/* (2, 3)-minor => [[1, 2], [7, 8]] */
minor(matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 2, 3);
/* Create the 7 x 6 matrix B from rearrangements of the elements of the 4 x 4
   matrix A (this is easiest to do with a MATLAB style notation):
   B = [A(1:3,2:4), A([1,2,4],[3,1,4]); A, [A(1:2,3:4); A([4,1],[3,2])]]
   => [[12 13 14|13 11 14],
       [22 23 24|23 21 24],
       [32 33 34|43 41 44],
       [--------+--+-----]
       [11 12 13 14|13 14],
       [21 22 23 24|23 24],
       [           +-----]
       [31 32 33 34|43 42],
       [41 42 43 44|13 12]].  See Michael James Wester, _Symbolic Calculation
   and Expression Swell Analysis of Matrix Determinants and Eigenstuff_, Ph.D.
   dissertation, University of New Mexico, Albuquerque, New Mexico, December
   1992, p. 89. */
A: matrix([11, 12, 13, 14],
          [21, 22, 23, 24],
          [31, 32, 33, 34],
          [41, 42, 43, 44])$
[< A[1..3,2..4], A[[1,2,4],[3,1,4]]; A, [A[1..2,3..4]; A[[4,1],[3,2]]] >];
remvalue(A)$
/* Create a block diagonal matrix */
mat_diag([matrix([a, 1],[0, a]), b, matrix([c, 1, 0],[0, c, 1],[0, 0, c])]);
/* => [[1 1], [1 0]] */
modulus: 2$
rat(matrix([7, 11], [3, 8]));
modulus: false$
/* => [[-cos t, -sin t], [sin t, -cos t]] */
matrix([cos(t), sin(t)], [-sin(t), cos(t)]);
diff(%, t, 2);
/* => [[(a + 7) x + (2 a - 8) y,   (3 a - 9) x + (4 a + 10) y,
        (5 a + 11) x + (6 a - 12) y]] */
matrix([x, y]) . (a*matrix([1, 3, 5], [2, 4, 6])
                  + matrix([7, -9, 11], [-8, 10, -12]));
/* Matrix norms: infinity norm => 7 */
mat_norm(matrix([1, -2*%i], [-3*%i, 4]), 'inf);
/* Frobenius norm => (a^2 + b^2 + c^2)/(|a| |b| |c|)   (a, b, c real) */
mat_norm(matrix([a/(b*c), 1/c, 1/b], [1/c, b/(a*c), 1/a], [1/b, 1/a, c/(a*b)]),
         'f);
scanmap(factor, %);
/* Hermitian (complex conjugate transpose) => [[1, f(4 + 5 i)], [2 - 3 i, 6]]
   (This assumes f is a real valued function.  In general, the (1, 2) entry
   will be conjugate[f(4 - 5 i)] = conjugate(f)(4 + 5 i).) */
hermitian(matrix([1, 2 + 3*%i], [f(4 - 5*%i), 6]));
m: matrix([a, b], [1, a*b]);
/* Invert the matrix => 1/(a^2 - 1) [[a, -1], [-1/b, a/b]] */
minv: m^^(-1);
m . minv;
remvalue(m, minv)$
/* Inverse of a triangular partitioned (or block) matrix
   => [[A_11^(-1), -A_11^(-1) A_12 A_22^(-1)], [0, A_22^(-1)]].
   See Charles G. Cullen, _Matrices and Linear Transformations_, Second
   Edition, Dover Publications Inc., 1990, p. 35. */
declare([A11, A12, A22], nonscalar)$
matrix([A11, A12], [0, A22])^^(-1);
remove([A11, A12, A22], nonscalar)$
/* LU decomposition of a symbolic matrix   [David Wood]
   [ 1    0   0] [1  x-2  x-3]   [ 1      x-2        x-3    ]
   [x-1   1   0] [0   4   x-5] = [x-1  x^2-3x+6   x^2-3x-2  ]
   [x-2  x-3  1] [0   0   x-7]   [x-2    x^2-8   2x^2-12x+14] */
matrix([ 1,     x-2,         x-3     ],
       [x-1, x^2-3*x+6,   x^2-3*x-2  ],
       [x-2,   x^2-8,   2*x^2-12*x+14])$
ratsimp(lu_decomp_symb(%));
lu_matrices(%, 1..mat_nrows(%));
/* Reduced row echelon form   [Cullen, p. 43]
   => [[1 0 -1 0 2], [0 1 2 0 -1], [0 0 0 1 3], [0 0 0 0 0]] */
matrix([1, 2, 3, 1, 3],
       [3, 2, 1, 1, 7],
       [0, 2, 4, 1, 1],
       [1, 1, 1, 1, 4])$
echelon(%);
/* => 2.  See Gerald L. Bradley, _A Primer of Linear Algebra_, Prentice-Hall,
      Inc., 1975, p. 135. */
rank(matrix([-1, 3, 7, -5], [4, -2, 1, 3], [2, 4, 15, -7]));
/* => 1 */
rank(matrix([2*sqrt(2), 8], [6*sqrt(6), 24*sqrt(3)]));
/* => 1 */
rank(matrix([sin(2*t), cos(2*t)],
            [2*(1 - cos(t)^2)*cos(t), (1 - 2*sin(t)^2)*sin(t)]));
/* Null space => [[2 4 1 0], [0 -3 0 1]]^T or variant   [Bradley, p. 207] */
nullspace(matrix([1, 0, -2, 0], [-2, 1, 0, 3], [-1, 2, -6, 6]));
/* Define a Vandermonde matrix (useful for doing polynomial interpolations) */
matrix([1,   1,   1,   1  ],
       [w,   x,   y,   z  ],
       [w^2, x^2, y^2, z^2],
       [w^3, x^3, y^3, z^3]);
determinant(%);
/* The following formula implies a general result:
   => (w - x) (w - y) (w - z) (x - y) (x - z) (y - z) */
factor(%);
/* Minimum polynomial => (lambda - 1)^2 (lambda + 1)   [Cullen, p. 181] */
matrix([17,  -8, -12, 14],
       [46, -22, -35, 41],
       [-2,   1,   4, -4],
       [ 4,  -2,  -2,  3])$
/* Compute the eigenvalues of a matrix from its characteristic polynomial
   => lambda = {1, -2, 3} */
m: matrix([ 5, -3, -7],
          [-2,  1,  2],
          [ 2, -3, -4]);
charpoly(m, lambda);
solve(% = 0, lambda);
remvalue(m)$
/* In theory, an easy eigenvalue problem! => lambda = {2 - a} for k = 1..100
   [Wester, p. 154] */
ratmx: false$
sparse: true$
eigenvalues((2 - a)*ident(100));
sparse: false$
ratmx: true$
/* => lambda = {4 sin^2(pi k/[2 (n + 1)])} for k = 1..n for an n x n matrix.
      For n = 5, lambda = {2 - sqrt(3), 1, 2, 3, 2 + sqrt(3)}
   See J. H. Wilkinson, _The Algebraic Eigenvalue Problem_, Oxford University
   Press, 1965, p. 307. */
matrix([2, 1, 0, 0, 0],
       [1, 2, 1, 0, 0],
       [0, 1, 2, 1, 0],
       [0, 0, 1, 2, 1],
       [0, 0, 0, 1, 2])$
eigenvalues(%);
/* Eigenvalues of the Rosser matrix.  This matrix is notorious for causing
   numerical eigenvalue routines to fail.   [Wester, p. 146 (Cleve Moler)]
   => {-10 sqrt(10405), 0, 510 - 100 sqrt(26), 1000, 1000,
       510 + 100 sqrt(26), 1020, 10 sqrt(10405)} =
      {-1020.049, 0, 0.098, 1000, 1000, 1019.902, 1020, 1020.049} */
rosser: matrix([ 611,  196, -192,  407,   -8,  -52,  -49,   29],
               [ 196,  899,  113, -192,  -71,  -43,   -8,  -44],
               [-192,  113,  899,  196,   61,   49,    8,   52],
               [ 407, -192,  196,  611,    8,   44,   59,  -23],
               [  -8,  -71,   61,    8,  411, -599,  208,  208],
               [ -52,  -43,   49,   44, -599,  411,  208,  208],
               [ -49,   -8,    8,   59,  208,  208,   99, -911],
               [  29,  -44,   52,  -23,  208,  208, -911,   99])$
eigenvalues(rosser);
errcatch(eigenvalues(sfloat(rosser)));
eigenvalues_by_schur(sfloat(rosser));
remvalue(rosser)$
/* Eigenvalues of the generalized hypercompanion matrix of
   (x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0)*(x^2 + x + 1)^2
   => {[-1 +- sqrt(3) i]/2, [-1 +- sqrt(3) i]/2,
       RootsOf(x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0)} */
matrix([-a4, -a3, -a2, -a1, -a0,  0,  0,  0,  0],
       [  1,   0,   0,   0,   0,  0,  0,  0,  0],
       [  0,   1,   0,   0,   0,  0,  0,  0,  0],
       [  0,   0,   1,   0,   0,  0,  0,  0,  0],
       [  0,   0,   0,   1,   0,  0,  0,  0,  0],
       [  0,   0,   0,   0,   0, -1, -1,  0,  0],
       [  0,   0,   0,   0,   0,  1,  0,  0,  0],
       [  0,   0,   0,   0,   0,  0,  1, -1, -1],
       [  0,   0,   0,   0,   0,  0,  0,  1,  0])$
eigenvalues(%);
/* Eigenvalues and eigenvectors => lambda = {a, a, a, 1 - i, 1 + i},
   eigenvectors = [[1 0 0 0 0], [0 0 1 0 0], [0 0 0 1 0],
                   [0, (1 + i)/2, 0, 0, 1], [0, (1 - i)/2, 0, 0, 1]]^T */
matrix([a,  0, 0, 0, 0],
       [0,  0, 0, 0, 1],
       [0,  0, a, 0, 0],
       [0,  0, 0, a, 0],
       [0, -2, 0, 0, 2]);
eigenvectors(%);
/* Eigenvalues and generalized eigenvectors   [Johnson and Riess, p. 193]
   => lambda = {1, 1, 1}, eigenvectors = [[4 -1 4], [1 -1 2], [3 -1 3]]^T */
matrix([-1,  -8, 1],
       [-1,  -3, 2],
       [-4, -16, 7])$
eigenvectors(%);
/* Eigenvalues and generalized eigenvectors   [Johnson and Riess, p. 199]
   => lambda = {1, 1, 1, 1, 2, 2}, eigenvectors =
      [[1 -1  0  0  0 0], [-1 0 0 1 0 0], [0 0 1 -1 0 -1],
       [0  0 -1 -2 -1 3], [ 0 2 0 0 0 0], [2 0 1  1 0  0]]^T */
matrix([1, 0, 1, 1, 0, 1],
       [1, 2, 0, 0, 0, 0],
       [0, 0, 2, 0, 1, 1],
       [0, 0, 1, 1, 0, 0],
       [0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1, 1])$
eigenvectors(%);
/* Jordan form => diag([[1 1],[0 1]], [[1 1],[0 1]], -1)   [Gantmacher, p. 172]
   */
matrix([1,  0,  0,  1, -1],
       [0,  1, -2,  3, -3],
       [0,  0, -1,  2, -2],
       [1, -1,  1,  0,  1],
       [1, -1,  1, -1,  2])$
jordan_form(%);
/* Smith normal form => [[1, 0], [0, x^4 - x^2 + 1]]   [Cullen, p. 230] */
matrix([x^2, x - 1], [x + 1, x^2]);
/* Matrix exponential => e [[cos 2, -sin 2], [sin 2, cos 2]] */
mat_expm(matrix([1, -2], [2, 1]));
rectform(%);
/* Matrix exponential   [Rick Niles] =>
   [[1, 4 sin(w t)/w - 3 t , 6 [w t - sin(w t)], 2/w [1 - cos(w t)] ],
    [0, 4 cos(w t) - 3     , 6 w [1 - cos(w t)], 2 sin(w t)         ],
    [0, -2/w [1 - cos(w t)], 4 - 3 cos(w t)    , sin(w t)/w         ],
    [0, -2 sin(w t)        , 3 w sin(w t)      , cos(w t)           ]] */
matrix([0, 1,    0,     0  ],
       [0, 0,    0,     2*w],
       [0, 0,    0,     1  ],
       [0, -2*w, 3*w^2, 0  ])$
rectform(matrix_exp(%, t));
/* Sine of a Jordan matrix => diag([[sin a, cos a],[0, sin a]], sin b,
   [[sin c, cos c, -sin(c)/2],[0, sin c, cos c],[0, 0, sin c]])
   See F. R. Gantmacher, _The Theory of Matrices_, Volume One, Chelsea
   Publishing Company, 1977, p. 100 to see how to do a general function. */
matrix([a, 1, 0, 0, 0, 0],
       [0, a, 0, 0, 0, 0],
       [0, 0, b, 0, 0, 0],
       [0, 0, 0, c, 1, 0],
       [0, 0, 0, 0, c, 1],
       [0, 0, 0, 0, 0, c])$
mat_sinm(%);
/* Sine of a matrix => [[1 0 0], [0 1 0], [0 0 1]]   [Cullen, p. 261] */
%pi/2*matrix([2, 1, 1], [2, 3, 2], [1, 1, 2])$
mat_sinm(%);
/* Matrix square root => {+-[[3 1], [1 4]], +-1/sqrt(5) [[-1 7], [7 6]]} */
matrix([10, 7], [7, 17]);
mat_sqrtm(%);
denest_sqrts(%);
/* Square root of a non-singular matrix   [Gantmacher, p. 233]
   => [[e, (e - n) v w + e/2, (n - e) v], [0, e, 0], [0, (e - n) w, n]
   for arbitrary v and w with arbitrary signs e and n = +-1 */
matrix([1, 1, 0], [0, 1, 0], [0, 0, 1])$
mat_sqrtm(%);
% . %;
/* Square root of a singular matrix   [Gantmacher, p. 239]
   => [[0 a b], [0 0 0], [0 1/b 0]] for arbitrary a and b */
matrix([0, 1, 0], [0, 0, 0], [0, 0, 0])$
errcatch(mat_sqrtm(%));
/* Singular value decomposition
   => [1/sqrt(14)  3/sqrt(10) 1/sqrt(35) ] [2 sqrt(7) 0] [1/sqrt(2)  1/sqrt(2)]
      [2/sqrt(14)  0           -sqrt(5/7)] [0         0] [1/sqrt(2) -1/sqrt(2)]
      [3/sqrt(14) -1/sqrt(10) 3/sqrt(35) ] [0         0]
      = U Sigma V^T --- singular values are [2 sqrt(7), 0] */
matrix([1, 1], [2, 2], [3, 3]);
[u, s, vT] <: svd_symb(%);
ratsimp(u . s . vT);
remvalue(u, s, vT)$
/* Jacobian of (r cos t, r sin t) => [[cos t, -r sin t], [sin t, r cos t]] */
jacobian([r*cos(t), r*sin(t)], [r, t]);
/* Hessian of r^2 sin t => [[2 sin t, 2 r cos t], [2 r cos t, -r^2 sin t]] */
r^2*sin(t);
/* Wronskian of (cos t, sin t) => [[cos t, sin t], [-sin t, cos t]] */
wronskian([cos(t), sin(t)], t);
/* How easy is it to define functions to do the last three operations?
   Jacobian of (r cos t, r sin t) => [[cos t, -r sin t], [sin t, r cos t]] */
MYjacobian(f, x):= block([n: length(x)],
   genmatrix(lambda([i, j], diff(f[i], x[j])), n, n))$
MYjacobian([r*cos(t), r*sin(t)], [r, t]);
/* Hessian of r^2 sin t => [[2 sin t, 2 r cos t], [2 r cos t, -r^2 sin t]] */
MYhessian(f, x):= block([n: length(x)],
   genmatrix(lambda([i, j], diff(f, x[i], 1, x[j], 1)), n, n))$
MYhessian(r^2*sin(t), [r, t]);
/* Wronskian of (cos t, sin t) => [[cos t, sin t], [-sin t, cos t]] */
MYwronskian(f, x):= block([n: length(f)],
   genmatrix(lambda([i, j], diff(f[j], x, i-1)), n, n))$
MYwronskian([cos(t), sin(t)], t);