%%  An Eigenvalue Sensitivity Example
%
% In this example, we construct a matrix whose eigenvalues are
% moderately  sensitive  to  perturbations  and  then  analyze that
% sensitivity. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B = [3 0 7 ; 
     0 2 0 ; 
     0 0 1 ]; 

[V,D]=eig(B);                   % eigenvalues are 1, 2 and 3 
cond(V)                         % well conditioned eigenvectors   
pause

% generate a similarity transformation to make eigenvalues more sensitive
 
L = [1 0 0 ; 
     2 1 0 ; 
    -3 4 1 ];
M = L\L';

det(M)                            % det(M) = 1 
cond(M)                           % moderately ill-conditioned
pause
   
% transform the matrix A - the eigenvalues remain the same

A = M*B/M;

[W,D]=eig(A);
diag(D)                           % computed in another order by EISPACK
cond(W)                           % moderately ill-conditioned eigenvectors    
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Because  det(M) = 1, the elements of  A  would be exact integers
% if there were no roundoff.  So, round the matrix.

A1 = round(A);                    % A1 = A + E, E - from rounding
norm(A-A1)                        % norm(E) = 1e-12
 
[X,D1] = eig(A1);
cond(X)                           % the conditioning is the same as for W
pause

% since eigenvectors were ill-conditioned, eigenvalues can be sensitive 

format long
diag(D1)                          % the last four significant digits are 
                                  % contaminated  by  roundoff  errors
format short
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analyze sensitivity of individual eigenvalues using left eigenvectors.

Y = inv(X');                      % left eigenvectors
 
% compute the condition number of individual eigenvalues

for j = 1:3, 
    c(j) = norm(Y(:,j))*norm(X(:,j)); 
end
   
c                                           % lambda1 is the most sensitive
pause


% consider a perturbation given by the directions x1 and y1
 
E = -1.e-6*Y(:,1)*X(:,1)'

eig(A + .4*E)                               % large change of eigenvalues
eig(A + .5*E)                               % complex eigenvalues!