%%%%%%    MATH 4242 -- LAB 1
%  The goal of this lab is to become familiar with the use of Matlab.  
%  Lines like this one which begin with a % sign are explanations.
%  Other lines show commands that you should type into Matlab

%  BASIC MATRIX OPERATIONS

% To enter a matrix you can use commas or
% spaces to separate entries in the rows
% and semicolons to go the next row. Enter the first
% line below and then hit the Enter or Return key to
% execute the command.  Then enter the second line and do the same.

a = [1,2,3;4,5,6;7,8,9]
a = [1 2 3;4 5 6;7 8 9]

% This enters the matrix and give it a name: a
% Now you try to enter the following matrices
%b =
%   -1
%     0
%     2

%c =
%     1     2
%     3     4
%    5     6

%d =
%   -1     0    -2
%   -3     0    -5
%   -6     0    -7

% You can easily add and multiply matrices
a + d
a*d
a+b
a*b
b*a

% You can read off or modify individual matrix entries as follows:
a(1,1)
a(1,2)
a(2,3)
a(3,3)=10

% You can specify a range of indices for rows or columns using a colon.
% For example to get the first two entries of the third row, type:
a(3,1:2)

% The entire third row is given by
a(3,:)


% You can reset the third row to some other vector as follows:
a(3,:) = [-7 -8 -9]

% See if you can figure out how to replace the second column by
%  [1;0;0]

% There are built in commands to produce an identity matrix, 
% zero matrix or random matrix
eye(5)
zeros(2,3)
rand(5,7)

% The random matrix has entries uniformly distributed between 0 and 1. This
% can be changed, for example to 0 and 10, by multiplication.  If you want a matrix 
% of integers you can apply the floor() command to throw away fractional parts.
r = 10*rand(5,5)
floor(r)

% Once you have entered a matrix you can find its reduced row ecehelon form
a=[1 2 3;4 5 6;7 8 9]
rref(a)

% For square matrices you can compute determinant and inverse as follows
det(a)
inv(a)

% The results are sometimes untrustworthy due to numerical imprecision.
% A singular matrix may become nonsingular due to roundoff errors of 
% floating point numbers.

% You can carry out your own computation of the inverse using the rref commmand.
% First make a big matrix [A I] then row reduce to get [I inv(A)]

big = [a eye(3)]
rref(big)

% From this you can see that the matrix is not invertible.
% On the other hand, if the matrix is invertible, the inv command
% works very quickly and efficiently.
a(3,3)=10
ai = inv(a)
% You can check the answer by matrix multiplication
a*ai

% Here is a larger example
r= rand(25)
ri =  inv(i)
r*ri



%%%%     SOLVING LINEAR EQUATIONS
% There are several way to get Matlab to solve AX=B
% First you can construct the augmented matrix [A B] and use rref
aug = [a b]
rref(aug)

% This shows that there is no solution (this b is not in the column space)
% Let's change the matrix so that it's invertible and try again
a(3,3)=10
aug = [a b]
rref(aug)

% the solution is in the last column
% In this invertible case you could also use the formula X=inv(A)*B
inv(a)*b

% Finally, Matlab has a "matrix division" command which solves 
% linear equations (sometimes with unexpecteed results)
% To solve AX=B you would use a backslash \ as follows
a\b

% Intuitively, you are dividing B by A, hence the backslash \ instead of the division symbol /
% This works fine in the invertible n x n case.  Let's try a few other cases.
% First we change A back to a singular matrix

a(3,3)=9
a\b

% As you can see, Matlab forges ahead and gives you an answer even though
% no solution exists.  The answer is the "best approximation" or "least squares"
% solution to the problem, i.e., a vector which comes as close as possible to 
% being a solution.  This is either convenient or misleading depending on 
% what you are trying to accomplish. But at least you get warned.

% Next we try a problem with infinitely many solutions.

p = [1 2 3]
q = [6]

% the solution set of pX=q is a plane in three-dimensional space

p\q

% Matrix division gives you a single solution.
% This is only one particular solution to the equation.  To find the general
% solution you need to add on a general vector in the null space.  There
% is a convenient built-in command to find a basis for the nullspace
basis = null(p)

% The two columns in this matrix are the basis.
% We can add any linear combination of the column to the previous solution
% to get another solution.  For example
xp = p\q
xn = 5*basis(:,1) + 7*basis(:,2)
x = xp+xn
p*x

% To demonstrate the impressive speed of the numerical equation solving
% software, try creating a random 100x100 matrix, a, and a random 100x1
% column vector,b, and then solving ax=b.
a = rand(100,100)
b = rand(100,1)
a\b

%%%%     MATRIX FACTORIZATION

% The PA = LU factorization is built-in to Mattlab.

a=[1,2,3;4,5,6;7,8,9]
[l,u,p]= lu(a)

% You can check the result by matrix multiplication

p*a
l*u

% You can convert LU into LDU by extracting the pivots from U
% The command diag applied to a matrix removes the diagonal and stores it
% as a vector.  The same command applied to a vector constructs a diagtonal
% matrix.  So if you use the diag command twice you extract the diagonal
% matrix part of the matrix you started with.

diag(u)
d = diag(diag(u))
 
% now if we replace u by inv(d)*u we get the ldu factorization

u1 = inv(d)*u

% now d*u1 is the original u and l*d*u1 is pa

d*u1
l*d*u1

%%%    COLUMN AND ROW SPACES
% The reduced row echelon form was our main tool for finding
% bases for the row and column spaces.  The nonzero rows of
% rref(a) give the "reduced" basis for the row space while
% the pivot columns of rref(a) show which of the columns of
% the original matrix can be used as a basis for C(A)

a = [1 0 1 0 1;0 1 0 1 0;1 -1 0 1 -1;2 0 1 2 0]
r = rref(a)
redrowbasis = r(1:3,:)
colbasis = a(:,1:3)

% On the other hand we can find a reduced basis for C(A)
% and an unreduced one for R(A) by using the transposed matrix.
% The symbol for a transpose of a in Matlab is a' 
a'
rtrans = rref(a')

% The first three rows are a basis for the row space of a' so
% if we transpose them back we have a reduced basis for C(a)
redcolbasis = rtrans(1:3,:)'

% Since the first three columns of a' are the pivot columns,
% the first three rows of a are a basis for R(a)
rowbasis = a(1:3,:)

% Thus we can get two different basis for each space.
% Of course one of the big theorems was that the dimensions of both
% spaces are equal -- the dimension is just the rank of a.
% Matlab has a command to compute the rank
rank(a)



%   SOME PROBLEMS TO SOLVE
% See if you can figure out how to make Matlab solve these problems. 
% Or make up some for yourself.  Have fun.

% 1. Recall that you can do Gauss elimination by multiplying your
%    matrix by so-called elementary matrices.  Enter the 4x4 matrix a below
%    and then find a 4x4 elementary matrix E21 such that E21*a will have a 
%	 zero in position (1,2).  Can you figure out how to enter the matrix
%	 e21 into Matlab without laboriously typing the whole thing ?  Once you have
%    e21 multiply e21*a and see if it works.  Continue in this way to
%    construct matrices e32 and e43 so that e43*e32*e21*a is upper triangular
%    Finally, by using the inverses of the elementary matrices, find the
%    matrix L of the LU decomposition.  
     2    -1     0     0
    -1     2    -1     0
     0    -1     2    -1
     0     0    -1     2

% 2. Show that the vectors below do not form a basis
%    for R4.  Find a basis for the subspace of
%    R4 which they span.  What is the dimension of this subspace ?
%    Since they don't form a basis they must be linearly dependent
%    Find some numbers c1,c2,c3,c4 such that
%					       c1 X1 +c2 X2 + c3 X3 + c4 X4 = 0
%
%         (1,2,3,-1)  (0,1,4,5)   (1,0,1,-1)   (3,-1,-7,-18)

% 3.  Using random square matrices, A and B, of some size, 
%    test the formulas below:
%    inv(AB) = inv(B)*inv(A)    
%    (AB)' = B'A'
%    det(AB)= det(A)det(B)
%    AB not equal BA (usually)

% 4.  Find a cubic polynomial f(x) = a + b x + c x^2 + d x^3
%     whose graph passes through the four given points below:
%     (-1,5)  (1,3)  (2,0)  (3,0)
%     In other words, you want to satisfy f(-1)=5, f(1)=3, etc.
%     Find out what equations (a,b,c,d) must satisfy and use Matlab
%     to solve them.