getType(r)="STR"
polyroot() uses Laguerre's algorithm for both the initial root-finding and the polishing, however, this is
not mandatory: you could use some other method for the polishing. The references at the end of this
tip discuss this more and give some suggestions.
Since polyroot() uses complex arithmetic even to find the real roots, we have a dilemma I have not
encountered before: when do we decide that a complex root is actually real, and only has a small
complex part because of round-off error? It may be that the root really does have a small complex part,
in that case we don't want to throw it away. Conversely, we don't really want to return a small complex
component when the root is real. Press and his co-authors use this criteria, for a complex root of x = a
+ b*i:
if b
[
2
$
eps
2
$
a then a
d
x
eps is the desired relative accuracy of f(x). Press gives no justification for this criteria, and I could not
derive it from several reasonable assumptions. So, I posted this as a question on the Usenet
newsgroup sci.math.num-analysis. Mr. Peter Spellucci was kind enough to answer with this:
"I think this is one of the typical ad hoc "equal to zero within the roundoff level" decisions
one often is forced to use. First of all, this makes sense only for a real polynomial, and I
assume that this decision is applied simultaneously for the conjugate complex value. In
Laguerre's method there appears the following:
denom= p'+sign(p')*sqrt((n-1)*((n-1)*p'^2-n*p*p'')) )
and if the radicand is negative, it branches into the complex plane. Hence a sound decision
would be based on the possible roundoff level in the evaluation of this expression, in the
first position on the possible roundoff errors in the evaluation of p, p' and p''. Bounds on this
can be computed within the evaluation for example using Wilkinson's techniques. But this
bounds may be too pessimistic and hence one usually relies on some rules of thumb. In the
real case with a complex zero z, also conj(z) will be a zero. Multiplying out we get a
quadratic factor
x^2-2*re(z)*x+abs(z)^2
and if this factor would not change within the computing precision by neglecting the
imaginary part I would set it to zero. This however would mean |b| < |a|*sqrt(eps) , if eps
represents the computing precision. Hence I wonder a bit about the settings you mention."
There is quite a difference between Peter's criteria and that used in NRIF: compare 1.4E-6 to 2E-24, if
we use eps = 2E-12. Peter's criteria, while theoretically quite sound, seems rather extreme in throwing
away complex roots. At this point, I discovered that the authors of Numerical Recipes run a forum for
this type of question, so I posted:
In the xroots() code on p367 of NRIF, a complex root x = a + bi is forced real if |b| <= 2*eps**2*|a|.
I cannot derive this from a few different assumptions; the closest I can come is b^2 <= 2*eps*|a|,
assuming that the root is forced real if |x| - |a| <= eps. How is this criteria derived?
One of the authors of Numerical Recipes, Saul Teukolsky, answered:
The criterion in the Recipe, using eps**2, is purely empirical. After all, you may well have a root
that has a very small imaginary part that is meaningful. But your criterion, with eps, is perfectly
6 - 122
To Next Page
Loading...
Need help?
General
