% 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 = vector of polynomial roots
%   e = maximum relative perturbation to make to each coefficient
%   m = number of random polynomials to generate
%
% Output:
%   Plot of perturbed roots
function polyplot(r,e,m)
%
% Generate polynomial coefficients
p=poly(r);
% 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(1), hold off, clf, 
% 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
   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,
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)])