#% FILE anton3toEV.mpl

interface(echo=2);
# dotprod(nVec1,nVec2);
# crossprod(threeVec1,threeVec2);

# AR6 Section 4.3 Example 8 p168; AR7 Section 5.2 Example 9
with(linalg,vector,augment,gausselim,backsub,col):
v1:= vector( [1,2,-1] ):                    v2:= vector( [6,4,2] ):
b:= vector( [b1,b2,b3] ):
# b is a linear combination of v1,v2 if there are scalars a1, a2 with b equal to
evalm(a1*v1 + a2*v2); 

# This can be done as follows
A:= augment(v1,v2):                         x:= vector( [a1,a2] ):
evalm(A&*x);

# Since this is the same as  a1*v1 + a2*v2,  solve the equation    A*x = b
U:= gausselim( augment(A,b),2 );

# where  (a)  b:= vector( [9,2,7] ):
subs (b1=9, b2=2, b3=7, evalm(U) );

# The row of zeroes shows that [9,2,7] is in span {v1,v2}.
# i.e. [9,2,7] = a1*v1+a2*v2, where a1, a2 are given by 
# with(linalg,col): 
x:= backsub(");

#        (b)  b:= vector( [4,-1,8] ):
subs (b1=4, b2=-1, b3=8, evalm(U) );

# The last row shows that [4,-1,8] is not in span {v1,v2}.
# i.e. [4,-1,8] is not a linear combination of v1,v2.



# Example 4 - AR6 Section 4.4 p174; AR7 Section 5.3 p223
with(linalg):
v1:= vector([1,-2,3]):    v2:= vector([5,6,-1]):    v3:= vector([3,2,1]):
# These are linearly dependent if there are numbers c1, c2, c3, not all zero,
# such that                      0 = c1*v1 + c2*v2 + c3*v3
# These can be found as follows
A:= augment( v1,v2,v3 ):
U:= gausselim( augment(A,vector([0,0,0])) );

x:=backsub(");

# This shows that the  v1,v2,v3  above are linearly dependent. For example
evalm( v1 + v2 - 2*v3 );

# In particular    number of pivot columns in U  <=  rowdim(U)  =  rowdim(A)
# So the equations must always have a solution for any vectors  v1, v2, v3


# Example 3 - AR6 Section 4.5 p184; AR7 Section 5.4 p246
with(linalg):
# Show that the vectors v1,v2,v3 form a basis for R^3, where
v1:=vector([1,2,1]):      v2:=vector([2,9,0]):      v3:=vector([3,3,4]):
#
# We need to show that v1, v2, v3 span R^3 and are linearly independent.
# To show they span R^3, consider an arbitrary vector in R^3
b:=vector([b1,b2,b2]):
# and solve the equation    Ac = b  where
A:= augment( v1,v2,v3 );

gausselim( augment(A,b) );

# Since the pivots on the diagonal of this matrix are all non-zero, the
#    equation    Ac = b  has a solution for any b1,b2,b2, and so they span R^3.
# In fact the solutions are unique and so, setting b1 = b2 = b3 = 0,
#    we can see that the equation    Ac = 0  has only one solution.
# Thus v1, v2, v3 are linearly independent.


with(linalg):
I3:= diag(1,1,1):                                  nullspace(I3);

# Output like  {}  indicates that the answer is the zero subspace.
#
# AR6 Section 4.6 Example 3 p195; AR7 Section 5.5 Example 4 p263
# Find the null space of the matrix
A:= matrix( [[ 2, 2,-1, 0, 1],  
               [-1,-1, 2,-3, 1],  
               [ 1, 1,-2, 0,-1],
               [ 0, 0, 1, 1, 1]] ):
nA:= nullspace(A);

# The displayed vectors are a basis for nullSpace(A), or kernel, of A.
#
# Alternatively, solve the equation  Ax = 0  by reducing A to Echelon form
U:= gausselim(A):                                nU:= nullspace(U);
basis([op(nA),op(nU)]);
# Thus nullSpace(U) = nullSpace(A), illustrating the theorem in [AR]

with(linalg):
A:= matrix( [[a11,a12,a12],  
             [a21,a22,a22]] ):
nullspace(A);
rowspace(A);
colspace(A);


# AR6 Section 4.6 Example 4 p196- ; AR7 Section 5.5 Example 7 p267-
with(linalg):
# Example 4:    Find a basis for the vector space spanned by
v1:= vector( [ 1,-2, 0, 0, 3] ):     v2:= vector( [ 2,-5,-3,-2, 6] ):
v3:= vector( [ 0, 5,15,10, 0] ):     v4:= vector( [ 2, 6,18, 8, 6] ):
#
# To do this, find a basis for the row space of the matrix A
A:= matrix( [v1, v2, v3, v4] ):                         rowspace(A);

# Now consider the rowspace of the Echelon form U of A.
U:= gausselim(A);
rowspace(U);
#    which shows that rowSpace(A) = rowSpace(U).


# Find bases for the column spaces colSpace(U) and colSpace(A):
basis([ col(U,1..5) ]); 
# These are just the pivot columns of U (i.e. the columns of U with pivots)
colspace(A);
# Because of the 1's in the last places in some of the last set of vectors,
#                       colSpace(A) <> colSpace(U)
#
# These two bases have the same number of vectors as the number of pivots of U.
# Thus  dim colSpace(A) = dim colSpace(U) = dim rowSpace(U) = dim rowSpace(A).

with(linalg):
A:= matrix( [[a11,a12,a12],  
             [a21,a22,a22]] ):
rank(A);

# Adapted from 
# AR6 Section 5.4 Example 2 p251; AR7 Section 6.5 Example 4 p343
with(linalg):
# Find the coordinates of the vector x  relative to the basis  u1,u2  where
x:= vector([x1,x2]):
u1:=vector([1,1]):        u2:=vector([2,1]):
#
# First check that it is a basis for R^2 using the  rank  function:
B:= augment( u1,u2 ): rank(B);

# Thus  dim( colSpace(B) ) = 2 and so the columns are a basis for R^2.
#
# To find the coordinates of  x  relative to  {u1,u2}, solve  as follows
xB:= linsolve(B,x);

# Find the coordinates of  x  relative to the basis  v1,v2  where
v1:=vector([2,1]):         v2:=vector([1,1]):
C:= augment( v1,v2 ):                                xC:= linsolve(C,x);
# Alternatively:
P:= evalm( inverse(C)&*B );
evalm( P&*xB );

# 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) );


# Here is an example for this theorem:
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) );


# Example 7 - AR6 Section 5.3 p242; AR7 Section 6.3 p320
with(linalg,vector,GramSchmidt,dotprod,norm):
# Use GramSchmidt to find an orthonormal basis for the spaces
# spanned by the list of vectors [v1,v2,v3] where
v1:= vector([1,1,1]):    v2:= vector([0,1,1]):    v3:= vector([0,0,1]):
# This can be done using the function  GramSchmidt  in the  linalg  package:
GramSchmidt( [v1,v2,v3] );
# Each of these should be divided by its norm to give an orthonormal basis
#
# To help understand the process, we also do it as follows:
w1:= v1:                                    u1:= evalm( w1/norm(w1,2) ):
  w2:= evalm( v2 - dotprod(v2,u1)*u1 ):       u2:= evalm( w2/norm(w2,2) ):
  w3:= evalm( v3 - dotprod(v3,u1)*u1 - dotprod(v3,u2)*u2 ):
  u3:= evalm( w3/norm(w3,2) ):
  print (u1,u2,u3);
# Apart from the norms of the vectors, this agrees with the above.

interface(echo=1);

