
% This code is to make a demonstration 
%  of the effect of a small perturbation 
%  on the eigenvalue spectrum of normal and non-normal matrices.

is_matrix_normal='no';        % Is the perturbed matrix normal or not?
                               % Answer 'yes' or 'no'.

                               
% Define the matrix, A, that will be perturbed.
%
M=20;                          % matrix dimension
x=[1:M];                       % auxiliary vector of discretization points
                               
if strcmp(is_matrix_normal,'yes') == 1    % As an example of a normal matrix, take one
                                            %  similar to the symmetric matrix from the 
                                            %  problem about the gold prospector.
    
    h=1/M;                      % Auxiliary parameter that is used to control
                                %  how ill-conditioned A is. It does not have to be
                                %  related to M.
                                %  The smaller h, the more ill-conditioned A is.
    xx=x*h;                     % auxiliary vector used in the next line
	A=1./( (xx'*ones(1,M) - ones(M,1)*xx).^2 + 1 );   
                                % matrix with entries 1/(1+(x_i-x_j)^2);
                                %  this is the matrix to be perturbed
                                
else         % As an example of a non-normal matrix, take the 
             %  "discretized first derivative" matrix.
    
    A=diag(ones(1,M))-diag(ones(1,M-1),1);
    
end


% Define the perturbing matrix, B.
EPSILON=10^(-13);               % the magnitude of the entry(ies) of the perturbing matrix
B=zeros(M);
% B(M,1)=EPSILON;
B=EPSILON*(2*rand(M)-1);


% Plot the eigenvalues sorted in ascending order.
figure(1);
plot(x,sort(eig(A)),x,sort(eig(A+B)),'ro:')   
xlabel('sequential number of eigenvalue','fontsize',12)
ylabel('\lambda  or  Re(\lambda)','fontsize',12)
title('Eigenvalues of the original and perturbed matrices','fontsize',13)
legend('original','perturbed (real part of \lambda only)')

% Plot the difference between the corresponding eigenvalues.
figure(2);
plot(sort(eig(A))-sort(eig(A+B)),'ro:')   
xlabel('Re(\lambda_{orig} - \lambda_{pert})','fontsize',12)
ylabel('Im(\lambda_{orig} - \lambda_{pert})','fontsize',12)
title('\lambda_{original} - \lambda_{perturbed}','fontsize',14)
axis('equal')

    