# file maple5

# Hill, Section 5.1, Determinants.
with(linalg):
A:= matrix([[a11,a12],
            [a21,a22]]):
det(A);                 # See Hill Section 5.1 Definition 5.3, p316
# 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

# Expansion of a 3 by 3 determinant along its third row
with(linalg):
A:= matrix(3,3):
A31:= submatrix(A,[1,2],[2,3]):
A32:= submatrix(A,[1,2],[1,3]):
A33:= submatrix(A,[1,2],[1,2]):
simplify( A[3,1]*det(A31)-A[3,2]*det(A32)+A[3,3]*det(A33) - det(A) );
quit

# Hill Section 5.1 Theorem 5.7 p321
with(linalg):
# If A is an n by n lower triangular matrix, then
# det(A) is the product of the diagonal entries of A.  For example
L:= matrix( [[a11,  0,  0],
             [a21,a22,  0],
             [a31,a32,a33]]):
det(L);
#
# Hill Section 2.3 Theorem 2.25 p106
# Let A be an n by n matrix.    Then det(transpose(A)) = det(A).  For example
#
A:= matrix(3,3):
simplify(det(transpose(A))-det(A));
#
# COROLLARY Any true statement about det(A) is still true if the words `row'
# and `column' are interchanged wherever they occur.  Thus theorem 5.7 yields
U:= matrix( [[a11,a12,a13],
             [  0,a22,a23],
             [  0,  0,a33]]):
det(U);
#
# THEOREM. det(A)=0 if and only if the rows of A are linearly dependent
#    (This includes theorems 2.7, 2.12 of Hill 2.2).    For example
A0:= matrix( [[a11,a12],
              [a11,a12]]):
det(A0);
#
# Hill Section 5.1 Theorem 5.9 p323
# Let A and B be n by n matrices.    Then det( evalm(A&*B) ) = det(A)*det(B).
#
B:= matrix(3,3):
simplify( det(evalm(A&*B)) - det(A)*det(B) );
quit

# Hill Section 5.2 Examples 4 and 5, pp. 328-330
with(linalg):
A:= matrix([ [5,-7,7],
             [4,-3,4],
             [4,-1,2] ]):
# The eigenvalues are the roots of the characteristic polynomial:
charpoly(A,lambda);                         # This is  det( lambda * Id - A)
solve(",lambda);
# Alternatively, use the function  eigenvals  in the  linalg package:
eigenvals(A);
B:= matrix([ [0,-1],
             [1, 0] ]):
eigenvals(B);

# Continuing Hill Section 5.2 Example 4 (p328):
Id:= diag(1,1,1):                        W5:= nullspace( evalm(A - 5*Id) );
W1:= nullspace( evalm(A - Id) );
W_2:= nullspace( evalm(A + 2*Id) );
# For example,  W1  is a basis for  NS(A-Id),  and so is a set of eigenvectors.
# To CHECK that the vectors of this set satisfy the equation  A x = x,  find
evalm(A &* W1[1]);  # The product of A with the first (and only) vector in W1
# You should CHECK the other eigenvectors for yourself.
quit

# Hill Section 5.3 Examples 1 and 2, p338.
with(linalg):
A:= matrix([ [2,3,0],
             [4,3,0],
             [0,0,6] ]):
evects:= eigenvects(A);
# The eigenvalues are  evects[1][1] and evects[2][1],  occurring
#   evects[1][2] and evects[2][2]  times respectively, and the corresponding
#   eigenvectors are the elements of the sets  evects[1][3] and evects[2][3].
#
# Since one of the eigenvalues is repeated, we force it to be the second one
if evects[1][2] = 2 then    evects:= evects[2],evects[1]    fi;
S:= augment( evects[1][3][1],evects[2][3][1],evects[2][3][2] );
evalsA:= evects[1][1], evects[2][1], evects[2][1]:
equal( evalm( inverse(S)&*A&*S ), diag(evalsA) );

# Hill Section 5.3 Example 3, p340.
B:= matrix([ [a,1], [0,a] ]):
evalsB:= eigenvals(B);
nullspace( B-a*diag(1,1) );
# Since the dimension of this eigenspace is only 1,   B is not diagonalizable

# An example of the power of a matrix (see Hill Section 5.3 Theorem 5.34)
#
# The k'th power of the diagonal form of the matrix A in example 1 (p338) is
diag( evalsA[1]^k,evalsA[2]^k,evalsA[3]^k );
# and the  k'th  power of A is
evalm( S &* " &* inverse(S) );
quit

with(linalg):
u:= vector([I,1+I]):                            v:= vector([-2*I,1+I]):
simplify(dotprod(u,v)),                         innerprod(u,v);
quit

# Hill Section 5.4 Example 3, p352.
with(linalg):
# Find an orthogonal matrix Q which diagonalizes the matrix A where
A:= matrix([ [0, 2, 2],
             [2, 0,-2],
             [2,-2, 0] ]):
# Use the same ideas as in Example 1 of Hill Section 5.3:
evects:= eigenvects(A);
# One of the eigenvalues is repeated:  ensure that it is the second one
if evects[1][2] = 2 then    evects:= evects[2],evects[1]    fi;
#
u1:= evects[1][3][1]:       u2:= evects[2][3][1]:      u3:= evects[2][3][2]:
print( dotprod(u1,u2), dotprod(u1,u3), dotprod(u2,u3) );
v1:= evalm( u1/norm(u1,2) ):
# Since  u2, u3  are not orthogonal, use Gram-Schmidt:
GS:= GramSchmidt( [u2,u3] ):
v2:= evalm( GS[1]/norm(GS[1],2) ):    v3:= evalm( GS[2]/norm(GS[2],2) ):
Q:= augment( v1,v2,v3 );
orthog(Q);
evals:= evects[1][1], evects[2][1], evects[2][1]:
equal( evalm( transpose(Q)&*A&*Q ), diag(evals) );
quit

# Consider the matrix Q produced by the housetrans procedure of Hill 4.6
with(linalg):
Q:= matrix([ [ -2/3,  -1/3,   2/3 ],
             [ -1/3, 14/15,  2/15 ],
             [  2/3,  2/15, 11/15 ] ]):
orthog(Q), det(Q);
# So  Q  reflects any vector  u  orthogonal to the plane P where
u:= nullspace( Q+diag(1,1,1) )[1]:
P:= dotprod( u,vector([x,y,z]) ) = 0;
# It also rotates  P  through an angle  theta  where
theta:= arccos( (trace(Q)+1)/2 );
# So (as expected)  Q  is simply a reflection in the plane P.
quit

# Hill Section 1.1 Example 4, p4  (and pp 356-358)
with(linalg):
T:= matrix([ [8/10, 1/10],
             [2/10, 9/10] ]):
eigenvects(T);
# The eigenvector associated with the eigenvalue  1  yields
evalm( 90000000 * "[1][3][1] );
# which agrees with the long-term population given in Chapter 1.
quit

# Hill Section 5.6 Example 3, p366.
with(linalg):
A:= matrix([ [5,-7,7],
             [4,-3,4],
             [4,-1,2] ]):
evalm( exponential(t*A) &* vector([c1-c3,c1+c2,c1+c2+c3]) );
quit

# Hill Section 5.7 Example 1, p369:    Graphs of some conic sections
#
# The following is the circle    4*x^2 + 4*y^2 = 1
plot( [cos(t)/2, sin(t)/2, t=-Pi..Pi] );
# The following is the ellipse    9*x^2 + 4*y^2 = 1
plot( [cos(t)/3, sin(t)/2, t=-Pi..Pi] );
# The following is the right hand-side of the hyperbola    9*x^2 - 4*y^2 = 1
plot( [sec(t)/3, tan(t)/2, t=-Pi/4..Pi/4] );
# A more familiar example is the following (rotated) hyperbola
plot( [x, 1/x, x=-infinity..infinity] );
quit

# Moving the origin to the centre  (h,k)  of the quadratic locus  Q2  where
Q2:= a11*X1^2 + a22*X2^2 + 2*a12*X1*X2 + 2*b1*X1 + 2*b2*X2 - c:
q2:= expand( subs(X1=x1-h,X2=x2-k, Q2) ):
#
# We want the coefficients of x1,x2 to be 0: collect q2 as a polynomial in x1
qx1:= collect( q2,x1 ):                        # Also as a polynomial in x2:
qx2:= collect( q2,x2 ):
coeff1:= coeff(qx1,x1);
coeff2:= coeff(qx2,x2);
# These include the unwanted term  2 a12 x1 x2:  extract the constant term
lprint();   solve( {coeff(coeff1,x2,0)=0,coeff(coeff2,x1,0)=0}, {h,k} );
# CHECK that these values give a quadratic in  x1,x2  where the coefficients
#    of x1,x2 are 0.    (This is not as easy as it seems!)
quit

with(linalg):
# Hill Section 5.7 Example 5, p374.
# a)
A:= matrix([ [  3 ,-2/3,-2/3],
             [-2/3, 8/3,  0 ],
             [-2/3,  0 ,10/3] ]):
# is the symmetric matrix for the quadratic form
lprint();    X:= vector( [x,y,z] ):   qxyz:= evalm( transpose(X) &* A &* X );

# b) Diagonalize A   (See Examples 1,3 of Hill Sections 5.3, 5.4 respectively)
evects:= eigenvects(A);
u1:= evalm(evects[1][3][1]/norm(evects[1][3][1],2)):
u2:= evalm(evects[2][3][1]/norm(evects[2][3][1],2)):
u3:= evalm(evects[3][3][1]/norm(evects[3][3][1],2)):
S:= augment( u1,u2,u3 );
evals:= evects[1][1], evects[2][1], evects[3][1]:
B:= evalm( transpose(S)&*A&*S ):                      equal( B,diag(evals) );
# Thus if the surface  qxyz = 12  is rotated so that the x1-, y1-, z1-axes
#  point along  u1, u2, u3  respectively, it becomes the ellipsoid
X:= vector( [x1,y1,z1] ):               evalm( transpose(X) &* B &* X ) = 12;
quit

# Examples from Hill Section 5.8
with(linalg):

# First load the qrfac2  procedure from chapter 4
housetrans:= proc(v)
    local Id,a,w;
    Id:= diag( 1$vectdim(v) ):          a:= - sign( evalf(v[1]) ):
    w:= copy(v):                        w[1]:= v[1] - a*norm(v,2):
    map( simplify, evalm( Id - 2*w&*transpose(w)/dotprod(w,w) ) ):
end:
qrfac2:= proc(A,Q)
    local m,n,R,j,cj,Qj:
    m:= rowdim(A):                    n:= coldim(A):
    Q:= housetrans( col(A,1) );
    R:= evalm( Q&*A );
    for j from 2 to n do
        cj:= subvector(R,j..m,j); # The last m-j+1 entries of the j'th column
        Qj:= diag( 1$(j-1), housetrans(cj) );
        Q:= map( simplify, evalm(Q&*Qj) ); R:= map( simplify, evalm(Qj&*R) );
    od;
end:

# The following matrix is used by Hill in all the examples of Section 5.8
A:= matrix([ [4,-2],  [3,-1] ]);
EV:= eigenvects(A);
# Make sure the eigenvectors are in the same order as Hill lists them:
if EV[1][1] = 1 then      EV:= EV[2],EV[1]      fi;
v1:= map( evalf, evalm(EV[1][3][1]/norm(EV[1][3][1],2)) );
v2:= map( evalf, evalm(EV[2][3][1]/norm(EV[2][3][1],2)) );
# Example 1 (p378)  Direct iteration
direct_it:= proc(A,w,m,lambda) #  w = initial value,  m = no. of iterations.
    local j,wj,Aw:
    wj:= w:
    for j from 1 to m do
        Aw:= map( evalf, evalm(A&*wj) ):
        lambda:= evalf( norm(Aw,2) ):
        wj:= map( evalf, evalm(Aw/lambda) ):
        print( eval(lambda),wj );  # we need  eval  because lambda is a name
    od:
    wj:
end:
direct_it( A,vector([1,0]),5,'lambda' ):
# this is converging (slowly) to the eigenvalue 2  and eigenvector v1.
#
# Example 3 (p381)  Inverse iteration
inverse_it:= proc(A,w,m,lambda)  #  w,m  as for direct iteration.
    local wj,y,j:
    wj:= w:
    for j from 1 to m do
        y:= linsolve( A,wj ):
        lambda:= evalf( 1/norm(y,2) ):
        wj:= map( evalf, evalm(lambda*y) ):
        print( 1/eval(lambda),wj );      # Hill lists the INVERSE of lambda
    od:
    wj:
end:
inverse_it( A,vector([1,0]),4,'lambda' ):
# lambda,y  are converging to the eigenvalue 1  and eigenvector v2.
#
# Example 4 (p382)  Shifted inverse iteration - this is faster than the above.
#   First find a reasonable approximation to use as the shift
w:= direct_it( A,vector([1,0]),3,'a' ):
inverse_it( A-diag(a,a),w,3,'lambda' ):
print(a-lambda);       # This is an approximation to the eigenvalue  2
# Example 5 (p383/384)  The QR Algorithm
#  -  the diagonal entries converge to the eigenvalues
Ak:= copy(A):
for k from 1 to 4 do
    Rk:= map( evalf, qrfac2(Ak,'Qk') ):         #See Section 4.6 for  qrfac2
    Ak:= map( evalf, evalm(Rk&*Qk) ):                 print(Ak);
od:
# The Shifted QR Algorithm  (starting with the above results)
n:= 2:                      a:= Ak[n,n]:
for k from 1 to 3 do
    Rk:= map( evalf, qrfac2(Ak-diag(a,a),'Qk') ):
    Ak:= map( evalf, evalm(Rk&*Qk+diag(a,a)) ):                 print(Ak);
od:
quit

