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

echo off
format compact;
format short;
disp(sprintf('**** Additional Examples for Chapter 04 ****'));


echo off
disp(sprintf('**** Example 4.x Eigenvectors of Triangular Matrix *'));
echo on
%%%% demonstrate how eigenvectors can be obtained from a triangular matrix
A=[2 3 4 5; 0 1 3 0; 0 0 7 2; 0 0 0 1]
[E V]=eig(A)	% MatLab eigenvalues and eigenvectors
D=(A-V(3, 3)*eye(4))	% derive eigenvector for 3th eigenvalue
y=D(1:2,1:2)\D(1:2,3)	%upper triangular block solved for upper part of column 3
x=[y',-1,0]'	% compose eigenvector corresponding to 3th eigenvalue
(E(:,3)'*x)/(norm(E(:,3))*norm(x))	% differs only for scalar factor.


echo off
disp(sprintf('**** Example 4.x Deflation *'));
echo on
%%%% if we know an eigenvalue and corresponding eigenvector, we can make the
%%%% problem smaller for 1 dimension by deflation method.
A=[1 0 2; 0 2 1; 2 1 1]	% symmetric matrix
[E V]= eig(A)
v1=V(1,1)	% first eigenvalue and
x1=E(:,1)	% corresponding eigenvector
I=eye(3);	% Identity matrix

%%%% annihilate the subdiagonal elements of the first column
v=x1;	% take the first column
n=norm(v);		% calculate norm and
v(1)= v(1)+ n;	% form vector v
v
H=I-2*v*v'/(v'*v)	% The elementary elimination matrix H1
A1=H*A*inv(H)	% H transforms A in block triangular form with v1=a11
B=A1(2:3,2:3)	% remaining two eigenvalues can be computed 
			% from B with dimension n-1, and eigenvectors
[Y2 V2]=eig(B)
a=(A1(1,2:3)*Y2(:,1))/(V2(1,1)-v1)	% (b'*y2)/v2-v1
x2=inv(H)*[a Y2(:,1)']'	% calculate eigenvector of A x2 from 
		% eigenvector of B y2, and compare with the second column in E


echo off
disp(sprintf('**** Example 4.x Preliminary reduction to Hessenberg Matrix *'));
echo on
%%%% reduces a general matrix to a nearly triangular form (Hessenberg)
%A=[1 0 2; 0 2 1; 2 1 1]	% symmetric matrix
A=[1 2 3 4;3 0 3 1; 7 2 5 3;3 1 1 4]	% general matrix
s=size(A);
s=s(2);		% maximal number of iterations

%%%% annihilate the subdiagonal elements of the first column (starting with 3th row)
v=A(2:s,1);	% take the first column, below diagonal
n=norm(v);		% calculate norm and
v(1)= v(1)+ n;	% form vector v
v
I=eye(s-1);	% identity matrix
H=I-2*v*v'/(v'*v)	% The elementary elimination matrix H
r1=zeros(s,1);	% enlarge H by first row of identity matrix
r1(1)=1;	% construct row of identity matrix
H1=[r1'; 0 H(1,:);0 H(2,:);0 H(3,:);]
A=H1*A*inv(H1)	% Hessenberg form (note that inv(H1)=H1)

%%%% annihilate the subdiagonal elements of the second column (starting with 4th row)
v=A(3:s,2);	% take the second column, below diagonal
n=norm(v);		% calculate norm and
v(1)= v(1)+ n;	% form vector v
v
I=eye(s-2);	% identity matrix
H=I-2*v*v'/(v'*v)	% The elementary elimination matrix H
r1=zeros(s,1);	% enlarge H by first row of identity matrix
r1(1)=1;	% construct first row of identity matrix
r2=zeros(s,1);	% enlarge H by second row of identity matrix
r2(2)=1;	% construct second row of identity matrix
H1=[r1'; r2'; 0 0 H(1,:); 0 0 H(2,:);]
H1*A*inv(H1)	% Hessenberg form (note that inv(H1)=H1)


echo off
disp(sprintf('**** Example 4.x Spectrum slicing by congruent transformation *'));
echo on
%%%% demonstrate how congruent transformation preserves inertia of a symmetric matrix
%A=[1 0 2; 0 2 1; 2 1 1];	% symmetric matrix
A=[2.9766 0.3945 0.4198 1.1159;...
      0.3945 2.7328 -0.3097 0.1129;...
      0.4198 -0.3097 2.5675 0.6079;...
      1.1159 0.1129 0.6079 1.7231;]	% symmetric matrix
s=size(A);
s=s(2);		% matrix size
eig(A)		% MatLab eigenvalues
ro=0.1;
B=A-ro*eye(s)
	% find the factorization LDL' =A sect 2.5.2
[R]=chol(B)
R'*eye(s)*R
      
      
echo off
disp(sprintf('**** Example 4.x Spectrum slicing by Sturm Sequence transformation *'));
echo on
%%%% demonstrate how eigenvalues of a symmetric matrix can be determined
%A=[1 0 2; 0 2 1; 2 1 1];	% symmetric matrix
A=[2.9766 0.3945 0.4198 1.1159;...
      0.3945 2.7328 -0.3097 0.1129;...
      0.4198 -0.3097 2.5675 0.6079;...
      1.1159 0.1129 0.6079 1.7231;]	% symmetric matrix
s=size(A);
s=s(2);		% matrix size
d=zeros(s,1);	% space for minor determinants
eig(A)		% MatLab eigenvalues
ro=2.1;
B=A-ro*eye(s);
echo off	% calculate all determinants of minors
for i=1:s;
   d(i)=det(B(1:i,1:i));	%determinant of minors
end;	
echo on
d		% the number of agreements in successive determinant's signs is equal 
		% to the number of eigenvalues strictly greater than ro.