Tue Feb 23 16:54:17 MST 1999 aquarius% macsyma Starting Macsyma math engine with no window system... /opt/local/macsyma_422/macsyma.422 local This is Macsyma 422.0 for Sparc (Solaris 2.x) computers. Copyright (c) 1982 - 1998 Macsyma Inc. All rights reserved. Portions copyright (c) 1982 Massachusetts Institute of Technology. All rights reserved. Type "DESCRIBE(TRADE_SECRET);" to see important legal notices. Type "HELP();" for more information. /aquarius/data2/opt/local/macsyma_422/system/init.lsp being loaded. /aquarius/home/wester/macsyma-init.lsp being loaded. (c1) (c2) /* ----------[ M a c s y m a ]---------- */ /* ---------- Initialization ---------- */ symbol_display_case: lower_case$ Time= 0 msecs (c3) showtime: all$ Time= 0 msecs (c4) prederror: false$ Time= 0 msecs (c5) /* ---------- Matrix Theory ---------- */ ratmx: true$ Time= 0 msecs (c6) keepfloat: true$ Time= 0 msecs (c7) /* Extract the superdiagonal => [2, 6] */ mat_diag(matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 1); /aquarius/data2/opt/local/macsyma_422/matrix/lang.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/supportl.so being loaded. Time= 100 msecs [ 2 ] (d7) [ ] [ 6 ] (c8) /* (2, 3)-minor => [[1, 2], [7, 8]] */ minor(matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]), 2, 3); Time= 0 msecs [ 1 2 ] (d8) [ ] [ 7 8 ] (c9) /* 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])$ Time= 0 msecs (c10) [< A[1..3,2..4], A[[1,2,4],[3,1,4]]; A, [A[1..2,3..4]; A[[4,1],[3,2]]] >]; /aquarius/data2/opt/local/macsyma_422/matrix/support.so being loaded. Time= 160 msecs [ 12 13 14 13 11 14 ] [ ] [ 22 23 24 23 21 24 ] [ ] [ 32 33 34 43 41 44 ] [ ] (d10) [ 11 12 13 14 13 14 ] [ ] [ 21 22 23 24 23 24 ] [ ] [ 31 32 33 34 43 42 ] [ ] [ 41 42 43 44 13 12 ] (c11) remvalue(A)$ Time= 0 msecs (c12) /* Create a block diagonal matrix */ mat_diag([matrix([a, 1],[0, a]), b, matrix([c, 1, 0],[0, c, 1],[0, 0, c])]); Time= 60 msecs [ a 1 0 0 0 0 ] [ ] [ 0 a 0 0 0 0 ] [ ] [ 0 0 b 0 0 0 ] (d12) [ ] [ 0 0 0 c 1 0 ] [ ] [ 0 0 0 0 c 1 ] [ ] [ 0 0 0 0 0 c ] (c13) /* => [[1 1], [1 0]] */ modulus: 2$ Time= 0 msecs (c14) rat(matrix([7, 11], [3, 8])); Time= 0 msecs [ 1 1 ] (d14)/R/ [ ] [ 1 0 ] (c15) modulus: false$ Time= 0 msecs (c16) /* => [[-cos t, -sin t], [sin t, -cos t]] */ matrix([cos(t), sin(t)], [-sin(t), cos(t)]); Time= 0 msecs [ cos(t) sin(t) ] (d16) [ ] [ - sin(t) cos(t) ] (c17) diff(%, t, 2); Time= 10 msecs [ - cos(t) - sin(t) ] (d17) [ ] [ sin(t) - cos(t) ] (c18) /* => [[(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])); Time= 10 msecs (d18)/R/ Col 1 = [ (2 a - 8) y + (a + 7) x ] Col 2 = [ (4 a + 10) y + (3 a - 9) x ] Col 3 = [ (6 a - 12) y + (5 a + 11) x ] (c19) /* Matrix norms: infinity norm => 7 */ mat_norm(matrix([1, -2*%i], [-3*%i, 4]), 'inf); /aquarius/data2/opt/local/macsyma_422/matrix/mathfuns.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/simpfuns.so being loaded. Time= 200 msecs (d19) 7 (c20) /* 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); Time= 390 msecs 2 2 2 c b a 2 2 2 (d20) sqrt(----- + ----- + ----- + -- + -- + --) 2 2 2 2 2 2 2 2 2 a b a c b c c b a (c21) scanmap(factor, %); Time= 340 msecs 2 2 2 c + b + a (d21) -------------------- abs(a) abs(b) abs(c) (c22) /* 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])); Time= 10 msecs [ 1 f(5 %i + 4) ] (d22) [ ] [ 2 - 3 %i 6 ] (c23) m: matrix([a, b], [1, a*b]); Time= 0 msecs [ a b ] (d23) [ ] [ 1 a b ] (c24) /* Invert the matrix => 1/(a^2 - 1) [[a, -1], [-1/b, a/b]] */ minv: m^^(-1); Time= 10 msecs [ a 1 ] [ ------ - ------ ] [ 2 2 ] [ a - 1 a - 1 ] (d24)/R/ [ ] [ 1 a ] [ - ---------- ---------- ] [ 2 2 ] [ (a - 1) b (a - 1) b ] (c25) m . minv; Time= 10 msecs [ 1 0 ] (d25)/R/ [ ] [ 0 1 ] (c26) remvalue(m, minv)$ Time= 0 msecs (c27) /* 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)$ Time= 0 msecs (c28) matrix([A11, A12], [0, A22])^^(-1); Time= 0 msecs [ 1 a12 ] [ --- - ------- ] [ a11 a11 a22 ] (d28)/R/ [ ] [ 1 ] [ 0 --- ] [ a22 ] (c29) remove([A11, A12, A22], nonscalar)$ Time= 0 msecs (c30) /* 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])$ Time= 0 msecs (c31) ratsimp(lu_decomp_symb(%)); /aquarius/data2/opt/local/macsyma_422/matrix/luaux.so being loaded. Time= 100 msecs [ 1 x - 2 x - 3 ] [ ] (d31) [ x - 1 4 x - 5 ] [ ] [ x - 2 x - 3 x - 7 ] (c32) lu_matrices(%, 1..mat_nrows(%)); Time= 10 msecs [ 1 0 0 ] [ 1 x - 2 x - 3 ] [ 1 0 0 ] [ 1 0 0 ] [ ] [ ] [ ] [ ] (d32) [[ x - 1 1 0 ], [ 0 4 x - 5 ], [ 0 1 0 ], [ 0 1 0 ]] [ ] [ ] [ ] [ ] [ x - 2 x - 3 1 ] [ 0 0 x - 7 ] [ 0 0 1 ] [ 0 0 1 ] (c33) /* 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])$ Time= 0 msecs (c34) echelon(%); Time= 0 msecs [ 1 2 3 1 3 ] [ ] [ 0 1 2 0 - 1 ] (d34)/R/ [ ] [ 0 0 0 1 3 ] [ ] [ 0 0 0 0 0 ] (c35) /* => 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])); Time= 0 msecs (d35) 2 (c36) /* => 1 */ rank(matrix([2*sqrt(2), 8], [6*sqrt(6), 24*sqrt(3)])); Time= 0 msecs (d36) 2 (c37) /* => 1 */ rank(matrix([sin(2*t), cos(2*t)], [2*(1 - cos(t)^2)*cos(t), (1 - 2*sin(t)^2)*sin(t)])); Time= 10 msecs (d37) 2 (c38) /* 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])); Time= 20 msecs 1 1 3 1 (d38)/R/ [[1, -, -, -], [0, 1, 0, - -]] 5 2 5 3 (c39) /* 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]); Time= 210 msecs [ 1 1 1 1 ] [ ] [ w x y z ] [ ] (d39) [ 2 2 2 2 ] [ w x y z ] [ ] [ 3 3 3 3 ] [ w x y z ] (c40) determinant(%); Time= 30 msecs 2 2 2 2 2 3 (d40) ((x - w) y + (w - x ) y + w x - w x) z 3 3 3 3 3 2 + ((w - x) y + (x - w ) y - w x + w x) z 2 2 3 3 3 2 2 3 3 2 2 2 3 + ((x - w ) y + (w - x ) y + w x - w x ) z + (w x - w x ) y 3 3 2 3 2 2 3 + (w x - w x) y + (w x - w x ) y (c41) /* The following formula implies a general result: => (w - x) (w - y) (w - z) (x - y) (x - z) (y - z) */ factor(%); Time= 110 msecs (d41) (x - w) (y - w) (y - x) (z - w) (z - x) (z - y) (c42) /* 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])$ Time= 0 msecs (c43) /* 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]); Time= 0 msecs [ 5 - 3 - 7 ] [ ] (d43) [ - 2 1 2 ] [ ] [ 2 - 3 - 4 ] (c44) charpoly(m, lambda); Time= 10 msecs 3 2 (d44) - lambda + 2 lambda + 5 lambda - 6 (c45) solve(% = 0, lambda); Time= 40 msecs (d45) [lambda = 3, lambda = - 2, lambda = 1] (c46) remvalue(m)$ Time= 0 msecs (c47) /* In theory, an easy eigenvalue problem! => lambda = {2 - a} for k = 1..100 [Wester, p. 154] */ ratmx: false$ Time= 0 msecs (c48) sparse: true$ Time= 0 msecs (c49) eigenvalues((2 - a)*ident(100)); /aquarius/data2/opt/local/macsyma_422/matrix/eigen.so being loaded. /aquarius/data2/opt/local/macsyma_422/library1/sprdet.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/meigen.so being loaded. Time= 3370 msecs (d49) [[2 - a], [100]] (c50) sparse: false$ Time= 0 msecs (c51) ratmx: true$ Time= 0 msecs (c52) /* => 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])$ Time= 0 msecs (c53) eigenvalues(%); Time= 90 msecs (d53) [[2, 3, 2 - sqrt(3), sqrt(3) + 2, 1], [1, 1, 1, 1, 1]] (c54) /* 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])$ Time= 0 msecs (c55) eigenvalues(rosser); Time= 430 msecs (d55) [[510 - 100 sqrt(26), 100 sqrt(26) + 510, 1020, - 10 sqrt(10405), 10 sqrt(10405), 1000, 0], [1, 1, 1, 1, 1, 2, 1]] (c56) errcatch(eigenvalues(sfloat(rosser))); Polynomial quotient is not exact Time= 50 msecs (d56) [] (c57) eigenvalues_by_schur(sfloat(rosser)); /aquarius/data2/opt/local/macsyma_422/share/nkaux.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/schur.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/schurl.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/bla_lu.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/blinalgl.so being loaded. Time= 2210 msecs (d57) [- 1020.0490184299985d0, 0, 0, 1020.0490184299975d0, 1020.000000000001d0, 1019.9019513592783d0, 999.99999999999966d0, 1000.0d0] (c58) remvalue(rosser)$ Time= 0 msecs (c59) /* 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])$ Time= 0 msecs (c60) eigenvalues(%); /aquarius/data2/opt/local/macsyma_422/library1/combin.so being loaded. Time= 390 msecs sqrt(3) %i + 1 sqrt(3) %i - 1 (d60) [[- --------------, --------------], [2, 2], 2 2 2 3 4 5 a0 + a1 %lambda + a2 %lambda + a3 %lambda + a4 %lambda + %lambda ] (c61) /* 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]); Time= 0 msecs [ a 0 0 0 0 ] [ ] [ 0 0 0 0 1 ] [ ] (d61) [ 0 0 a 0 0 ] [ ] [ 0 0 0 a 0 ] [ ] [ 0 - 2 0 0 2 ] (c62) eigenvectors(%); Time= 940 msecs (d62) [[a, 3, [[1, 0, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 1, 0, 0]]], %i + 1 %i - 1 [1 - %i, 1, [[0, ------, 0, 0, 1]]], [%i + 1, 1, [[0, - ------, 0, 0, 1]]]] 2 2 (c63) /* 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])$ Time= 0 msecs (c64) eigenvectors(%); Time= 60 msecs 1 (d64) [[1, 3, [[1, - -, 0]]]] 4 (c65) /* 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])$ Time= 0 msecs (c66) eigenvectors(%); Time= 170 msecs (d66) [[2, 2, [[0, 1, 0, 0, 0, 0]]], [1, 4, [[1, - 1, 0, 0, 0, 0]]]] (c67) /* 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])$ Time= 10 msecs (c68) jordan_form(%); /aquarius/data2/opt/local/macsyma_422/matrix/jordan.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/jordanl.so being loaded. /aquarius/data2/opt/local/macsyma_422/share/powers.so being loaded. Time= 370 msecs [ - 1 0 0 0 0 ] [ ] [ 0 1 1 0 0 ] [ ] (d68) [ 0 0 1 0 0 ] [ ] [ 0 0 0 1 1 ] [ ] [ 0 0 0 0 1 ] (c69) /* Smith normal form => [[1, 0], [0, x^4 - x^2 + 1]] [Cullen, p. 230] */ matrix([x^2, x - 1], [x + 1, x^2]); Time= 10 msecs [ 2 ] [ x x - 1 ] (d69) [ ] [ 2 ] [ x + 1 x ] (c70) /* Matrix exponential => e [[cos 2, -sin 2], [sin 2, cos 2]] */ mat_expm(matrix([1, -2], [2, 1])); /aquarius/data2/opt/local/macsyma_422/matrix/phase1.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/matexp.so being loaded. /aquarius/data2/opt/local/macsyma_422/library1/result.so being loaded. /aquarius/data2/opt/local/macsyma_422/library1/laplac.so being loaded. Time= 860 msecs [ 2 %i + 1 1 - 2 %i 2 %i + 1 1 - 2 %i ] [ %e %e %i %e %e %i ] [ ---------- + ---------- ------------- - ------------- ] [ 2 2 2 2 ] (d70) [ ] [ 1 - 2 %i 2 %i + 1 2 %i + 1 1 - 2 %i ] [ %e %i %i %e %e %e ] [ ------------- - ------------- ---------- + ---------- ] [ 2 2 2 2 ] (c71) rectform(%); Time= 30 msecs [ %e cos(2) - %e sin(2) ] (d71) [ ] [ %e sin(2) %e cos(2) ] (c72) /* 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 ])$ Time= 0 msecs (c73) rectform(matrix_exp(%, t)); /aquarius/data2/opt/local/macsyma_422/library1/binoml.so being loaded. Time= 1690 msecs [ 4 sin(t w) 2 2 cos(t w) ] [ 1 ---------- - 3 t 6 t w - 6 sin(t w) - - ---------- ] [ w w w ] [ ] [ 0 4 cos(t w) - 3 6 w - 6 w cos(t w) 2 sin(t w) ] (d73) [ ] [ 2 cos(t w) 2 sin(t w) ] [ 0 ---------- - - 4 - 3 cos(t w) -------- ] [ w w w ] [ ] [ 0 - 2 sin(t w) 3 w sin(t w) cos(t w) ] (c74) /* 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])$ Time= 0 msecs (c75) mat_sinm(%); /aquarius/data2/opt/local/macsyma_422/matrix/na.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/matsolve.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/struct.so being loaded. /aquarius/data2/opt/local/macsyma_422/matrix/linalg.so being loaded. Time= 3870 msecs [ sin(a) cos(a) 0 0 0 0 ] [ ] [ 0 sin(a) 0 0 0 0 ] [ ] [ 0 0 sin(b) 0 0 0 ] [ ] (d75)/R/ [ sin(c) ] [ 0 0 0 sin(c) cos(c) - ------ ] [ 2 ] [ ] [ 0 0 0 0 sin(c) cos(c) ] [ ] [ 0 0 0 0 0 sin(c) ] (c76) /* 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])$ Time= 10 msecs (c77) mat_sinm(%); Time= 750 msecs [ 1 0 0 ] [ ] (d77)/R/ [ 0 1 0 ] [ ] [ 0 0 1 ] (c78) /* Matrix square root => {+-[[3 1], [1 4]], +-1/sqrt(5) [[-1 7], [7 6]]} */ matrix([10, 7], [7, 17]); Time= 0 msecs [ 10 7 ] (d78) [ ] [ 7 17 ] (c79) mat_sqrtm(%); Time= 520 msecs (d79)/R/ matrix([((sqrt(7 sqrt(5) + 27) + sqrt(27 - 7 sqrt(5))) sqrt(5) - sqrt(7 sqrt(5) + 27) + sqrt(27 - 7 sqrt(5)))/(2 sqrt(2) sqrt(5)), 2 ((sqrt(7 sqrt(5) + 27) - sqrt(27 - 7 sqrt(5))) sqrt(5) - sqrt(7 sqrt(5) + 27) + sqrt(27 - 7 sqrt(5)))/(4 sqrt(2) sqrt(5))], sqrt(7 sqrt(5) + 27) - sqrt(27 - 7 sqrt(5)) [-------------------------------------------, sqrt(2) sqrt(5) ((sqrt(7 sqrt(5) + 27) + sqrt(27 - 7 sqrt(5))) sqrt(5) + sqrt(7 sqrt(5) + 27) - sqrt(27 - 7 sqrt(5)))/(2 sqrt(2) sqrt(5))]) (c80) denest_sqrts(%); /aquarius/data2/opt/local/macsyma_422/share/sqdnst.so being loaded. Time= 820 msecs [ 3 1 ] (d80) [ ] [ 1 4 ] (c81) /* 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])$ Time= 0 msecs (c82) mat_sqrtm(%); Time= 250 msecs [ %r14 %r16 %r14 ] [ 1 - ----------- ------ ] [ 2 %r18 %r15 2 %r18 ] [ ] [ %r16 - 2 %r18 %r15 ] (d82)/R/ [ 0 - ------------- ------ ] [ 2 %r18 2 %r18 ] [ ] [ 2 ] [ %r16 %r16 + 2 %r18 ] [ 0 - ----------- ------------- ] [ 2 %r18 %r15 2 %r18 ] (c83) % . %; Time= 60 msecs [ %r14 %r16 %r14 ] [ 1 - --------- ---- ] [ %r15 %r18 %r18 ] [ ] [ %r18 - %r16 %r15 ] (d83)/R/ [ 0 ----------- ---- ] [ %r18 %r18 ] [ ] [ 2 ] [ %r16 %r18 + %r16 ] [ 0 - --------- ----------- ] [ %r15 %r18 %r18 ] (c84) /* 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])$ Time= 0 msecs (c85) errcatch(mat_sqrtm(%)); Division by 0 Time= 110 msecs (d85) [] (c86) /* 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]); Time= 0 msecs [ 1 1 ] [ ] (d86) [ 2 2 ] [ ] [ 3 3 ] (c87) [u, s, vT] <: svd_symb(%); Time= 790 msecs [ 1 3 sqrt(5) ] [ --------------- -------- --------- ] [ sqrt(7) sqrt(2) sqrt(10) 5 sqrt(7) ] [ ] [ 2 sqrt(7) 0 ] [ 2 sqrt(5) ] [ ] (d87)/R/ [[ --------------- 0 - ------- ], [ 0 0 ], [ sqrt(7) sqrt(2) sqrt(7) ] [ ] [ ] [ 0 0 ] [ 3 1 3 sqrt(5) ] [ --------------- - -------- --------- ] [ sqrt(7) sqrt(2) sqrt(10) 5 sqrt(7) ] [ 1 1 ] [ ------- ------- ] [ sqrt(2) sqrt(2) ] [ ]] [ 1 1 ] [ ------- - ------- ] [ sqrt(2) sqrt(2) ] (c88) ratsimp(u . s . vT); Time= 90 msecs [ 1 1 ] [ ] (d88) [ 2 2 ] [ ] [ 3 3 ] (c89) remvalue(u, s, vT)$ Time= 0 msecs (c90) /* 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]); /aquarius/data2/opt/local/macsyma_422/matrix/matfuncs.so being loaded. Time= 120 msecs [ cos(t) - r sin(t) ] (d90) [ ] [ sin(t) r cos(t) ] (c91) /* 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); Time= 0 msecs 2 (d91) r sin(t) (c92) /* Wronskian of (cos t, sin t) => [[cos t, sin t], [-sin t, cos t]] */ wronskian([cos(t), sin(t)], t); Time= 0 msecs [ cos(t) sin(t) ] (d92) [ ] [ - sin(t) cos(t) ] (c93) /* 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))$ Time= 0 msecs (c94) MYjacobian([r*cos(t), r*sin(t)], [r, t]); Time= 10 msecs [ cos(t) - r sin(t) ] (d94) [ ] [ sin(t) r cos(t) ] (c95) /* 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))$ Time= 0 msecs (c96) MYhessian(r^2*sin(t), [r, t]); Time= 10 msecs [ 2 sin(t) 2 r cos(t) ] (d96) [ ] [ 2 ] [ 2 r cos(t) - r sin(t) ] (c97) /* 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))$ Time= 0 msecs (c98) MYwronskian([cos(t), sin(t)], t); Time= 0 msecs [ cos(t) sin(t) ] (d98) [ ] [ - sin(t) cos(t) ] (c99) /* ---------- Quit ---------- */ quit(); Bye. real 29.64 user 22.13 sys 2.52