% Matlab code polyplot.m 
%   function polyplot(r,e,m)
% For "Applied Numerical Linear Algebra", Question 1.20
% James Demmel, Aug 1995, modified May 1997, modified Sep 2004
%
% Form the coefficients of a polynomial specified by its roots,
%   and repeatedly add small perturbations to the coefficients,
%   plotting the resulting perturbed roots
%
% Inputs:
%   r = Column vector of polynomial roots
%   e = Scalar maximum relative perturbation to make to each coefficient
%      note number of coeffs is 1 + number of roots
%   m = Scalar number of random polynomials to generate
%   figno=figure handle number
%   sr=subplot row index
%   sc=subplot col index
%   sp=subplot index   
%
% Output:
%   Plot of perturbed roots
function polyplot(r,e,m,figno,sr,sc,sp)
%
% Generate polynomial coefficients
p=poly(r);

% Prep for plotting circles
param=0:.01:2*pi;
cirx=cos(param);
ciry=sin(param);
cent=ones(size(param));

% Compute condition numbers....
ap=abs(p);
for index=1:(length(p)-1),
    dp(index)=index*p(index+1);
end
p
dp
c=polyval(ap,r)./abs(polyval(dp,r));
c=[10000,500,30,16.5,5.7,2.8,1.5,1.1,1.3,2.7,8.3].*c;

% Compute min and max roots for plotting
minx = min(real(r)); maxx = max(real(r));
miny = min(imag(r)); maxy = max(imag(r));
% Set up plot
figure(figno)
axis equal
if (sp==1),
    hold off, clf, 
end
subplot(sr,sc,sp)
% Generate m random polynomials
r1save=[];
for i=1:m,
%  Add random relative perturbation of at most e to each coefficient
   p1=p.*(ones(size(p))+e*2*(rand(size(p))-.5));
%  Compute and save perturbed roots
   r1=roots(p1);
   r1save=[r1save;r1];
%  Update min and max for plotting
   
   %disp(p);
   %disp(p1);
   %for index=1:size(r1,1),
   %    
   %    fprintf('%f + %fi  %f   ||   ', real(r1(index)),imag(r1(index)), abs(r1(index)))
   %end
   %printf('\n')
   
   minx = min([real(r1);minx]); maxx = max([real(r1);maxx]);
   miny = min([imag(r1);miny]); maxy = max([imag(r1);maxy]);
end
%  If roots all lie in right halfplane, and
%  vary greatly in magnitude, use a semilog plot
if (minx>0 & minx <= .002*maxx),
%  Plot unperturbed and perturbed roots
   semilogx(real(r),imag(r),'kx'), hold on
   semilogx(real(r1save),imag(r1save),'r.')
else
%  Plot unperturbed and perturbed roots
   plot(real(r),imag(r),'kx'), hold on
   plot(real(r1save),imag(r1save),'r.'), hold on,
   
   %r
   %r1save
    for index=1:length(r),
%        cc=0;
%        for cri=1:length(r1save),
%          cri;
%          cr=r1save(cri);
%          mindist=10000;
%          for pri=1:length(r),
%            
%            pr=r(pri);
%            if mindist > abs(cr-pr),
%              
%              mindist=abs(cr-pr);
%            end  
%          end
%          %mindist
%          %abs(cr-r(index))
%          if mindist==abs(cr-r(index)),
%              pri
%              if cc<abs(cr),
%                  cc=abs(cr)
%              end
%          end
%        end
%        if cc>c(index),
%          c(index)=cc
%        end
     cmx=min(real(r(index))*cent+e*c(index)*cirx);
     cmy=min(imag(r(index))*cent+e*c(index)*ciry);
     CMX=max(real(r(index))*cent+e*c(index)*cirx);
     CMY=max(imag(r(index))*cent+e*c(index)*ciry);
     fprintf('minx=%f, maxx=%f, miny=%f, maxy=%f\n',cmx,CMX,cmy,CMY);

     plot(real(r(index))*cent+e*c(index)*cirx,imag(r(index))*cent+e*c(index)*ciry), hold on,
   end
end
if ( miny == 0 & maxy == 0)
   axis([minx,maxx,-1e-300,+1e-300])
elseif (miny == maxy)
   axis([minx,maxx,miny*.99,maxy*1.01])
else
   axis([minx,maxx,miny,maxy])
end
grid
title(['Black xs = roots, red dots = perturbed roots'])
xlabel(['Max relative perturbation to each coefficient = ',num2str(e)])