# file maple1

# Example on pp81-83 of Hill Section 1.7	takes forever !

# Hill Section 1.1 Example 4, p4
# An example of a Markov process, allowing determination of long-term trends
with(linalg):
# Suppose that the population of U.S.A. remains stable at 270 million. Also
# each year     20 per cent of people living in Texas move out of Texas
#               10 per cent of those outside Texas move into Texas.
#
# (a) If y1, y2 are the populations inside and outside Texas respectively
#     at the beginning of a year,  and  u1, u2 are the populations inside and
#     outside Texas respectively at the end of that year, then we can write
#                       u = T y,      where  u = [u1,u2],  y = [y1,y2]  and
T:= matrix([ [80/100, 10/100],
             [20/100, 90/100] ]):
#
# (b) At the beginning of 1991, the populations are given by y where
y:= vector( [ 50000000, 220000000] ):
#
#     At the end of 1991, the populations are u1, u2 where u = [u1,u2] and
u:= evalm( T&*y );
# (c) To find the populations at the beginning of 1990, solve  T x = y:
x:= linsolve( T,y);
# (d) If the populations  z  remain the same for two consecutive years, then
#                        z = T z
#               (Id - T) z = 0,         where
Id:= diag(1,1):
z:= linsolve( Id-T, vector([0,0]) );
#    (In Maple V,  _t[1]  was just  t1 ).
# Since the total population is 270 million  =  _t[1] + 2*_t[1]  =  3*_t[1],
#    the populations are _t[1] = 90 million and 2*_t[1] = 180 million.
#
# (e) To find out what happens over the six year period from 1990 to 1995,
#     repeat (b) six times  or use the following loop:
for i from 1 to 6 do    x:= evalm( T&*x );    od;
# These answers are tending (slowly!) toward the answer  z  in (d).
quit

# Hill Section 1.1 Example 5, p5
with(linalg):
# The coefficient and augmented matrices of the following row-echelon system
#        2 x + 3 y - 2 z = -7,
#              3 y + 2 z =  3,
#                    5 z = 15
# are given by
Coeff:= matrix([ [2, 3,-2],
                 [0, 3, 2],
                 [0, 0, 5] ]):
linsys:= augment( Coeff, vector([-7,3,15]) );
# and the solution is given by
backsub(");
quit

# Hill Section 1.2 Example 4, p14
# This concerns the "Gaussian elimination" method of solving linear systems. 
with(linalg):
# Consider the linear system of equations
#  2*x1 -   x2 -   x3 =  3 ,
# -6*x1 + 6*x2 + 5*x3 = -3 ,
#  4*x1 + 4*x2 + 7*x3 =  3 .
# This system can be written using an augmented matrix A_b = [A | b] where
A:= matrix([ [ 2,-1,-1],
             [-6, 6, 5],
             [ 4, 4, 7] ]):
b:= vector( [3,-3,3] ):                                 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.
# 
pivot:= A[1,1]:          # The pivot is the first non-zero entry in the row
# 
# The command to add  -A[2,1]/pivot  times the first row to the second row is
A1:= addrow( A_b,1,2,-A[2,1]/pivot ):
# and similarly to add  -A[3,1]/pivot  times the first row to the third row:
A2:= addrow( A1, 1,3,-A[3,1]/pivot ); 
pivot:= A2[2,2]:                     Ech:= addrow( A2, 2,3,-A2[3,2]/pivot ); 
# 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 pivot,A1,A2,Ech;
  pivot:= A[1,1]:
  A1:=addrow(  A,1,2,-A[2,1]/pivot ):   A2:=addrow( A1,1,3,-A1[3,1]/pivot ):
  pivot:= A2[2,2]:          Ech:= addrow( A2, 2,3,-A2[3,2]/pivot ); 
end:
#
equal( noswaps(A_b),Ech );
# showing that the result of  noswaps  is  Ech.  The answer to the system is
backsub(Ech);
quit

# Hill Section 1.2 Example 5, p16
with(linalg):
# 
# Consider the linear system of equations
#               6*x3 + 2*x4 - 4*x5 - 8*x6 = 8 ,
#               3*x3 +   x4 - 2*x5 - 4*x6 = 4 ,
# 2*x1 - 3*x2 +   x3 + 4*x4 - 7*x5 +   x6 = 2 ,
# 6*x1 - 9*x2        +11*x4 -19*x5 + 3*x6 = 1 .
# The associated augmented matrix is
A_b:= matrix([ [ 0, 0, 6,  2, -4,-8, 8],
               [ 0, 0, 3,  1, -2,-4, 4],
               [ 2,-3, 1,  4, -7, 1, 2],
               [ 6,-9, 0, 11,-19, 3, 1] ]):
# In this case, the pivot is in the third row, and so we start by
# interchanging the first and third rows:
A1:= swaprow(A_b,1,3);
# and then proceed as before:
A2:= addrow(A1,1,4,-A1[4,1]/A1[1,1]):   A3:= addrow(A2,2,3,-A2[3,3]/A2[2,3]):
A4:= addrow(A3,2,4,-A3[4,3]/A3[2,3]);
# Since all the zero rows of an Echelon form must be at the bottom,
# we finally interchange the last two rows:
Ech:= swaprow(A4,3,4);
# The  gausselim  function in the linalg package also gives this matrix
equal( gausselim(A_b),Ech );
quit

# Hill Section 1.3 Examples 8 - 10 (pp. 29-32).
with(linalg):
# These examples deal with the system of equations
#        2x - 3y + 4z =  5
#        4x      -  z = -2
# whose coefficient matrix is
A:= matrix([ [2,-3, 4],
             [4, 0,-1] ]):
# Let
X:= vector( [x,y,z] ):       AX:= evalm( A&*X ):       B:= vector( [5,-2] ):
evalm(AX) = evalm(B);
# Thus  A*X = B  is the matrix equation of the system of equations.
#
# Now make the substitution
lprint();  subs( x = 3*s-t, y = 2*s+5*t, z = -s+2*t, " );
# The coefficient matrix of this substitution is
C:= matrix([ [ 3,-1],
             [ 2, 5],
             [-1, 2] ]):
# and the coefficient matrix of the new system (in s and t) is
AC:= evalm( A&*C );
#
# The relationship between matrix multiplication and the columns of a matrix
#
# Example 9 (p 31).
A1:= col(A,1):              A2:= col(A,2):              A3:= col(A,3):
evalm( x*A1 + y*A2 + z*A3 - AX );
# This shows that  AX  is the linear combination  x*A1 + y*A2 + z*A3.
#
# Example 10 (p 32).
equal( evalm( augment( A&*col(C,1), A&*col(C,2) ) ), AC );
# This shows that the j'th column of A*C equals A times (j'th column of C)
#
row1:= submatrix(A,1..1,1..3):                 row2:= submatrix(A,2..2,1..3):
equal( evalm( stack( row1&*C, row2&*C) ), AC );
# This shows that the i'th row of A*C equals (i'th row of A) times C
quit

# Hill Section 1.4 Example 5, p41
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 );
quit

# Hill Section 1.4 Example 9, p46/47
with(linalg):
# Use the function  augment  to construct the matrix  [A : I]
A:= matrix([ [ 2, 1, 0],
             [-4,-1,-3],
             [ 3, 1, 2] ]):
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 
AI2:= addrow(AI1,3,2,6):                        AI3:= addrow(AI2,2,1,-1);
# Now divide each row by the diagonal element.
AI4:= mulrow(AI3,1,1/2):                         AI5:= mulrow(AI4,3,2);
# The inverse of the matrix A is now the right hand half of this matrix:
Inverse:= submatrix( AI5,1..3,4..6 ):        equal( evalm(Inverse&*A), Id );

# Hill Section 1.4 Example 10, p48
A_I[3,3]:= 3/2:     # Otherwise, the entries are the same as in Example 9
#
# Reduce the matrix to Echelon form as above
Ech:= gausselim(A_I);
# As the left hand half of the last row of this matrix consists of
#   a row of zeroes, the matrix A has no inverse.
quit

# Hill Section 1.5 Example 3, p59/60
with(linalg):
# Find the LU decomposition of the matrix
A:= matrix([ [ 2,  4, -4],
             [ 1, -4,  3],
             [-6, -9, 10] ]):
# Start with the diagonal entries of  L  equal to 1 and  U  equal to A
L:= diag(1,1,1):                               U:= A:
#
# Now reduce  U  to row-echelon form, storing the INVERSES of the row
# operations used in the matrix  L
L[2,1]:= U[2,1]/U[1,1]:               U:= addrow( U,1,2,-L[2,1]):
L[3,1]:= U[3,1]/U[1,1]:               U:= addrow( U,1,3,-L[3,1]):
L[3,2]:= U[3,2]/U[2,2]:               U:= addrow( U,2,3,-L[3,2]):
print( L,U );
# Verify that    A = L*U
equal( A,evalm(L&*U) );
# Hill Section 1.5 Example 4, p61
# Solve the equation  A*X = B  (i.e. L*U*X = B)  where
B:= vector( [12,-21,-24] ):
#
# Let  U*X = Y  (and so  L*Y = B)  where Y is found by forward substitution:
y1:= B[1]:                                 y2:= B[2] - L[2,1]*y1:
y3:= B[3] - L[3,1]*y1 - L[3,2]*y2:         Y:= vector( [y1,y2,y3] );
# The solution is now given by
X:= backsub( augment(U,Y) );
quit

# Hill Section 1.5 Example 6, pp 63-65
# An example of the  PA = LU  decomposition of a matrix
with(linalg):
A:= matrix([ [0,0, 3],
             [2,4,-1],
             [6,5, 5] ]):
# The following row operations will reduce  A  to Echelon form:
U:= swaprow(A,1,2):
U:= addrow(U,1,3,-U[3,1]/U[1,1]):
U:= swaprow(U,2,3):         # (The value of this matrix U is given later)
#
# Alternatively, perform the above row interchanges on the identity matrix
P:= diag(1,1,1):           P:= swaprow(P,1,2):        P:= swaprow(P,2,3);
PA:= evalm( P&*A );
# Now perform an LU decomposition on PA, as in Example 3 above
L:= diag(1,1,1):           L[2,1]:= PA[2,1]/PA[1,1]:
print( L,U );
equal( PA,evalm(L&*U) );

# Hill Section 1.5 Example 7, p 65
# The above can be used (as in example 4) to solve the equation  A*X = B where
B:= vector( [6,-10,-7] ):
# by solving  P*A*X = P*B  where
PB:= evalm( P&*B );
quit

# Hill Section 1.6 Theorem 1.69, p70.
A:= array( symmetric,1..3,1..3 ):
with(linalg):
equal( A,transpose(A) );
quit

# Hill Section 1.6 Example 5, p73.
with(linalg):
# The LU decomposition of  the matrix
A:= matrix([ [ 1, 1,  0, 0],
             [ 2, 4,  1, 0],
             [ 0, 6,  7, 1],
             [ 0, 0, 12, 2] ]):
# is given by
L:= diag(1,1,1,1):                               U:= A:
L[2,1]:= U[2,1]/U[1,1]:               U:= addrow( U,1,2,-L[2,1]):
L[3,2]:= U[3,2]/U[2,2]:               U:= addrow( U,2,3,-L[3,2]):
L[4,3]:= U[4,3]/U[3,3]:               U:= addrow( U,3,4,-L[4,3]):
#
print( L,U );          # these are both banded matrices.
inverse(A);             # this is NOT a banded matrix.
quit

# Hill Section 1.7 Example 1, p75
# We set the value of Digits to 5, and try to solve the equations
eqn1:= 10^(-8)*x+y=1;                                   eqn2:= x+y=2;
vars:= {x,y}:            Digits:= 5:            fsolve( {eqn1,eqn2},vars );
# If we pivot on the first coefficient of  eqn2,  we get
eqn:= evalf(eqn1-eqn2/10^8);                      fsolve( {eqn,eqn2},vars );
# which is the correct answer to five significant figures.
#
# However, if we try to pivot on the first coefficient of  eqn1,  we get
eqn:= evalf(eqn2-10^8*eqn1);                     fsolve( {eqn1,eqn},vars );
# which is not the correct answer.
# Note that Maple gives the correct answer (in both cases) with the value
Digits:= 10:
quit

# Example 6 of Hill Section 1.7, p79:      Solve
eqn1:= x + 2*y = 3;                             eqn2:= 1.001*x + 2*y = 3+f;
vars:= {x,y}:                  Digits:= 5:
for f from 1/1000 by 1/1000 to 3/1000 do
    print( evalf(eqn2),solve( {eqn1,eqn2},vars) );
od:
quit

with(linalg):
A:= matrix([ [ 1,2], [ 1.001,2] ]):
cond(A);
# Note that Maple's function  cond  may have a second parameter
# which yields different answers.
quit

