Rootfinding (of polynomial in one variable) is one of the classic problems
of numerical analysis.
Commercial macsyma has a few extra contributions over maxima to solutions to
this problem;
One of them is to do a root "refinement" using newton iteration.
The main issue there is that if there are 10 roots, crudely computed,
refining them might not find 10 DIFFERENT roots. The refinements can wander
off.
Because of the fractal nature of newton iteration it is tricky to prevent
this, but I think there are ways by "deflating" to fix this, but deflation
has its problems too.
I'm not an expert on this, but I've tried to get students to work on it.
Here's a citation to what I think is a good paper
http://citeseer.ist.psu.edu/469859.html
> -----Original Message-----
> From: maxima-bounces at math.utexas.edu
> [mailto:maxima-bounces at math.utexas.edu] On Behalf Of Stavros Macrakis
> Sent: Tuesday, March 20, 2007 9:23 AM
> To: sen1 at math.msu.edu
> Cc: maxima list; Alasdair McAndrew; Robert Dodier
> Subject: Re: [Maxima] Solving (systems of) non-linear equations?
>
> 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)),f
> pprec: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
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>