# file maple2

# Hill, Section 2.1, Example 8, p94
with(linalg):
A:= matrix([ [a11,a12],
             [a21,a22] ]):
det(A);
# The formula for the determinant of a 3 by 3 matrix can also be displayed.
A:= matrix([ [a11,a12,a13],
             [a21,a22,a23],
             [a31,a32,a33] ]):
det(A);  # This is left to the reader to run.
quit

# Hill, Section 2.2, Evaluating Determinants.
with(linalg):
# Instead of checking some of the numerical examples in Hill, we illustrate
# the results of some of his theorems using general 2 by 2 matrices.
#
# We follow our earlier style to deal with matrices with general entries:
A:=matrix([ [a11,a12],
            [a21,a22] ]):
#
# THEOREM 2.7 (p97). If a square matrix A has a row of all zeros then det(A)=0.
#
# Example
A0:=matrix([ [a11,a12],
             [0  ,  0] ]):
det(A0);
# THEOREM 2.8 (p98). Let A be an n by n (upper or lower) triangular matrix.
#   Then det(A) is the product of the diagonal entries of A.
#
# Example
U:=matrix([ [a11,a12,a13],
            [  0,a22,a23],
            [  0,  0,a33] ]):
det(U);
# In particular, determinants of elementary matrices are as in Hill Theorem 2.9
#
# THEOREM 2.10 (p99). If E is an elementary n by n matrix and
#   A is any n by n matrix,  then det(evalm(E&*A)) = det(E)*det(A).
#
# COROLLARY 2.11 (p100). Let A be an n by n matrix.
# (a) If B is obtained from A by adding a multiple of one row of A
#     to another row, then  det(B) = det(A).
# (b) If B is obtained from A by multiplying a row of A
#     by the number m, then  det(B) = m*det(A).
# (c) If B is obtained from A by interchanging two rows of A,
#     then  det(B) = -det(A).
#
# Examples
B:=addrow( A,1,2,m):                        simplify( det(B)-det(A) );
B:=mulrow(A, 1, m):                         simplify( det(B)-m*det(A) );
B:=swaprow( A,1,2 ):                        simplify( det(B)+det(A) ); 
# Several applications of Corollary 2.11(b) yield
# THEOREM 2.27 (p108). If A is an n by n matrix, then det(k*A) = k^n*det(A).
#
simplify( det(k*A) - k^2*det(A) );

# Hill, Section 2.3, Properties of the Determinant
# We continue with the general 2 by 2 matrices of Hill, Section 2.2
#
# THEOREM 2.13 (p104). A square matrix A is invertible
#   if and only if det(A) <> 0.
#
B:= subs( a22=a12*a21/a11, op(A) ):                  inverse(B);
########################################################################
# As mentioned earlier, the computation of det involves numerical dangers
# and so (as mentioned on page 105 of Hill) the above theorem is never
# used to test a matrix for invertibility.
########################################################################
#
# THEOREM 2.20 (p105). If A and B are n by n matrices,
#   then  det(evalm(A&*B)) = det(A)*det(B).
#
B:=matrix([ [b11,b12],
            [b21,b22] ]):
simplify( det(evalm(A&*B)) - det(A)*det(B) );
# THEOREM 2.25 (p106). If A is an n by n matrix, det(transpose(A)) = det(A).
#
simplify( det(transpose(A))-det(A) );
quit

# The expansion of a 3 by 3 determinant along its third row, and its adjoint
with(linalg):
A:=matrix([ [a11,a12,a13],
            [a21,a22,a23],
            [a31,a32,a33] ]):
A31:= submatrix(A,[1,2],[2,3]):        A32:= submatrix(A,[1,2],[1,3]):
A33:= submatrix(A,[1,2],[1,2]):
simplify( a31*det(A31)-a32*det(A32)+a33*det(A33) - det(A) );
adjA:= adjoint(A);
# This is the transpose of the matrix of cofactors.  For example
adjA[1,3] - det(A31);
# THEOREM 2.33 (p116). If A is an n by n matrix, then  A adj(A) = det(A) Id.
# For example,
equal( evalm( A&*adjA), det(A)*diag(1,1,1) );
quit

# The simplest use of determinants is an easily-remembered formula for the
#    cross product of two vectors:
with(linalg):
a:= vector( [a1,a2,a3] ):                       b:= vector( [b1,b2,b3] ):
basis3:= vector( [i,j,k] ):                     ab:= augment(a,b):
aXb:= crossprod (a,b):                          dotprod(aXb,basis3);
simplify( " - det(augment(ab,basis3)) );

# A similar formula can be used to find the volume of a tetrahedron:
c:= vector( [c1,c2,c3] ):
with(geom3d):
point3d(A,a1,a2,a3),      point3d(B,b1,b2,b3),      point3d(C,c1,c2,c3):
point3d(Q,0,0,0):           #  Note that O  is a reserved word in Maple
tetrahedron(QABC, [Q, A, B, C]):
simplify( volume(QABC) + det(augment(ab,c))/6 );

# Maple has a  jacobian  command which can be used to evaluate integrals
#    using, for example, polar co-ordinates:
jacobian( [r*cos(theta),r*sin(theta)], [r,theta] );

# and Maple V (from Release 2 on) has a  Wronskian  command:
Wr:= Wronskian( vector([exp(x),exp(-x)]), x );
# If you only have a copy of Maple V or earlier,
# a wronskian command can easily be written:
wronskian:= proc(f,x)    #  f  should be a vector of functions of  x
    local Wr,Wi,i;
    Wi:= f;        Wr:= f;
    for i from 3 to nops(f) do
        Wi:= map(diff,Wi,x);        Wr:= augment(Wr,Wi);
    od:
    evalm(Wr);
end:
Wr:= wronskian( vector([exp(x),exp(-x)]), x );
# Hill (and many other authors) use `Wronskian' for the determinant of this:
simplify( det(Wr) );
quit

