Mon Feb 16 11:31:13 MST 1998 aquarius% mupad *----* MuPAD 1.4.0 -- Multi Processing Algebra Data Tool /| /| *----* | Copyright (c) 1997 - 98 by SciFace Software GmbH | *--|-* All rights reserved. |/ |/ *----* Licensed to: Michael Wester >> # ----------[ M u P A D ]---------- # >> # ---------- Initialization ---------- # >> TEXTWIDTH:= 80: >> read("../../Time.mupad"): >> # ---------- Matrix Theory ---------- # >> export(Dom): export(linalg): >> matrix:= Matrix(ExpressionField(normal)): >> # Extract the superdiagonal => [2, 6] # >> matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]): >> [%[j - 1, j] $ j = 2..ncols(%)]; [2, 6] Time: 4220 msec Type: DOM_LIST >> # (2, 3)-minor => [[1, 2], [7, 8]] # >> delCol(delRow(matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), 2), 3); +- -+ | 1, 2 | | | | 7, 8 | +- -+ Time: 3810 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # 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]]): >> stackMatrix(A[1..3, 2..4] . extractMatrix(A,[1,2,4],[3,1,4]), &> A . stackMatrix(A[1..2, 3..4], extractMatrix(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 | +- -+ Time: 790 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> A:= NIL: >> # Create a block diagonal matrix # >> #diag([matrix([[a, 1],[0, a]]), b, matrix([[c, 1, 0],[0, c, 1],[0, 0, c]])]);# >> # => [[1 1], [1 0]] # >> matrix([[7, 11], [3, 8]]) mod 2; Error: Domain attribute "_mod" missing >> map(matrix([[7, 11], [3, 8]]), modp, 2); Error: Domain attribute "modp" missing >> # => [[-cos t, -sin t], [sin t, -cos t]] # >> matrix([[cos(t), sin(t)], [-sin(t), cos(t)]]); +- -+ | cos(t), sin(t) | | | | -sin(t), cos(t) | +- -+ Time: 1230 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> diff(%, t$2); +- -+ | -cos(t), -sin(t) | | | | sin(t), -cos(t) | +- -+ Time: 220 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # => [[(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]])); array(1..1, 1..3, (1, 1) = 7 x - 8 y + a x + 2 a y, (1, 2) = - 9 x + 10 y + 3 a x + 4 a y , (1, 3) = 11 x - 12 y + 5 a x + 6 a y ) Time: 770 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Matrix norms: infinity norm => 7 # >> norm(matrix([[1, -2*I], [-3*I, 4]])); 7 Time: 360 msec Type: Dom::ExpressionField(normal, iszero@normal) >> # Frobenius norm => (a^2 + b^2 + c^2)/(|a| |b| |c|) (a, b, c real) # >> assume(a, Type::RealNum): >> assume(b, Type::RealNum): >> assume(c, Type::RealNum): >> norm(matrix([[a/(b*c), 1/c, 1/b], [1/c, b/(a*c), 1/a], [1/b, 1/a, c/(a*b)]]), &> Frobenius); / / 2 2 2 2 2 2 | | 2 abs(a) abs(b) + 2 abs(a) abs(c) + 2 abs(b) abs(c) + | \ \ 2 2 2 / a \2 2 2 2 / b \2 abs(a) abs(b) abs(c) abs| --- | + abs(a) abs(b) abs(c) abs| --- | + \ b c / \ a c / 2 2 2 / c \2 \ 2 2 2 \ abs(a) abs(b) abs(c) abs| --- | | / (abs(a) abs(b) abs(c) ) |^(1/2) \ a b / / | / Time: 5690 msec Type: Dom::ExpressionField(normal, iszero@normal) >> simplify(%); / / 2 2 2 2 2 2 | | 2 abs(a) abs(b) + 2 abs(a) abs(c) + 2 abs(b) abs(c) + | \ \ 2 2 2 / a \2 2 2 2 / b \2 abs(a) abs(b) abs(c) abs| --- | + abs(a) abs(b) abs(c) abs| --- | + \ b c / \ a c / 2 2 2 / c \2 \ 2 2 2 \ abs(a) abs(b) abs(c) abs| --- | | / (abs(a) abs(b) abs(c) ) |^(1/2) \ a b / / | / Time: 99200 msec Type: Dom::ExpressionField(normal, iszero@normal) >> unassume(a, Type::RealNum): >> unassume(b, Type::RealNum): >> unassume(c, Type::RealNum): >> # 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).) # >> conjugate(transpose(matrix([[1, 2 + 3*I], [f(4 - 5*I), 6]]))); +- -+ | 1, conjugate(f(4 - 5 I)) | | | | 2 - 3 I, 6 | +- -+ Time: 1300 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> m:= matrix([[a, b], [1, a*b]]); +- -+ | a, b | | | | 1, a b | +- -+ Time: 320 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Invert the matrix => 1/(a^2 - 1) [[a, -1], [-1/b, a/b]] # >> minv:= 1/m; +- -+ | a 1 | | ------, - ------ | | 2 2 | | a - 1 a - 1 | | | | 1 a | | - -----------, ----------- | | 2 2 | | - b + a b - b + a b | +- -+ Time: 1080 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> m * minv; +- -+ | 1, 0 | | | | 0, 1 | +- -+ Time: 580 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> m:= NIL: >> minv:= NIL: >> # 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. # >> matrix([[A11, A12], [0, A22]])^(-1); +- -+ | 1 A12 | | ---, - ------- | | A11 A11 A22 | | | | 1 | | 0, --- | | A22 | +- -+ Time: 1060 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # 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]]): >> # 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]]): >> gaussJordan(%); +- -+ | 1, 0, -1, 0, 2 | | | | 0, 1, 2, 0, -1 | | | | 0, 0, 0, 1, 3 | | | | 0, 0, 0, 0, 0 | +- -+ Time: 1650 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # => 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]])); 2 Time: 590 msec Type: DOM_INT >> # => 1 # >> rank(matrix([[2*sqrt(2), 8], [6*sqrt(6), 24*sqrt(3)]])); 2 Time: 610 msec Type: DOM_INT >> # => 1 # >> rank(matrix([[sin(2*t), cos(2*t)], &> [2*(1 - cos(t)^2)*cos(t), (1 - 2*sin(t)^2)*sin(t)]])); 2 Time: 770 msec Type: DOM_INT >> # 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]])); -- +- -+ +- -+ -- | | 2 | | 0 | | | | | | | | | | 4 | | -3 | | | | |, | | | | | 1 | | 0 | | | | | | | | | | 0 | | 1 | | -- +- -+ +- -+ -- Time: 640 msec Type: DOM_LIST >> # 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]]); +- -+ | 1, 1, 1, 1 | | | | w, x, y, z | | | | 2 2 2 2 | | w , x , y , z | | | | 3 3 3 3 | | w , x , y , z | +- -+ Time: 470 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> det(%); 2 3 3 2 2 3 2 3 3 2 3 2 2 3 3 2 - w x y + w x y + w x y - w x y - w x y + w x y + w x z - w x z 2 3 2 3 3 2 3 2 2 3 3 2 2 3 - w x z + w x z + w x z - w x z - w y z + w y z + w y z - 2 3 3 2 3 2 2 3 3 2 2 3 2 3 3 2 w y z - w y z + w y z + x y z - x y z - x y z + x y z + x y z 3 2 - x y z Time: 3180 msec Type: Dom::ExpressionField(normal, iszero@normal) >> # The following formula implies a general result: &> => (w - x) (w - y) (w - z) (x - y) (x - z) (y - z) # >> Factor(%); (w - x) (- w + y ) (x - y) (- w + z ) (x - z) (y - z) Time: 950 msec Type: Dom::ExpressionField(normal, iszero@normal) >> # 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]]); +- -+ | 5, -3, -7 | | | | -2, 1, 2 | | | | 2, -3, -4 | +- -+ Time: 760 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> charPolynomial(m, lambda); 3 2 poly(lambda + (-2) lambda + (-5) lambda + 6, [lambda], Dom::ExpressionField(normal, iszero@normal)) Time: 1990 msec Type: DOM_POLY >> solve(% = 0, lambda); {1, -2, 3} Time: 1380 msec Type: DOM_SET >> m:= NIL: >> # In theory, an easy eigenvalue problem! => lambda = {2 - a} for k = 1..100 &> [Wester, p. 154] # >> eigenValues((2 - a)*matrix(100, 100, fun(1), Diagonal), All); Terminated real 6113.84 user 6091.87 sys 11.03 ------------------------------------------------------------------------------- Mon Feb 23 16:08:09 MST 1998 aquarius% mupad *----* MuPAD 1.4.0 -- Multi Processing Algebra Data Tool /| /| *----* | Copyright (c) 1997 - 98 by SciFace Software GmbH | *--|-* All rights reserved. |/ |/ *----* Licensed to: Michael Wester >> # ----------[ M u P A D ]---------- # >> # ---------- Initialization ---------- # >> TEXTWIDTH:= 80: >> read("../../Time.mupad"): >> # ---------- Matrix Theory ---------- # >> export(Dom): export(linalg): >> matrix:= Matrix(ExpressionField(normal)): >> # => 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(%, All); 1/2 1/2 {[1, 1], [2, 1], [3 + 2, 1], [3, 1], [- 3 + 2 , 1]} Time: 12760 msec Type: Dom::Multiset >> # 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, All); 1/2 1/2 {[1020, 1], [1000, 2], [- 100 26 + 510 , 1], [100 26 + 510, 1], 1/2 1/2 [- 10 10405 , 1], [10 10405 , 1], [0, 1]} Time: 14580 msec Type: Dom::Multiset >> eigenValues(float(rosser), All); {[-1020.049018, 1], [-0.00000000009232777908, 1], [0.09804864081, 1], [999.9989606, 1], [1000.001039, 1], [1019.90053, 1], [1020.004551, 1], [1020.045887, 1]} Time: 253500 msec Type: Dom::Multiset >> numeric::eigenvalues(float(rosser)); Error: non-numeric entries not handled [numeric::eigenvalues] >> rosser:= NIL: >> # 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(%, All); 1/2 1/2 {[- 1/2 I 3 - 1/2 , 2], [1/2 I 3 - 1/2, 2], [RootOf(a0 + 5 2 a1 ____eigenValue_x + ____eigenValue_x + a2 ____eigenValue_x + 3 4 a3 ____eigenValue_x + a4 ____eigenValue_x , ____eigenValue_x), 1]} Time: 44180 msec Type: Dom::Multiset >> # 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]]); +- -+ | 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: 810 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> eigenVectors(%); -- -- -- +- -+ +- -+ +- -+ -- -- | | | | 1 | | 0 | | 0 | | | | | | | | | | | | | | | | | | 0 | | 0 | | 0 | | | | | | | | | | | | | | | | a, 3, | | 0 |, | 1 |, | 0 | | |, | | | | | | | | | | | | | | | 0 | | 0 | | 1 | | | | | | | | | | | | | | | | | | 0 | | 0 | | 0 | | | -- -- -- +- -+ +- -+ +- -+ -- -- -- -- +- -+ -- -- | | | 0 | | | | | | | | | | | | 1/2 - 1/2 I | | | | | | | | | | 1 + I, 1, | | 0 | | |, | | | | | | | | | 0 | | | | | | | | | | | | 1 | | | -- -- +- -+ -- -- -- -- +- -+ -- -- -- | | | 0 | | | | | | | | | | | | | | 1/2 + 1/2 I | | | | | | | | | | | | 1 - I, 1, | | 0 | | | | | | | | | | | | | | 0 | | | | | | | | | | | | | | 1 | | | | -- -- +- -+ -- -- -- Time: 6400 msec Type: DOM_LIST >> # 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(%); -- -- -- +- -+ -- -- -- | | | | -4 | | | | | | | | | | | | | | 1, 3, | | 1 | | | | | | | | | | | | | | | | 0 | | | | -- -- -- +- -+ -- -- -- Time: 1770 msec Type: DOM_LIST >> # 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(%); -- -- -- +- -+ -- -- -- -- +- -+ -- -- -- | | | | 0 | | | | | | -1 | | | | | | | | | | | | | | | | | | | | | | 1 | | | | | | 1 | | | | | | | | | | | | | | | | | | | | | | 0 | | | | | | 0 | | | | | | 2, 2, | | | | |, | 1, 4, | | | | | | | | | | 0 | | | | | | 0 | | | | | | | | | | | | | | | | | | | | | | 0 | | | | | | 0 | | | | | | | | | | | | | | | | | | | | | | 0 | | | | | | 0 | | | | -- -- -- +- -+ -- -- -- -- +- -+ -- -- -- Time: 5420 msec Type: DOM_LIST >> # 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]]): >> jordanForm(%); +- -+ | -1, 0, 0, 0, 0 | | | | 0, 1, 0, 0, 0 | | | | 0, 1, 1, 0, 0 | | | | 0, 0, 0, 1, 0 | | | | 0, 0, 0, 1, 1 | +- -+ Time: 8490 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Smith normal form => [[1, 0], [0, x^4 - x^2 + 1]] [Cullen, p. 230] # >> matrix([[x^2, x - 1], [x + 1, x^2]]); +- -+ | 2 | | x , x - 1 | | | | 2 | | x + 1, x | +- -+ Time: 610 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Matrix exponential => e [[cos 2, -sin 2], [sin 2, cos 2]] # >> exp(matrix([[1, -2], [2, 1]])); array(1..2, 1..2, exp(1 - 2 I) exp(1 + 2 I) (1, 1) = ------------ + ------------, 2 2 (1, 2) = - 1/2 I exp(1 - 2 I) + 1/2 I exp(1 + 2 I) , (2, 1) = 1/2 I exp(1 - 2 I) - 1/2 I exp(1 + 2 I), exp(1 - 2 I) exp(1 + 2 I) (2, 2) = ------------ + ------------ 2 2 ) Time: 1850 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> rectform(%); array(1..2, 1..2, exp(1 - 2 I) exp(1 + 2 I) (1, 1) = ------------ + ------------, 2 2 (1, 2) = - 1/2 I exp(1 - 2 I) + 1/2 I exp(1 + 2 I) , (2, 1) = 1/2 I exp(1 - 2 I) - 1/2 I exp(1 + 2 I), exp(1 - 2 I) exp(1 + 2 I) (2, 2) = ------------ + ------------ 2 2 ) Time: 810 msec Type: 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 ]]): >> exp(%, t); array(1..4, 1..4, (1, 1) = 1, - 3 t w + 2 I exp(-I t w) - 2 I exp(I t w) (1, 2) = -------------------------------------------, w (1, 3) = 6 t w - 3 I exp(-I t w) + 3 I exp(I t w), - exp(-I t w) - exp(I t w) + 2 (1, 4) = -------------------------------, w (2, 1) = 0, (2, 2) = 2 exp(-I t w) + 2 exp(I t w) - 3, (2, 3) = 6 w - 3 w exp(-I t w) - 3 w exp(I t w), (2, 4) = I exp(-I t w) - I exp(I t w), (3, 1) = 0, exp(-I t w) + exp(I t w) - 2 (3, 2) = ----------------------------, w 3 exp(-I t w) 3 exp(I t w) (3, 3) = - ------------- - ------------ + 4 , 2 2 I exp(-I t w) - I exp(I t w) (3, 4) = ----------------------------, 2 w (4, 1) = 0, (4, 2) = - I exp(-I t w) + I exp(I t w) , (4, 3) = 3/2 I w exp(-I t w) - 3/2 I w exp(I t w), exp(-I t w) exp(I t w) (4, 4) = ----------- + ---------- 2 2 ) Time: 8540 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # 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. # >> m:= 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]]): >> sin(%); / +- -+ \ | | a, 1, 0, 0, 0, 0 | | | | | | | | 0, a, 0, 0, 0, 0 | | | | | | | | 0, 0, b, 0, 0, 0 | | sin| | | | | | 0, 0, 0, c, 1, 0 | | | | | | | | 0, 0, 0, 0, c, 1 | | | | | | | | 0, 0, 0, 0, 0, c | | \ +- -+ / Time: 1380 msec Type: "sin" >> Im(exp(m*I)); Im(FAIL) Time: 31540 msec Type: "Im" >> # Sine of a matrix => [[1 0 0], [0 1 0], [0 0 1]] [Cullen, p. 261] # >> m:= PI/2*matrix([[2, 1, 1], [2, 3, 2], [1, 1, 2]]): >> sin(%); / +- -+ \ | | PI PI | | | | PI, --, -- | | | | 2 2 | | | | | | | | 3 PI | | sin| | PI, ----, PI | | | | 2 | | | | | | | | PI PI | | | | --, --, PI | | | | 2 2 | | \ +- -+ / Time: 1410 msec Type: "sin" >> Im(exp(m*I)); Im(FAIL) Time: 109040 msec Type: "Im" >> m:= NIL: >> # Matrix square root => {+-[[3 1], [1 4]], +-1/sqrt(5) [[-1 7], [7 6]]} # >> matrix([[10, 7], [7, 17]]); +- -+ | 10, 7 | | | | 7, 17 | +- -+ Time: 1370 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> sqrt(%); FAIL Time: 680 msec Type: DOM_FAIL >> # 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]]): >> sqrt(%); FAIL Time: 1380 msec Type: DOM_FAIL >> # 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]]): >> sqrt(%); FAIL Time: 1400 msec Type: DOM_FAIL >> # 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]]); +- -+ | 1, 1 | | | | 2, 2 | | | | 3, 3 | +- -+ Time: 720 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> #u * s * vT; &> u:= NIL: &> s:= NIL: &> vT:= NIL:# >> # 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]); +- -+ | cos(t), -r sin(t) | | | | sin(t), r cos(t) | +- -+ Time: 1240 msec Type: Dom::Matrix(Dom::ExpressionField(id, iszero)) >> # Hessian of r^2 sin t => [[2 sin t, 2 r cos t], [2 r cos t, -r^2 sin t]] # >> hessian(r^2*sin(t), [r, t]); +- -+ | 2 sin(t), 2 r cos(t) | | | | 2 | | 2 r cos(t), - r sin(t) | +- -+ Time: 800 msec Type: Dom::Matrix(Dom::ExpressionField(id, iszero)) >> # Wronskian of (cos t, sin t) => [[cos t, sin t], [-sin t, cos t]] # >> [cos(t), sin(t)]; [cos(t), sin(t)] Time: 680 msec Type: DOM_LIST >> # 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:= proc(f, x) local n; &> begin &> n:= nops(x); &> matrix([[diff(f[i], x[j]) $ j = 1..n] $ i = 1..n]) &> end_proc: >> MYjacobian([r*cos(t), r*sin(t)], [r, t]); +- -+ | cos(t), -r sin(t) | | | | sin(t), r cos(t) | +- -+ Time: 1470 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Hessian of r^2 sin t => [[2 sin t, 2 r cos t], [2 r cos t, -r^2 sin t]] # >> MYhessian:= proc(f, x) local n; &> begin &> n:= nops(x); &> matrix([[diff(f, x[i], x[j]) $ j = 1..n] $ i = 1..n]) &> end_proc: >> MYhessian(r^2*sin(t), [r, t]); +- -+ | 2 sin(t), 2 r cos(t) | | | | 2 | | 2 r cos(t), - r sin(t) | +- -+ Time: 1460 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # Wronskian of (cos t, sin t) => [[cos t, sin t], [-sin t, cos t]] # >> MYwronskian:= proc(f, x) local n; &> begin &> n:= nops(f); &> matrix([[diff(f[j], x$(i - 1)) $ j = 1..n] $ i = 1..n]) &> end_proc: >> MYwronskian([cos(t), sin(t)], t); +- -+ | cos(t), sin(t) | | | | -sin(t), cos(t) | +- -+ Time: 1470 msec Type: Dom::Matrix(Dom::ExpressionField(normal, iszero@normal)) >> # ---------- Quit ---------- # >> quit real 521.01 user 518.65 sys 1.45