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

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


echo off
disp(sprintf('**** Example 4.3 Characteristic Polynomial *'));
echo on
%%%% calculate det(A-rI)=0
%A=[1 0 2; 0 2 1; 2 1 1] % symmetric matrix 
A=[3 1; 1 3]	% symmetric matrix
p=poly(A)	% coefficients of characteristic polynomial
r = roots(p)	% eigenvalues 


echo off
disp(sprintf('**** Example 4.7 Eigenvalue Sensitivity *'));
echo on
%%%% demonstrate the eigenvalue sensitivity
A=[-149 -50 -154; 537 180 546; -27 -9 -25]
p=poly(A)	% coefficients of characteristic polynomial
r = roots(p)	% distinct eigenvalues -> A id nondefective and diagonalizable 
A'*A
A*A'	% but A is not normal -> left and right eigenvectors are different
[X,RR]=eig(A)	% MatLab right eigenvectors and eigenvalues
[Y,LR]=eig(A')	% MatLab left (right of A') eigenvectors and eigenvalues
cond(X)	% overall condition number of the eigenvalues
Y(:,1)'*X(:,1)	% and therefore left and right eigenvectors corresponding 
Y(:,2)'*X(:,2)	% the same eigenvalue are 
Y(:,3)'*X(:,3)	
A(2,2)=179.99	% changing slightly an element of A 
p=poly(A);
r = roots(p)	% gives totally different (even complex) eigenvalues


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.11 Power method with normalization*'));
echo on
%%%% derives approximate normalized eigenvector with maximal modulus and eigenvalue
%A=[1 0 2; 0 2 1; 2 1 1];	% symmetric matrix
A=[3 1; 1 3] 	% symmetric matrix
s=size(A);
s=s(2);
%xi=rand(s,1)	%random initial approximation of eigenvector
xi=[0 1]'
xi=A*xi	% generate next vector
v=norm(xi,inf)	% approximate maximal eigenvalue
xi=xi/v		% approximate normalized eigenvector
echo off;	% after 8 additional iterations
for i=1:8
  xi=A*xi;	% generate next vector
  v=norm(xi,inf);	% approximate maximal eigenvalue
  xi=xi/v;		% approximate normalized eigenvector
end
echo on;
v		% approximate maximal eigenvalue after 9 iterations
xi		% approximate normalized eigenvector


echo off
disp(sprintf('**** Example 4.13 Power method-inverse Iteration, Reyleigh Quotient*'));
echo on
%%%% derives approximate normalized eigenvector with minimal modulus and eigenvalue
%A=[1 0 2; 0 2 1; 2 1 1]	% symmetric matrix
A=[3 1; 1 3] 	% symmetric matrix
s=size(A);
s=s(2);
%xi=rand(s,1)	%random initial approximation of eigenvector
xi=[0 1]';
xi=A\xi	% solve the system at each iteration(LU factorization only once)
rv=(xi'*A*xi)/(xi'*xi)	% Reyleigh Quotient for approximation of eigenvalue
v=norm(xi,inf)	% approximate minimal eigenvalue
xi=xi/v		% normalized approximate eigenvector
echo off;	% after 5 additional iterations
for i=1:5
 xi=A\xi;	% solve the system at each iteration(LU factorization only once)
 rv=(xi'*A*xi)/(xi'*xi);	% Reyleigh Quotient for approximation of eigenvalue
 v=norm(xi,inf);	% approximate minimal eigenvalue
 xi=xi/v;		% approximate eigenvector
end
echo on;
v		% approximate minimal eigenvalue 6 iterations
rv		% approximate minimal eigenvalue with Reyleigh Quotient
xi		% approximate normalized eigenvector with minimal modulus


echo off
disp(sprintf('**** Example 4.14 Reyleigh Quotient iteration*'));
echo on
%%%% combines all mentioned futures to derive eigenvalue and corresponding 
%%%% eigenvector that is closer to the initial approximate vector in directions. 
A=[1 0 2; 0 2 1; 2 1 1]	% symmetric matrix
xi=[.5 -1 .2]';	% this approximation will yield to the 1. eigenvalue
%A=[1.5 -0.5; -0.5 1.5] 	% symmetric matrix
%A=[3 1; 1 3] 	% symmetric matrix
%xi=[.807 .397]';	% this approximation will yield to the 2. eigenvalue
s=size(A);
s=s(2);
%xi=rand(s,1)	%random initial approximation of eigenvector

rv=(xi'*A*xi)/(xi'*xi)	% Reyleigh Quotient for shift
xi=(A-(rv*eye(s)))\xi	% solve the system at each iteration(LU factorization only once)
v=norm(xi,inf)	% approximate eigenvalue
xi=xi/v		% and corresponding eigenvector
echo off;	% after 2 additional iterations
for i=1:4
 rv=(xi'*A*xi)/(xi'*xi);	% Reyleigh Quotient for approximation of eigenvalue
 xi=(A-(rv*eye(s)))\xi;	% solve the system at each iteration(LU factorization only once)
 v=norm(xi,inf);	% approximate eigenvalue
 xi=xi/v;		% and corresponding eigenvector
end
echo on;
v		% approximate eigenvalue after 3 iterations (large->the exact solution)
rv		% approximate eigenvalue with Reyleigh Quotient
xi		% corresponding normalized eigenvector

echo off
disp(sprintf('**** Example 4.?? Orthogonal Iteration *'));
echo on
%%%% reduces the matrix A into triangular or diagonal form
%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);

X=eye(s);	% Place for eigenvectors
[Qi R]=qr(X)		% orthogonal-triangular decomposition of A=Q*R
X=A*Qi		% form next iteration
echo off;	% after 5 additional iterations
for i=1:5	
 [Qi R]=qr(X);		% orthogonal-triangular decomposition of A=Q*R
 X=A*Qi;		% form next iteration
end
echo on;	
		% the resulting Q after 6 iterations 
Qi		% App. eigenvectors
R       % Triangular matrix with app. eigenvalues on the diagonal
[QM lM]=eig(A)  % MatLab eigenvetors and eigenvalues


echo off
disp(sprintf('**** Example 4.15 QR Iteration *'));
echo on
%%%% reduces the matrix A into triangular or diagonal form
%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);

Q=eye(s);	% Place for eigenvectors
[Qi R]=qr(A)		% orthogonal-triangular decomposition of A=Q*R
A=R*Qi		% form reverse product for QR iteration
Q=Q*Qi;
echo off;	% after 2 additional iterations
for i=1:2	
 [Qi R]=qr(A);		% orthogonal-triangular decomposition of A=Q*R
 A=R*Qi;		% form reverse product for QR iteration 
 Q=Q*Qi;		% accumulate approximate eigenvectors
end
echo on;	
		% the resulting A and Q after 3 iterations 
A		% diagonal entries are equal to the approximate eigenvalues
Q		% columns are equal to the approximate eigenvectors


echo off
disp(sprintf('**** Example 4.16 QR Iteration with Shifts *'));
echo on
%%%% reduces matrix A into triangular or diagonal form
% for shift take 1 that is equal to the lower right diagonal element of A
%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);

Q=eye(s);	% Place for eigenvectors
As=A-A(s,s)*eye(s)	% symmetric shifted matrix
[Qi R]=qr(As)	% orthogonal-triangular decomposition of A=Q*R
As=R*Qi+A(s,s)*eye(s)	% form reverse product for QR iteration, add shift
A=As;
Q=Q*Qi;
echo off;	% after 2 additional iterations
for i=1:2
 As=A-A(s,s)*eye(s);	% symmetric shifted matrix
 [Qi R]=qr(As);	% orthogonal-triangular decomposition of A=Q*R
 As=R*Qi+A(s,s)*eye(s);	% form reverse product for QR iteration
 A=As;
Q=Q*Qi;		% accumulate approximate eigenvectors
end
echo on;	
		% the resulting A and Q are closer to the exact values after 3 iterations 
A		% diagonal entries are equal to the approximate eigenvalues
Q		% columns are equal to the approximate eigenvectors



echo off
disp(sprintf('**** Example 4.19 Jacobi Method *'));
echo on
%%%% reduce symmetric matrix to a near diagonal matrix
A=[1 0 2; 0 2 1; 2 1 1];	% symmetric matrix
[E,R]=eig(A)	% MatLab eigenvectors and eigenvalues
J=eye(3);	% Place for complete Jacobi rotation

%% anihilate simetricaly placed entries (1,3) and (3,1) using plane rotation
Ji=eye(3);	% Place for step Jacobi rotation
t=sym('t');
k=(A(1,1)-A(3,3))/A(3,1);
p=1+t*k - t^2
q=solve(p);
a=q(1); b=q(2);
if (abs('a') <= abs('b')) 
   a=b;		% take smaller absolute value of roots
end;
c=1/sqrt(1+a^2);
s=c*a;
Ji(1,1)=c; Ji(3,3)=c;
Ji(3,1)=-s; Ji(1,3)=s;
A=Ji'*A*Ji
J=J*Ji
%% Annihilate symmetrically placed entries (1,2) and (2,1) using plane rotation
echo off;
Ji=eye(3);	% Place for step Jacobi rotation
k=(A(1,1)-A(2,2))/A(2,1);
p=1+t*k - t^2;
q=solve(p);
a=q(1); b=q(2);
if (abs('a') <= abs('b')) 
   a=b;		% take smaller absolute value of roots
end;
c=1/sqrt(1+a^2);
s=c*a;
Ji(1,1)=c; Ji(2,2)=c;
Ji(2,1)=-s; Ji(1,2)=s;
echo on;
A=Ji'*A*Ji
J=J*Ji
%% Annihilate symmetrically placed entries (2,3) and (3,2) using plane rotation
echo off;
Ji=eye(3);	% Place for step Jacobi rotation
p=1+t*(A(2,2)-A(3,3))/A(3,2) - t^2;
q=solve(p);
a=q(1); b=q(2);
if (abs('a') <= abs('b')) 
   a=b;		% take smaller absolute value of roots
end;
c=1/sqrt(1+a^2);
s=c*a;
Ji(2,2)=c; Ji(3,3)=c;
Ji(3,2)=-s; Ji(2,3)=s;
echo on;
A=Ji'*A*Ji
J=J*Ji
%% Annihilate symmetrically placed entries (1,3) and (3,1) using plane rotation
echo off;
Ji=eye(3);	% Place for step Jacobi rotation
p=1+t*(A(1,1)-A(3,3))/A(3,1) - t^2;
q=solve(p);
a=q(1); b=q(2);
if (abs('a') <= abs('b')) 
   a=b;		% take smaller absolute value of roots
end;
c=1/sqrt(1+a^2);
s=c*a;
Ji(1,1)=c; Ji(3,3)=c;
Ji(3,1)=-s; Ji(1,3)=s;
echo on;
A=Ji'*A*Ji
J=J*Ji
%%% The annihilation can proceed resulting in a diagonal matrix