% Matlab code eigscat.m
% For "Applied Numerical Linear Algebra",  Question 4.14
% Written by James Demmel, Oct 25, 1995
%                Modified, Jun  5, 1997
%                Modified, Oct 15, 2004
%
%  Plot eigenvalues of a matrix and eigenvalues of
%  random perturbations of the same matrix. 
%  This illustrates the sensitivity of eigenvalues to perturbations.
%
%  Inputs:
%     a = input matrix
%     err = size of perturbation
%     m = number of perturbed matrices to compute
%
%  Outputs:
%     Three plots where each symbol is the location of an eigenvalue of 
%     a perturbed matrix:
%        'o' marks the location of each unperturbed eigenvalue. 
%        'x' marks the location of each perturbed eigenvalue, where a
%                  real perturbation matrix of norm err is added to a. 
%        '.' marks the location of each perturbed eigenvalue, where a
%                  complex perturbation
%

%  Compute data
%  ea contains eigenvalues of unperturbed matrix a
%  er contains eigenvalues of matrix with real perturbations
%  ei contains eigenvalues of matrix with complex perturbations
function eigscat(a,err,m)

ea=eig(a).'; er=ea; ei=ea;
n=size(a,1);
da = zeros(n*m^2,1);

for i=1:m, 
    r=rand(size(a))-.5; 
    r=(err/norm(r))*r; 
    da(i) = norm(r,2);
    er=[er,eig(a+r).']; 
end
j=1i; 
for i=1:m, 
    r=(rand(size(a))-.5)+j*(rand(size(a))-.5); 
    r=(err/norm(r))*r;
    da(i+m) = norm(r,2);
    ei=[ei,eig(a+r).']; 
end

%
%  Plot data
axis('square'),
figure
hold off,plot(real(er),imag(er),'xr'),hold on,grid,
plot(real(ea),imag(ea),'ok'),title(['Real Perturbations, pertubation: ' num2str(err)]),
% disp('pause (hit return to continue)'),pause
figure
hold off,plot(real(ei),imag(ei),'.b'),hold on,grid,
plot(real(ea),imag(ea),'ok'),title(['Complex Perturbations, pertubation: ' num2str(err)]),
% disp('pause (hit return to continue)'),pause
figure
hold off,plot(real(ei),imag(ei),'.b'),hold on,grid,
plot(real(er),imag(er),'xr'),plot(real(ea),imag(ea),'ok'),
title(['Real and Complex Perturbations, pertubation: ' num2str(err)])
%
%  Compute and print condition numbers
[n,n]=size(a); cnd=1; format short e,
[v,d]=eig(a);
l = inv(v);
for i=1:n, 
    v(:,i)=v(:,i)/norm(v(:,i)); 
    l(:,i)=l(:,i)/norm(l(:,i)); 
end
vi=inv(v); 
for i=1:n, 
    cnd(i)=norm(vi(i,:)); 
    b(i) = (max(da)*n/abs((l(i,:))*v(:,i)));   
end
disp('       Eigenvalues      Condition Numbers'), [diag(d),cnd']
%b
% plot circles
ang=0:0.01:2*pi; 
for i = 1:length(ea)
    if b(i) < 10^3
        xp=min(b(i))*cos(ang);
        yp=min(b(i))*sin(ang);
        plot(real(ea(i))+xp,imag(ea(i))+yp);
hold on;
    end
end
