echo on
%%%%% ch02_exa1.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;
format short;
disp(sprintf('**** Examples from Chapter 02 ****'));
disp(sprintf('**** Example 2.9 Permutations *'));
echo on
%%%% Demonstrate the effects of permutations.

b=[1 2 3]'	% Define a vector
P=[0 0 1; 1 0 0; 0 1 0]	% and a permutation matrix.
P*b		% Left multiplication reorders the rows of a vector or matrix.
inv(P)		% Inverse P is equal to the transpose of P.
P'
A=[1 2 3; 2 5 7; 0 1 3];
x=A\b		% The solution of a system does not change if 
x=(P*A)\(P*b)	% it is multiplied from left (pre-multiply) with a nonsingular 
		% matrix (P is non singular!)

		% Multiplication form right (post-multiply) reorders the columns and 
		% permute also the solution.
x1=(A*P)\b	% x1=(A*P)^-1*b=P^-1*(A^-1*b), and hence  
x=P*x1		% the solution to the original system is x=P*x1


echo off
disp(sprintf('**** Example 2.10 Diagonal scaling *'));
echo on
%%%% Demonstrate the effects of diagonal scaling.

b=[1 2 3]'	% Define a vector
D=diag([3 3 3])	% and the diagonal matrix.
D*b		% Left multiplication scales the rows of a vector or matrix.
inv(D)		% Inverse matrix of D is not equal to the original matrix and to its transpose.
D'
A=[1 2 3; 2 5 7; 0 1 3];
x=A\b		% The solution of a system does not change if 
x=(D*A)\(D*b)	% it is multiplied from left (pre-multiply) with a nonsingular 
		% matrix (D is non singular!)

		% Multiplication form right (post-multiply) scales the columns and does change the solution.
x1=(A*D)\b	% x1=(A*D)^-1*b=D^-1*(A^-1*b), and hence  
x=D*x1		% the solution to the original system is x=D*x1.


echo off
disp(sprintf('**** Example 2.11 Triangular linear system *'));
echo on
%%%% The solution of a triangular system.

b=[2 4 8]'	% Define a right-hand vector
A=[2 4 -2; 0 1 1; 0 0 4]	% and a triangular matrix.
x=A\b		% MatLab solution of the system.	
x3=8/4	% The solution of the system can be generated in a simple way, by components.
x2=(4-(2*1))/1
x1=(2-((4*2)+(-2*2)))/2


echo off
disp(sprintf('**** Example 2.13 Gaussian elimination *'));
echo on;
%%%% Demonstrates how Gaussian elimination works.

b=[2 8 10]'	% Define a system Ax=b.
A=[2 4 -2; 4 9 -3; -2 -3 7]'	% First pivot is 2
M1=[1 -(4/2) -(-2/2);0 1 0;0 0 1]'	%m1_1,j = - a_1,j/a_1,1 for j=1+1,..,3
A1=M1*A		% We should not need to implement the complete matrix-matrix multiplication
		% only the lower submatrix changes (second pivot is 1),
b1=M1*b		% only the lower part of right-hand vector changes.
M2=[1 0 0;0 1 -(1/1);0 0 1]'		%m2_2,j = - a1_2,j/a1_2,2 for j=2+1,..,3
U=M2*A1		% Final result is an upper triangular system
b2=M2*b1	% and modified right-hand vector.
M=M2*M1		% Matrix M is not needed explicitly (2 elementary elimination matrices).
x=U\b2		% The solution x of modified upper triangular system
x=A\b		% is equal to the solution of the original system.
L=inv(M1)*inv(M2)	% Note that M^-1=(M2*M1)^-1 = M1^-1*M2^-1=L1*L2=L
L=inv(M)	% Lower triangular matrix is inverse of the transformation matrix M
A=L*U		% So L*U*x = b, A is factorized on L and U and 
y=L\b		% the solution of lower triangular system (forward substitution) followed by 
x=U\y		% the solution of upper triangular system (backward substitution) again gives x. 


echo off
disp(sprintf('**** Example 2.15 Small pivots *'));
echo on;
%%%% Using finite-precision arithmetic we must avoid small pivots. 

format short e;
e=eps/10		% eps/10 is a value smaller than machine precision
A=[e 1; 1 1]
M=[1 0;-(1/e) 1]
L=[1,0;(1/e),1]
U=M*A		% Upper triangular form.
L*U		% Is not equal to A, so we would have to interchange rows.
%%% If we interchange the rows:
A=[1,1; e, 1]	
M=[1,0;-(e/1),1]
L=[1,0;(e/1),1]
U=M*A		% Upper triangular form.
L*U		% Now L*U is equal to A, what is correct.
format short;


echo off
disp(sprintf('**** Example 2.16 Gaussian elimination with partial pivoting *'));
echo on;
%%%% Partial pivoting improves numerical stability.

b=[2 8 10]'	% Define a system Ax=b.
A=[2 4 -2; 4 9 -3; -2 -3 7]'	% First pivot is originally 2 but 	
P1=[0 1 0; 1 0 0; 0 0 1];		% after interchanging of 1. and 2. rows
A1=P1*A		% the new pivot becomes 4
M1=[1 -(2/4) (2/4);0 1 0;0 0 1]'	%m1_1,j = - a_1,j/a_1,1 for j=1+1,..,3
A1=M1*A1		% We should not need to implement the complete matrix-matrix multiplication
		% only the lower submatrix changes. The second pivot is would be -1/2,
P2=[1 0 0; 0 0 1; 0 1 0];		% but after interchanging of 2. and 3. rows
A2=P2*A1		% the new pivot becomes 3/2      
M2=[1 0 0;0 1 -((-1/2)/(3/2));0 0 1]'		%m2_2,j = - a2_2,j/a2_2,2 for j=2+1,..,3
M=M2*P2*M1*P1		% Matrix M is not needed explicitly 
						% (2 elementary eliminations, 2 permutations).
U=M2*A2			% Final result is an upper triangular system
bm=M*b			% and modified right-hand vector.
L=inv(M)		% L is P1'*L1*p2'*L2, and is not lower triangular.
L*U			% The product L*U is equal to the original matrix A.
P=P2*P1		% Alternatively, we may calculate permutation matrix.
L=P*L			% Is lower triangular, but now
L*U			% L*U is equal to P*A


echo off
disp(sprintf('**** Example 2.17 Small residual *'));
echo on
%%%% Suppose a system Ax=b and 4-digit decimal arithmetic. 
% Small residual does not always signal a good solution.

A=[0.913 0.659; 0.457 0.330]; 
b=[0.254 0.127]';
x=A\b			% The correct MatLab solution      

M1=[1 0; -.457/.913 1]	% Gaussian elimination matrix
U=M1*A			% Triangular system
bm=M1*b			% Modified right-hand side vector
x1=-0.0001/0.0001		% The solution is [x1 x2]'
x2=(0.254 - 0.913*x1)/0.659	% what is wrong.
r=b-A*[x1 x2]'		% Small residual but bad solution.
log10(cond(A))		% log_10(cond(A)) digits were lost in the accuracy of the solution 
			% so the result is inaccurate (ill-conditioned problem).



echo off
disp(sprintf('**** Example 2.18 Gauss-Jordan elimination - no pivoting *'));
echo on;
%%%% Such factorization transforms matrix A in an diagonal form. We should 
% use M*A*x=M*b, and M*A is diagonal. With a diagonal scaling we can get 
% identity matrix, and the solution is equal to the modified right-hand vector
% x=D*M*b.

b=[2 8 10]'	% Define a system Ax=b.
A=[2 4 -2; 4 9 -3; -2 -3 7]'	% First pivot is 2
M1=[1 -(4/2) -(-2/2);0 1 0;0 0 1]'	%m1_1,j = - a_1,j/a_1,1 for j=1,..,3
A1=M1*A		% We should not need to implement the complete matrix-matrix multiplication
			% only the lower submatrix changes (second pivot is 1),
M2=[1 0 0;-(4/1) 1 -(1/1);0 0 1]'		%m2_2,j = - a1_2,j/a1_2,2 for j=1,..,3
A2=M2*A1		% The last pivot is 4
M3=[1 0 0;0 1 0;-(-6/4) -(1/4) 1]'		%m3_3,j = - a2_3,j/a2_3,3 for j=1,..,3
D=M3*A2			% Final result is a diagonal system
M=M3*M2*M1		% Transformation matrix M (3 elementary eliminations).
x=inv(D)*M*b		% The solution x of modified diagonal system
x=A\b			% is equal to the MatLab solution of the original system.