

%Gershgorin Discs Plot
%by Mario Berges

%16 Feb 2007 (Updated 16 Feb 2007)
%This function plots the Gershgorin Discs for the matrix A passed as an argument. It will also plot t

% gerschdisc.m
% 
% This function plots the Gershgorin Discs for the matrix A passed as an argument.
% It will also plot the centers of such discs, and the actual eigenvalues
% of the matrix.

function gershdisc(A)

error(nargchk(nargin,1,1));
if size(A,1) ~= size(A,2)
    error('Matrix should be square');
    return;
end

% For each row (one for each circle), we say:
for i=1:size(A,1)

    % (1 - obtaining center (h,k) of circle)
    % The circle has center in (h,k) where h is the real part of A(i,i) and
    % k is the imaginary part of A(i,i)   :
    h=real(A(i,i)); k=imag(A(i,i)); 

    % (2 - obtaining radius r of circle)
    % Now we try to compute the radius of the circle, which is nothing more
    % than the sum of norm of the elements in the row where i != j
    r=0;
    for j=1:size(A,1)
       if i ~= j 
           r=r+(norm(A(i,j)));
       end    
    end 

    fprintf (' circle %5.0f: center (%5.2e, %5.2e) radius  %5.2e  \n', i, h, k, r );

    % (3 - plot points for the circle using obtained center (h,k) and radius r )    
    % We try to make a vector of points for the circle:
    N=256;
    t=(0:N)*2*pi/N;
    
    % Now we're able to map each of the elements of this vector into a
    % circle:
    hold on;
    plot( r*cos(t)+h, r*sin(t)+k ,'-');

    % We also plot the center of the circle for better undesrtanding:
    hold on;
    % axis square ; hold on % point ploting is faster!
    plot( h, k,'+');

end
%(end for i=1:size(A,1))

% For the circles to be better graphed, we would like to have equal axis:
axis equal;

% Now we plot the actual eigenvalues of the matrix:
ev=eig(A) 
for i=1:size(ev)
    rev=plot(real(ev(i)),imag(ev(i)),'ro') ;
end
hold off;
legend(rev,'Actual Eigenvalues');

end