1⁻12→eps © Set complex part limit
newlist(m)→roots © Create list to hold roots
seq(a[i],i,dim(a),1,⁻1)→a © reverse coefficients for simpler looping
a→ad © Copy coefficients for deflation
for j,m,1,⁻1 © Loop to find each root
0→x © Set guess to zero to find smallest remaining root
root(ad,x)→x © Find the root
if gettype(x)="STR":return x © Check for error from root()
if abs(imag(x)) ≤ 2*eps^2*abs(real(x)) © Remove imaginary part of root if very small
real(x)→x
x→roots[j] © Save the root
ad[j+1]→b © Perform forward deflation on remaining coefficients
for jj,j,1,⁻1
ad[jj]→c
b→ad[jj]
x*b+c→b
endfor
left(ad,j)→ad © Keep only deflated coefficients
endfor
for j,1,m © Polish roots with undeflated coefficients
root(a,roots[j])→roots[j]
if gettype(roots[j])="STR":return roots © On error, return roots so far
endfor
return roots
EndFunc
polyroot() will also work for polynomials with complex coefficients.
I will not describe Laguerre's algoithm in detail, but you can find additional description in the references
for Acton and Press below.
polyroot() first finds estimates for all the roots, then the roots are 'polished' to improve accuracy. For
the estimates, polyroot() finds the smallest remaining root, then creates a new polynomial by dividing
that root out of the polynomial. This process is called deflation. The new polynomial has degree of one
less than the original, since one root is divided out. Finding the smallest roots first reduces errors from
round-off; otherwise the small roots might disappear completely. Since deflating the polynomial can
introduce rounding errors by changing the coefficients, the second pass of polishing the roots finds the
roots again, this time using the original, undeflated coefficients, with the roots just found as initial
guesses.
As with any feedback-based iterative method, it is possible for the algorithm to get stuck, oscillating
between two solutions which do not meet the termination criteria. This oscillation is called a limit cycle,
and polyroot() attempts to break limit cycles by periodically perturbing the current solution by some
fractional amount. mt, mr and frac implement the cycle-breaking process such that a fractional step
size from frac is taken every mt steps. As shown, the total number of iterations to find a root is 80. If a
root isn't found in 80 iterations, polyroot() gives up and returns an error message. This would be
unusual but may occur with complex roots. If this error occurs, polyroot() returns a list of the roots
found up to the point of error, and the last element is the string "polyroot too many its". For example, if
the program finds three roots 1, 2 and 3, but fails on the fourth root, the returned list is
{1,2,3,"polyroot too many its"}
If you call polyroot() from another program, you can test for this condition with
root[dim(root)]→r
6 - 121
To Next Page
Loading...
Need help?
General
