On 3/19/07, sen1 at math.msu.edu <sen1 at math.msu.edu> wrote:
> There have been several recent developments on the "root finding
> problem" (variations of Newton's method for approximate roots).
...
> Perhaps it would be interesting for someone to implement some of the
> new methods in maxima, say the Shub-Smale method or the
> Hubbard-Schleicher-Sutherland method.
One thing worth improving would be complex rootfinding for univariate
polynomials: allroots is pretty bad. Try the following, for example:
poly: expand( product( (x - (1+i/10)), i, 1,20) )$
realroots(poly) gets all the roots within 1.8e-7 (with rootsepsilon =
1.0e-7 -- so not quite within error bounds), which is excellent.
On the other hand, allroots gives poor results (recall that the input
was using *exact* coefficients, not floats): 6 of the 20 roots have
significant imaginary parts (> 0.05) and almost all are far from the
true roots. Here's a list of the root to which each returned root is
closest and its distance from it:
[[ 1, 0.000], [ 2, 0.004], [ 3, 0.014], [ 4, 0.004],
[ 5, 0.111], [ 5, 0.111], [ 7, 0.046], [ 7, 0.239],
[ 7, 0.239], [10, 0.339], [10, 0.339], [13, 0.366],
[13, 0.366], [16, 0.021], [16, 0.310], [16, 0.310],
[18, 0.196], [18, 0.196], [20, 0.055], [20, 0.055]]
The median error is 0.15 and the maximum error is 0.37! And as far as
I can tell there is no way of requesting higher precision.
The basic problem appears to be that allroots works in floating-point,
not in rats or bigfloats. It turns out that realroots(float(poly))
finds only 4 roots (with ratepsilon=1.0e-17). That is, simply
converting the coefficients to floats already perturbs the roots
significantly, as can be seen by:
makelist(length(float(map(rhs,realroots(ev(rat(bfloat(poly)),fpprec:k,ratepsilon:0.1^k))))),k,2,25);
=> [2, 2, 2, 4, 2, 2, 0, 4, 2, 2, 2, 2, 2, 2, 6, 6, 10, 12, 20,
20, 20, 20, 20, 20]
It takes 20 digits to have 20 roots, and 22 digits to have them within 5e-5.
Of course, in this particular case, factor finds the roots exactly,
but that's neither here nor there....
-s