% 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)])