# FILE anton2.mpl
interface(echo=2);

# Section 2.1, Example 7 - AR6 p71; AR7 p83
with(linalg,matrix,det):
A:= matrix([ [a11,a12],
               [a21,a22] ]):
det(A);
# The formula for the det 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.


# Anton, Section 2.2, Evaluating Determinants.
with(linalg,det,matrix,mulrow,swaprow,addrow,
              transpose,multiply,inverse,add):
# Instead of merely checking some of the numerical examples in the text,
# we illustrate the results of some of 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 (AR 2.2.1). If A is any square matrix that contains 
# a row of zeros,   then det(A)=0.
#
# Example
A0:=matrix([ [a11,a12],
               [0  ,  0] ]):
det(A0);

# THEOREM (AR 2.2.2). If A is an n by n triangular matrix, then
# det(A) is the product of the entries on the main diagonal of A.
#
# Example
U:=matrix([ [a11,a12,a13],
            [  0,a22,a23],
            [  0,  0,a33] ]):
det(U);

# THEOREM (AR 2.2.3). Let A be any n by n matrix.
# (a) If B is obtained from A by multiplying a row of A
#     by the number k, then  det(B) = k*det(A).
# (b) If B is obtained from A by interchanging two rows of A,
#     then  det(B) = -det(A).
# (c) If B is obtained from A by adding a multiple of one row of A
#     to another row, then  det(B) = det(A).
#
# Examples
B:=mulrow(A, 1, k):                         simplify( det(B)-k*det(A) );

B:=swaprow( A,1,2 ):                        simplify( det(B)+det(A) ); 

B:=addrow( A,1,2,k):                        simplify( det(B)-det(A) );

A3p:=matrix([[0, 1,5],
             [3,-6,9],
             [2, 6,1]]);
det(A3p);
# Answer is 165


# AR, Section 2.3, Properties of the Determinant
# We continue with general 2 by 2 matrices
#
# THEOREM (AR6 2.3.1; AR7 2.2.1b). Let A be an n by n matrix.
# Then det(transpose(A))-det(A)=0.
#
simplify(det(transpose(A))-det(A));

# THEOREM (AR6 2.3.3; AR7 2.3.4). Let A and B be n by n matrices.
# Then det(multiply(A,B))-(det(A))*(det(B))=0.
#
B:=matrix([[b11,b12],
             [b21,b22]]);
simplify(det(multiply(A,B))-(det(A))*(det(B)));

#
# THEOREM (AR6 2.3.4; AR7 2.3.3.). A square matrix A is invertible
#   if and only if det(A) <> 0.
#
B:= subs( a22=a12*a21/a11, op(A) ):                  inverse(B);
# Error, (in inverse) singular matrix
#
########################################################################
# As mentioned earlier, the computation of det involves numerical dangers
# and so the above theorem is never used to test 
# a matrix of floats for invertibility.
########################################################################

# An example of the expansion of a 3 by 3 determinant along its third row
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) );

with(linalg,det):
# A general conic section through 5 points

prepare_row:= proc(xj,yj):
 [xj^2,xj*yj,yj^2,xj,yj,1]
end:

make_array:= proc(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5):
array(1..6,1..6,
 [prepare_row(x,y),
 prepare_row(x1,y1),
 prepare_row(x2,y2),
 prepare_row(x3,y3),
 prepare_row(x4,y4),
 prepare_row(x5,y5)]):
 end:

# AR S11.1 Example 3 has an orbit of an asteroid example.
# The numbers are nicer in AR S11.1 Q3.
# The equation is lhs=0 where lhs is the following expression
lhEq:=det(make_array(0,0,0,-1,2,0,2,-5,4,-1));
#
# by looking at lhs/160 one sees the quadratic terms are a perfect
square.
# This means that the equation lhs=0 represents a parabola.
expand(lhEq/160-(x+y)^2+(x+y)/2+3*(x-y)/2);

# 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):
A:= 'A': point3d(A,a1,a2,a3);      
B:= 'B': point3d(B,b1,b2,b3);      
C:= 'C': point3d(C,c1,c2,c3):
O:= 'O': point3d(O,0,0,0):                        
tetrahedron(OABC, [O, A, B, C]):
simplify( volume(OABC) + det(augment(ab,c))/6 );

# Maple has a  jacobian  command which can be used to evaluate integrals
# when these require the use, for example, of polar co-ordinates:
jacobian( [r*cos(theta),r*sin(theta)], [r,theta] );
interface(echo=1);
