% Matlab code polyperturb.m 
% function ebnd = polyperturb(r,e,m)
%
% 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:
%   ebnd = error bounds for each root (to be supplied by students!)
%   Plot of perturbed roots
%
function ebnd = polyperturb(r,e,m)
%
% Generate polynomial coefficients
p=poly(r);
%
% Generate m random polynomials and compute their roots
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];
end
%
%  Compute min and max roots to determine axes for plotting
   minx = min(min([real(r1save);real(r')]));
   maxx = max(max([real(r1save);real(r')]));
   miny = min(min([imag(r1save);imag(r')]));
   maxy = max(max([imag(r1save);imag(r')]));
%
%  If all roots lie in the right halfplane, and
%  vary greatly in magnitude, use a semilog plot
figure(1), hold off, clf, 
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
%
% Adjust axes on plot to make it look nice
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)])
ebnd = ones(size(r))./0;