% Using Moler's examples related to eigenvalues/ vectors, from his
% chapter 10.  Eli Tziperman
close all; clear all


disp('evaluate sensitivity of eigenvalues in two different ways:')
A=gallery(3)
pause
disp('eigenvectors and eigenvalues:')
[X,lambda]=eig(A)
pause
disp('first way, a rough upper bound via cond(X):')
cond_X=condest(X)
pause
disp('second way, a more accurate estimate for each eigenvalues separately:')
disp('via left eigenvectors:')
cond_eig=condeig(A)
pause

disp('perturb A and see the effect on eigenvalues:')
format long
delta=1.e-6
perturbed_eigenvalues=eig(A+delta*randn(3,3))
pause
error_in_eigenvalues=perturbed_eigenvalues-(1:3)'
pause
estimated_error_from_condeig=delta*condeig(A)
format short

pause
disp('sensitivity of a matrix with a multiple degenerate eigenvalue:')
disp('unperturbed matrix:')
A=[2 1 0 0 0; 0 2 1 0 0 ; 0 0 2 1 0 ; 0 0 0 2 1; 0 0 0 0 2]
pause
disp('the conditioning number of A is just perfect:')
conditioning_number_A=cond(A)
pause
disp('but the eigenvalues are still very sensitive:')
eigenvalue_cond_number_A=condeig(A)
pause
disp('eigenvalues and eigenvectors of unperturbed matrix:')
[V,E]=eig(A); V,eigenvalues=diag(E)'
pause
disp('perturbation amplitude and perturbed matrix:')
delta=1.e-8
A=A+[0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; delta 0 0 0 0]
pause
disp('perturbed eigenvalues and eigenvectors:')
[V,E]=eig(A); 
perturbed_eigenvectors=V
pause
perturbed_eigenvalues=diag(E)'
pause

disp('another example of a matrix with sensitive eigenvalues:')
A=gallery(5)
pause
disp('characteristic polynomial is lambda^5=0, but matlab gives:')
eigenvalues=eig(A)
plot(real(eigenvalues),imag(eigenvalues),'r*',0,0,'ko')
axis(0.1*[-1 1 -1 1])
axis square
pause

disp('now repeat with the perturbed matrix, to see that this')
disp('large error is indeed due to roundoff error:')
for k=1:10
  eigenvalues=eig(A+eps*randn(5,5).*A)'; plot(real(eigenvalues),imag(eigenvalues),'r*',0,0,'ko')
  pause
end