echo on
%%%%% ch02_exa2.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.xx Inversion versus Solution system *'));
echo on;
%%%% Suppose that we have to compute the product A^-1*B.

n=4;
A=rand(n,n)
B=rand(n,n)
%B=eye(n)	% If we want to compute the inverse of A we take B=I.
S=inv(A)*B	%The direct result of the product A^-1*B

%%% There is a better solution: LU factorization of A, followed by forward- and 
% backward substitutions using columns of B.
S1=zeros(n);
[L,U]=lu(A);
echo off;
for i=1:n
   S1(:,i)=L\B(:,i);
   S1(:,i)=U\S1(:,i)
end
echo on;
S1		% The result of the product A^-1*B solved with the consecutive solutions of
		% the linear system is same as the direct result.
      

echo off
disp(sprintf('**** Example 2.19 Rank-one updating *'));
echo on;
%%%% Solve rank-one modification of the linear system.

b=[2 8 10]'	% Define a system Ax=b.
A=[2 4 -2; 4 9 -3; -2 -3 7]'	% Original matrix for which the system is solved.
Am=[2 4 -2; 4 9 -1; -2 -3 7]'	% The element A(2,3) is changed from -3 to -1.
		% We can select vectors u and v as
u=[0 0 -2]'
v=[0 1 0]'
u*v'		% so that u*v' represents the change in element A(2,3). A-u*v'=b is still valid. 
z=A\u		% This solution has to be done O(n^2) with already computed LU. 
y=A\b		% We had already solved this system. The updated solution is
x=y+((v'*y)/(1-v'*z))*z	% with some additional work of O(n), so the added
			% complexity is O(n^2)+O(n)
 
      
echo off
disp(sprintf('**** Example 2.20 Poor scaling *'));
echo on;
%%%% Sometimes the reason for bad solution is poor scaling.

format short e;
e=eps/10		% eps/10 is a value smaller than machine precision
A=[1 0; 0 e]
b=[1 e]'
x=A\b			% The solution of original system.
cond(A)		% Very ill-conditioned, 
b=b-[0 -e]'	% the perturbation in b for [0 -e]'
x=A\b		% completely changes the solution. 
A=[1 0; 0 1]	% If second row is multiplied by 1/e, in A 
b=[1 1]'		% and in b,
cond(A)		% The system becomes well-conditioned and 
b=b-[0 -e]'	% the same perturbation in b for [0 -e]'
x=A\b		% produces very small change in the solution.
format short;


echo off
disp(sprintf('**** Example 2.21 Cholesky factorization *'));
echo on;
%%%% Cholesky factorization gives us A=LL^T. Try to compute L.

A=[5.0 0 2.5; 0 2.5 0; 2.5 0 2.125]'	% Symmetric positive definite matrix.
% Take only the lover triangular part because A is symmetric. The first 
% column has no prior columns, scale it by square root of the 
% diagonal entry A(1,1).
L=[A(:,1)/sqrt(A(1,1)) [0 2.5 0]' [0 0 2.125]']
% The second column must be updated by subtracting the previous column
% multiplied by its second entry. This entry is 0 so no updating necessary,
% only scaling by the square root of diagonal entry L(2,2).
L=[L(:,1) L(:,2)/sqrt(L(2,2)) [0 0 2.125]']
% The last column has to be updated by subtracting multiples of the 
% previous columns (1.1180 and 0) so only the first column is used.
L=[L(:,1) L(:,2) [0 0 (L(3,3)-L(3,1)*L(3,1))]']
% Finally, the last column is scaled by square root of L(3,3).
L=[L(:,1) L(:,2) [0 0 (L(3,3)/sqrt(L(3,3)))]']