Roots of Equations - Numerical Methods in Engineering with MATLAB - page 176

One nice property of Laguerre's methodisthatconverges to aroot, with very few
exceptions, from any starting valueof x .

polyRoots

The function polyRoots in this module computes all the roots of P n ( x ) = 0,

where the polynomial P n ( x ) definedbyits coefficient array a a n + 1 ].

After the first root iscomputedbythe subfunction laguerre , the polynomial is de-
flatedusing deflPoly and the nextzerocomputed by applying laguerre to the
deflatedpolynomial. This process is repeateduntil all n roots have been found.
If a computedroot has avery small imaginary part, it is very likely that it rep-
resents roundoff error. Therefore, polyRoots replaces a tiny imaginary part by
zero.

= [ a 1 , a 2 , a 3 ,...,

function root = polyroots(a,tol)
% Returns all the roots of the polynomial
% a(1)*xˆn + a(2)*xˆ(n-1) + ... + a(n+1).
% USAGE: root = polyroots(a,tol).
% tol = error tolerance (default is 1.0e4*eps).

if nargin == 1; tol = 1.0e-6; end
n = length(a) - 1;
root = zeros(n,1);
fori=1:n
x = laguerre(a,tol);
ifabs(imag(x))<tol;x=real(x);end
root(i) = x;
a = deflpoly(a,x);

end

functionx=laguerre(a,tol)
% Returns a root of the polynomial
% a(1)*xˆn + a(2)*xˆ(n-1) + ... + a(n+1).
x = randn; % Start with random number

n = length(a) - 1;
fori=1:30
[p,dp,ddp] = evalpoly(a,x);
if abs(p) < tol; return; end
g=dp/p;h=g*g-ddp/p;
f = sqrt((n - 1)*(n*h - g*g));