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

# AR6 and AR7  Section 1.1 Example 3 p5-
# 
# Consider the linear system of equations
#   x +   y + 2*z = 9 ,
# 2*x + 4*y - 3*z = 1 ,
# 3*x + 6*y - 5*z = 0 .
#
# If the preceding system is written A x = b, the matrix A
# can be defined to Maple
with(linalg,linsolve):
A  :=array([[1,1, 2],
            [2,4,-3],
            [3,6,-5]]);
b  :=array([9,1,0]);
# The solution is obtained
linsolve(A,b);

# AR6 and AR7 Section 1.1 Example 3, leads to the following 
# system in row-echelon form.
with(linalg,matrix,vector,augment,backsub,col):
# The coefficient and augmented matrices 
# of the following row-echelon system
#          x + y +  2  z =     9  ,
#              y - 7/2 z =  -17/2 ,
#                - 1/2 z =   -3/2 ,
#                    5 z = 15
# are given by
Coeff:= matrix([ [1, 1, 2], [0, 1, -7/2], [0, 0, -1/2] ]):

linSystem:= augment( Coeff, vector([9,-17/2,-3/2]) );
# and the answer is given by
backsub(");

# This concerns the "Gaussian elimination" method of solving linear systems. 
with(linalg,matrix,vector,augment,addrow):
# Consider the linear system of equations as above
# This system can be written using an augmented matrix A_b = [A | b] where
A:=matrix([[1,1,2],[2,4,-3],[3,6,-5]]):
b:= vector( [9,1,0] ):                                 A_b:= augment(A,b);
# The example begins by reducing A_b to row-echelon form with row operations.
# It should be emphasised that these operations are likely to be more useful
# in learning some matrix manipulative processes than they are in actually
# solving equations, as the latter is now built into Maple's  linsolve.
# 
# The command to add -A[2,1]/A[1,1] times the first row to the second row is
A1:=addrow(A_b,1,2,-A[2,1]/A[1,1]);
# and similarly to add  -A[3,1]/A[1,1] times the first row to the third row:
A2:=addrow(A1,1,3,-A[3,1]/A[1,1]);
Ech:=addrow(A2,2,3,-A2[3,2]/A2[2,2]);
# If you have had some experience with Maple procedures,
# you could put the above row operations into a procedure.
# (If you are just beginning your Maple you should ignore this next part.)
#
noswaps:=proc(A)
   local A1,A2,Ech;
   A1:=addrow(  A,1,2,-A[2,1]/A[1,1] ):
   A2:=addrow( A1,1,3,-A1[3,1]/A[1,1]):
   Ech:= addrow( A2, 2,3,-A2[3,2]/A2[2,2]); 
  end:
#
with(linalg,equal,backsub,pivot,col);
equal( noswaps(A_b),Ech );
# showing that the result of  noswaps  is  A3.  The answer to the system is
backsub(Ech);

# AR6 and AR7 Section 1.2 Example 3, p13
with(linalg):
# Consider the linear system of equations
#   x1 + 3*x2 - 2*x3        + 2*x5        = 0 ,
# 2*x1 + 6*x2 - 5*x3 - 2*x4 + 4*x5 - 3*x6 = -1 ,
#               5*x3 +10*x4        +15*x6 = 5 ,
# 2*x1 + 6*x2        + 8*x4 + 4*x5 +18*x6 = 6 .
# The associated augmented matrix is
A_b:=matrix([[1,3,-2, 0,2, 0, 0],
               [2,6,-5,-2,4,-3,-1],
               [0,0, 5,10,0,15, 5],
               [2,6, 0, 8,4,18, 6]]):

addrow(A_b,1,2,-A_b[2,1]/A_b[1,1]):
A2:=addrow(",1,4,-A_b[4,1]/A_b[1,1]);
addrow(A2,2,3,-A2[3,3]/A2[2,3]):
A3:=addrow(",2,4,-A2[4,3]/A2[2,3]);
# In this case, the pivot is in the fourth row, and so proceed by
# interchanging the third and fourth rows:
# (All the zero rows of an Echelon form must be at the bottom.)
with(linalg,swaprow,gausselim):
Ech:=swaprow(A3,3,4);

# The  gausselim  function in the linalg package also gives this matrix
equal( gausselim(A_b),Ech );

# There are numerous items of this general kind in AR.
# A treatment of issues involved in floating point approximate
# computations in this is given in AR Chapter 9.


R:=array([[-1,0],
            [ 0,1]]); 
S:=array([[ 0,1],
            [-1,0]]); 


# AR6 Section 1.4 p30; AR7 Section 1.3 p30
with(linalg,matrix,vector,col):
A:= matrix(2,3, (i,j) -> a.i.j):
X:= vector(3, j -> x.j): AX:= evalm( A&*X ):  B:= vector(2, j -> b.j):
print( AX[1] = B[1] );                  print( AX[2] = B[2] );
# Thus the system may be represented by the matrix equation  A*X = B.
# Now make the substitution
#
# The relationship between matrix multiplication and the columns of a matrix
A1:= col(A,1):              A2:= col(A,2):              A3:= col(A,3):
evalm( x1*A1 + x2*A2 + x3*A3 - AX );

# This shows that  AX  is the linear combination  x1*A1 + x2*A2 + x3*A3.


# 
with(linalg):
I3:= diag( 1,1,1):                                   M1:= mulrow(I3,1,c):
P23:= swaprow(I3,2,3):                             E13:= addrow(I3,1,3,m):
print( M1,P23,E13 );


# Example 4  - AR6 Section 1.6, p51; AR7 Section 1.5, p55
with(linalg,matrix,diag,augment,gausselim):
A:= matrix([ [ 1, 2, 3],
               [ 2, 5, 3],
               [ 1, 0, 8] ]):
# Use the function  augment  to construct the matrix  [A : I]
#
Id:= diag( 1,1,1):                       A_I:= augment(A,Id):
#
# Reduce this matrix to Echelon form
AI1:= gausselim(A_I);

# Do a formal back-substitution on this using Maple's row operations 
with(linalg,addrow,mulrow,submatrix,equal):
addrow(AI1,3,2,-3): AI2:=addrow(",3,1,3);
AI3:=addrow(AI2,2,1,-2);
# Now divide each row by the diagonal element.
# The inverse of the matrix A is now the right hand half of this matrix:
AI4:=mulrow(AI3,3,-1):  
  AInverse:=submatrix(AI4,1..3,4..6);
equal( evalm(AInverse&*A), Id );


#\headq{EXERCISE.} At each stage of Example 4 above, let\par
#  \verb$ A.i:= submatrix(AI.i,1..3,1..3);$\par
#  \verb$ E.i:= submatrix(AI.i,1..3,4..6);$
# Verify that\verb$ evalm( E.i &* A ) = A.i$.


A:= array( symmetric,1..3,1..3 ):
with(linalg,transpose,equal):
equal( A,transpose(A) );

interface(echo=1);
