% PS1_2.m is the m-file answer for Problem set #1, question 2.
% Written by T. Kenna 9809041213.

% Note: See PS1_1.m for explanation of answer format

clc
ans1=['Start by loading the data and obtaining the Determinant and Rank by'];
ans2=['using the det(A1) and rank(A1) commands. <return>'];
disp(ans1),disp(ans2)			% display the answer strings
pause					% pause so I can read it

clc
load A1.dat       			% load the data
load b1.dat       			% load the data
Determinant=det(A1)			% find the determinant of A1
Rank=rank(A1)				% find the rank of A1
ans1=['A determinant of zero is bad news, it means that the matrix is'];
ans2=['singular. The rank of 6 is less than the full dimensionality (7) of'];
ans3=['the matrix. What happens when you try and solve A1*x=b1 directly'];
ans4=['using the command A1\b1? <return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4)
pause

clc
A1\b1
ans1=['Ugh! Just as the determinant and rank suggested, it can not be solved'];
ans2=['uniquely because we have 6 equations and 7 unknowns. As a final nail'];
ans3=['in the coffin, we can check the condition of A1 by using cond(A1).'];
ans4=['<return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4)
pause

clc
cond(A1)
ans1=['Talk about ill-conditioned!  While we can not come up with a unique'];
ans2=['solution, by doing a singular value decomposition using the'];
ans3=['[U S V]=svd(A1,0) command, we are on our way to obtaining  one of'];
ans4=['many solutions that will work. Then, by using the command'];
ans5=['diagonal_of_matrix_S=diag(S), we can obtain the diagonal of the S'];
ans6=['matrix. <return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4),disp(ans5),disp(ans6)
pause

clc
[U S V]=svd(A1,0);
diagonal_of_matrix_S=diag(S)
ans1=['An examination of the diagonal of S reveals that the 7th element is'];
ans2=['zero.  Since the inverse of 0 is infinity, we must take a few extra'];
ans3=['steps in order to solve for x. <return>'];
disp(ans1),disp(ans2),disp(ans3)
pause

clc
diagonal_of_matrix_S(7,1)=1
ans1=['First, set the value of the last element in the diagonal of matrix S'];
ans2=['to 1 by using the command diagonal_of_matrix_S(7,1)=1, notice the 1'];
ans3=['in place of the zero for the last element.  <return>'];
disp(ans1),disp(ans2),disp(ans3)
pause

clc
inverse_of_diagonal_of_matrix_S=1./diagonal_of_matrix_S
ans1=['Next, take the inverse of the new and improved diagonal of matrix S'];
ans2=['by using the command inverse_of_diagonal_of_matrix_S=1./diagonal_of_'];
ans3=['matrix_S. <return>'];
disp(ans1),disp(ans2),disp(ans3)
pause

clc
inverse_of_diagonal_of_matrix_S(7)=0
ans1=['Then make the 7th element equal zero by using the command'];
ans2=['inverse_of_diagonal_of_matrix_S(7)=0.  <return>'];
disp(ans1),disp(ans2)
pause

clc
W=diag(inverse_of_diagonal_of_matrix_S)
ans1=['Now, convert the inverse matrix into a diagonal matrix by using the'];
ans2=['command W=diag(inverse_of_diagonal_of_matrix_S). <return>'];
disp(ans1),disp(ans2)
pause

clc
Best_guess_inverse_A1=V*W*U'
ans1=['Finally, with the offending element in matrix S taken care of, we can'];
ans2=['now calculate a best guess inverse of A1 by solving for the triple'];
ans3=['product using the command Best_guess_inverse_A1=V*W*U''. Note we are'];
ans4=['using W instead of S. <return>'];
disp(ans1),disp(ans2),disp(ans3),disp(ans4)
pause

clc
ans1=['Remembering that any matrix times its inverse is the identity matrix,'];
ans2=['what do you get when you multiply Best_guess_inverse_A1*A1? <return>'];
disp(ans1),disp(ans2)
pause

clc
Best_guess_inverse_A1*A1
ans1=['Use the command round(Best_guess_inverse_A1*A1) for a simpler view.'];
ans2=['<return>'];
disp(ans1),disp(ans2)
pause

clc
round(Best_guess_inverse_A1*A1)
ans1=['Get the idea? Now use the calculated inverse of A1 to solve for x.'];
ans2=['by using the command x=Best_guess_inverse_A1*b1. <return>'];
disp(ans1),disp(ans2)
pause

clc
x=Best_guess_inverse_A1*b1
ans1=['As a check, you should now solve for b1 using the command A1*x.'];
ans2=['Compare this to the original b1 matrix. <return>'];
disp(ans1),disp(ans2)
pause

clc
calculated_b1_vs_the_original_b1=[A1*x,b1]
ans1=['A side by side comparison reveals that they are the same. <return>'];
disp(ans1)
pause

clc
ans1=['End of Problem 2'];
disp(ans1)