Finding the roots of a simple polynomial with mpmath

Question

I'm trying to use mpmath.polyroots to find the roots of a simple polynomial with integer coefficients x*(x-4)**3, which when expanded has the coefficient vector of [1, -12, 48, 64, 0]. The following code fails:

import mpmath
p = [  1, -12,  48, -64,   0]
print mpmath.polyroots(p,maxsteps=2000)

with the error:

Traceback (most recent call last):
  File "poly.py", line 3, in <module>
    print mpmath.polyroots(p,maxsteps=2000)
  File "/usr/local/lib/python2.7/dist-packages/mpmath/calculus/polynomials.py", line 188, in polyroots
    % maxsteps)
mpmath.libmp.libhyper.NoConvergence: Didn't converge in maxsteps=2000 steps.

Increasing the number of steps does not help. The expected answer is obviously [0,4,4,4].

Is mpmath unable to find the roots of a polynomial if there is some multiplicity? How can I resolve this?

A triple root leads to an error level of `1e-5` in a numerical solver due to perturbations in the floating point evaluation. With a typical default level of `1e-6` or smaller this may be impossible to resolve. Cluster analysis to find multiple or close by roots might prevent this. — Lutz Lehmann, Aug 22 '15 at 15:21
@LutzL I think about the word cluster analysis in probably a different light (e.g. my use case for it would be something like principal component analysis). How does cluster analysis relate to finding multiple roots? — Hooked, Aug 22 '15 at 18:32
Multiple roots are closely related in numerical behavior to clusters of roots. I do not know where it is done in software, but from a mathematical perspective it would be advantageous to identify relatively small disks in the complex plane that contain several or multiple roots and find the corresponding polynomial factors. This intermediate factorization does not suffer from the loss of precision as in the case of a root with multiplicity. However, the small-degree factors might still lead to unstable roots like in the example. — Lutz Lehmann, Aug 23 '15 at 18:23

score 4 · Accepted Answer · edited Oct 08 '20 at 21:08

The documentation mentions that increasing extraprec may be necessary to achieve convergence. This is the case for your example, due to root multiplicity.

>>> mpmath.polyroots(p, maxsteps=100, extraprec=110)
[mpf('0.0'), mpf('4.0'), mpf('4.0'), mpf('4.0')]

The extraprec parameter means the extra number of digits to be used in the process of calculation, compared the number of digits needed in the result (which is 15 by default, set globally in mpmath.mp.dps). Its default value is 10; increasing it to 110 achieves convergence in 100 steps. This actually takes less time than trying (and failing) to find the roots with maxsteps=2000 at the default extraprec level.

Finding the roots of a simple polynomial with mpmath

1 Answers1