# file maple4

# See Hill Section 4.1 Figure 4.2 p237.
with(linalg):
# Instead of using the race car given by Hill, we will use the following
surface:= [x*sin(x)*cos(y), x*cos(x)*cos(y), x*sin(y) ]:
domain:= x=0..Pi, y=0..Pi:
plot3d( surface,domain );
Q:= matrix([ [cos(Pi/4),sin(Pi/4),0], [-sin(Pi/4),cos(Pi/4),0], [0,0,1] ]):
rotated:= convert( evalm(Q&*surface), list):   # To make it possible to plot
plot3d( rotated,domain );
Q:= matrix([ [-1,0,0], [0,-1,0], [0,0,1] ]):  # Rotation of Pi in the xy-plane
rotated:= convert( evalm(Q&*surface), list):
plot3d( rotated,domain );
Q1:= matrix([ [1,0,0], [0,cos(Pi/3),sin(Pi/3)], [0,-sin(Pi/3),cos(Pi/3)] ]):
Q2:= matrix([ [cos(Pi/3),0,sin(Pi/3)], [0,1,0], [-sin(Pi/3),0,cos(Pi/3)] ]):
rotated:= convert(evalm(Q1&*Q2&*surface),list): #Rotations in xz and yz-planes
plot3d( rotated,domain );
Q:= matrix([ [1,0,0], [0,-1,0], [0,0,1] ]):   # Reflection through xz-plane
reflected:= convert( evalm(Q&*surface), list):
plot3d( reflected,domain );
quit

with(linalg):
# Hill Section 4.1 Theorem 4.3 p241.
# Show that the induced function of  A  is a linear transformation where
A:= matrix( [[a11,a12,a13],  
             [a21,a22,a23]] ):
u:= vector( [u1,u2,u3] ):                         v:= vector( [v1,v2,v3] ):
evalm( A&*(u+v) - (A&*u+A&*v) ); 
evalm( A&*(r*u) - r*A&*u ); 

# Hill Section 4.1 Theorem 4.4 p242.
# Consider a linear transformation  T  from  R^3  to  R^2
# Ensure that  T  is a linear transformation:
define( Linear(T) );            # This is the Maple V Release 3 version
define( linear(T) );            # This is the version prior to Release 3
# Write
T([1,0,0]):= vector( [a11,a21] ):         T([0,1,0]):= vector( [a12,a22] ):
T([0,0,1]):= vector( [a13,a23] ):
# Show that the matrix A above is the matrix induced by  T.
evalm( T(Pi*[1,0,0]+E*[0,1,0]+gamma*[0,0,1]) - A&*vector([Pi,E,gamma]) );
# Pi, E and gamma must be used since define/Linear only works with constants
#
# Hill Section 4.1 Example 8 p244.
# Show that differentiation is a linear transformation
diff(f(x),x) + diff(g(x),x) - diff(f(x)+g(x),x);
diff(r*f(x),x) - r*diff(f(x),x);
quit

# Hill Section 4.2 Theorem 4.7 p247.
# Show that the following is an inner product:
dot:= proc(fx,gx,lims)        int( expand(fx*gx),x=lims )        end:
combine( dot(f(x),h(x),a..b) +dot(g(x),h(x),a..b) -dot((f+g)(x),h(x),a..b) );
combine( dot(f(x),g(x),a..b) - dot(g(x),f(x),a..b) );
combine( dot(r*f(x),g(x),a..b) - r*dot(f(x),g(x),a..b) );
# combine  was needed above to force the linearity of integration.
#
# An example of this inner product was given in Hill chapter 3:
f:= x->x:                               # Hill Section 3.3 Example 10 (p157)
dot(f(x),sin(x),0..2*Pi);
# The norm of  sin,  and the projection of  f  onto  sin  are given by
norm_sin:= sqrt( dot(sin(x),sin(x),0..2*Pi) );
f_proj_sin:= (dot(f(x),sin(x),0..2*Pi)/dot(sin(x),sin(x),0..2*Pi)) * sin;
quit

# Hill Section 4.2 Theorems 4.15 and 4.16 pp 250,251.
with(linalg):
u:= vector( [u1,u2] ):                             v:= vector( [v1,v2] ):
dotprod(u,v) - evalm( transpose(u)&*v );
# Note that Maple treats vectors as COLUMN vectors (i.e. m by 1 matrices). So
evalm( u&*transpose(v) );
A:= matrix(2,2):
dotprod( evalm(A&*u),v ) - dotprod( u,evalm(transpose(A)&*v) ):
simplify(");
quit

# Hill Section 4.3 Example 6, p265
with(linalg):
# Find the equation  y = m x + b  for the line of best fit for the points
#   (40,482), (45,467), (50,452), (55,433), (60,421)
A:= matrix([ [40,1], [45,1], [50,1], [55,1], [60,1] ]):
y:= vector( [482, 467, 452, 433, 421] ):
xbar:= leastsqrs(A,y):                               map( evalf,xbar );
# So the line of best fit is  y = -3.12 x + 607
#
# The same problem can be done with functions from the stats package.
# This has some advantages in connection with built-in graphics.
with(stats):
#
# If you are using Maple V Release 3,
# remove the (first)  #  from the next 6 lines:
# dat:= [ [40,45,50,55,60], [482,467,452,433,421] ]:
# map( evalf, fit[leastsquare[[x,y],y=a+b*x, {a,b}]](dat) );
# # Other sorts of least-squares fits are possible.
# # Hill S4.3 Q24 p269 has an example of fitting a quadratic curve
# dat:= [ [1,2,3,4], [-1,0,-1,-1] ]:
# map( evalf, fit[leastsquare[[x,y],y=a0+a1*x+a2*x^2 ]](dat) );
#
# If you are using Maple V Release 1 or 2,
# remove the (first)  #  from the next 10 lines:
# dat:= array([ [y,x], [482,40], [467,45], [452,50], [433,55], [421,60] ]):
# v:= regression( dat,y=m*x+b );
# assign(v);              # This assigns the above values of b and m to them
# statplot( dat,y=m*x+b ):
# # Other sorts of least-square fits are possible with the regression function.
# # Hill S4.3 Q24 p269 has an example of fitting a quadratic curve
# dat:= array([ [y,x],[-1,1],[0,2],[-1,3],[-1,4] ]):
# v:= regression( dat, y = a0 + a1*x + a2*x^2 );
# #  Run the above for answers, below for a plot.
# assign(v);                               statplot( dat, y = a0+a1*x+a2*x^2 ):
#

# Continuing Example 6 of Hill Section 4.3
Atr:= transpose(A):                            p:= evalm( A&*xbar):
equal( convert(evalm(Atr &* p),matrix), convert(evalm(Atr &* y),matrix) );
# verifying Hill's  normal equations  4.36.
#
#  p  is the projection of  y  onto  CS(A).
dotprod( y-p, col(A,1) ),                     dotprod( y-p, col(A,2) );
# showing that  y - p  is orthogonal to both the columns of A - and so to CS(A)
#
# The matrix of the perpendicular projection function is
P:= evalm( A &* inverse(Atr&*A) &* Atr );
quit

# Hill Section 4.4 Example 7, p276
with(linalg):
# Find an orthonormal basis for the space spanned by
v1:= vector([1,1,1,1]):    v2:= vector([0,1,1,1]):    v3:= vector([0,0,1,1]):
# This can be done using the function  GramSchmidt  in the  linalg  package:
GramSchmidt( [v1,v2,v3] );
# Warning: Although these vectors are orthogonal, they are not orthonormal,
# and so each should be divided by its norm to give an orthonormal basis
#
# To help understand the process, we also do it using projections:
w1:= v1:                                    u1:= evalm( w1/norm(w1,2) ):
w2:= evalm( v2 - dotprod(u1,v2)*u1 ):       u2:= evalm( w2/norm(w2,2) ):
w3:= evalm( v3 - dotprod(u1,v3)*u1 - dotprod(u2,v3)*u2 ):
u3:= evalm( w3/norm(w3,2) ):
print( u1,u2,u3 );
# Apart from the norms of the vectors, this agrees with the GramSchmidt result.
quit

# Hill Section 4.5 Example 2 p282.
with(linalg):
# Show that the matrix Q is orthogonal, where
Q:= matrix([ [cos(theta),sin(theta)],
             [-sin(theta),cos(theta)] ]):
orthog(Q);
# Alternatively,
equal( map( simplify,evalm(transpose(Q)&*Q) ), diag(1,1) );
# verifying that  transpose(Q)  is  inverse(Q)
quit

# A procedure to find the QR decomposition of a matrix using  GramSchmidt
with(linalg):
qrfac:= proc(A,Q)              # See Hill Theorem 4.64, pp 283,284
    local n,Qcols,R,i;
    n:= coldim(A);
    Qcols:= GramSchmidt( [col(A,1..n)] );
    Q:= transpose( matrix(Qcols) );
#   Normalise Q by dividing each  col(Q,i)  by its norm.
    for i from 1 to n do
        Q:= mulcol( Q,i,1/norm(col(Q,i),2) )
    od;
#   Since  Q &* R = A  and  transpose(Q) &* Q = Id,
    R:= evalm(transpose(Q)&*A);
end:
# Hill Section 4.5 Example 4, p284
# Find the QR decomposition of the matrix
A:= matrix([ [1,0,0],
             [1,1,0],
             [1,1,1],
             [1,1,1] ]):
R:= qrfac( A,'Q'):                                 print ( Q,R);
equal( A, evalm(Q&*R) );

# Hill Section 4.5 Example 5, p285
# Use QR decomposition to find the line  y = mx + b  of best fit for the points
#   (2,3), (3,4), (4,4), (7,7)
#
# We want the least-squares solution for   A c = y   where c = [b,m] and
A:= matrix([ [1,2], [1,3], [1,4], [1,7] ]):
y:= vector( [3, 4, 4, 7] ):
R:= qrfac(A,'Q'):            # See the function  qrfac  given above
c:= linsolve( R, evalm(transpose(Q)&*y) ):
#
# The line of best fit is
y = c[2]*x + c[1];
quit

# Examples on Householder transformations and QR-decomposition based on them
with(linalg):
housetrans:= proc(v)       # See Hill Theorem 4.75, p292
    local Id,a,w;
    Id:= diag( 1$vectdim(v) ):
    a:= - sign( evalf(v[1]) ):
#     (earlier versions of Maple had trouble with signs of square roots)
    w:= copy(v):
    w[1]:= v[1] - a*norm(v,2):
    map( simplify, evalm( Id - 2*w&*transpose(w)/dotprod(w,w) ) ):
end:
# Hill 4.6, Example 2(a), p292/293
Q:= housetrans( vector([2,1,-2]) );

qrfac2:= proc(A,Q)               # See Hill pp 294,295 - returns the matrix R
    local m,n,R,j,cj,Qj:
    m:= rowdim(A):                    n:= coldim(A):
#   Assume that m < n for what follows
    Q:= housetrans( col(A,1) );
    R:= evalm( Q&*A );         # R is used in place of Hill's Aj
    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) );        # update Q
        R:= map( simplify, evalm(Qj&*R) );        # update R
    od;
end:
# Hill Section 4.5 Example 4, p284, (using Householder transformations)
A:= matrix([ [1,0,0],
             [1,1,0],
             [1,1,1],
             [1,1,1] ]):
R:= qrfac2( A,'Q'):
#
# Because the terms get messy, we look at approximations to  Q,R:
map( evalf,Q );
map( evalf,R );
map( evalf,evalm(Q&*R) );
# which is the same as  A  (to within the order of accuracy).
quit

# Hill Section 4.7 Example 1, p298 and Example 5, p303
with(linalg):
# Consider the linear transformation  T  from R^2 to R^2 where
T:= proc(v) vector([ 2*v[1]+3*v[2], 4*v[1]+3*v[2] ]) end: # v[1]:=x: v[2]:=y: 
#
# Example 1:  Find the matrix of  T  relative to the (standard) basis  B  where
e1:= vector([1,0]):      e2:= vector([0,1]):     B:= augment( e1,e2 ):
#
#   This can be done using the coordinates of  T(e1), T(e2)  as in Section 3.8:
T_BB:= augment( linsolve(B,T(e1)), linsolve(B,T(e2)) );

# Now find the matrix of  T  relative to the basis  C  where
v1:= vector([1,-1]):      v2:= vector([3,4]):     C:= augment( v1,v2 ):
T_CC:= augment( linsolve(C,T(v1)), linsolve(C,T(v2)) );
# Example 5:    The transition matrix from  C  to  B  is given by
gaussjord( augment(B,C) ):
P:= submatrix( ", 1..2, 3..4 );
equal( evalm(P^(-1) &* T_BB &* P), T_CC );
quit

# Hill Section 4.7 Example 6, p304
with(linalg):
# Consider the linear transformation  T  from P_2 to P_2 where
T:= proc(f)      local a,b,c;
    a:= coeff(f,x,0);    b:= coeff(f,x,1);    c:= coeff(f,x,2);
    (5*a+b+c) + (b-a-c)*x + (a+b)*x^2;
end:
# The matrix of  T  relative to the standard basis is
T_S:= matrix([ [ 5, 1, 1],
               [-1, 1,-1],
               [ 1, 1, 0] ]):
# If
f1:= 1-x:        f2:= 1+x:       f3:= 1+x+x^2:      C:= [ f1,f2,f3 ]:
# then the transition matrix from  C  to the standard basis is given by
P:= matrix([ [ 1, 1, 1],
             [-1, 1, 1],
             [ 0, 0, 1] ]):
#
# and the matrix of  T  relative to the basis  C  is
T_C:= evalm( P^(-1) &* T_S &* P );
#   This can be checked using the coordinates of
T(f1),`    `, T(f2),`    `, T(f3);
simplify( T_C[1,1]*f1+T_C[2,1]*f2+T_C[3,1]*f3 - T(f1) ),
simplify( T_C[1,2]*f1+T_C[2,2]*f2+T_C[3,2]*f3 - T(f2) ),
simplify( T_C[1,3]*f1+T_C[2,3]*f2+T_C[3,3]*f3 - T(f3) );
quit

