[SOLVED] numpy polynomial roots incorrect compared to Matlab?

numpy polynomial roots incorrect compared to Matlab?

I would like to figure out the roots of the polynomial x^2 - 800x + 160000, which has a multiple root of 400.

In Matlab, if I type roots([1 -800 160000), I correctly get ans = 400.0000, 400.0000 whereas numpy np.roots([1, -800, 160000]) gives me array([400.+2.69940682e-06j, 400.-2.69940682e-06j]).

Why does numpy show some sort of floating point precision issue when matlab does not? Do they not both use an algorithm based on the companion matrix? How can I reliably determine, in Python, that the "actual" roots are 400, as Matlab seems to be doing?

Edit: This is a "toy" example that I am facing in a larger, more complex algorithm that is not necessarily faced with only real roots and multiple roots i.e. I do not actually know in advance, in my more complicated case, that the root is multiple and real.

Solution

I run this slightly ugly test:

max([
    np.abs(
        (np.roots(
            [1, -2 * i, i ** 2]
        ) / i - 1)
    ).max() 
    for i in range(1, 1000000, 10)
])

and got result 1.91e-08

You can do the same for your version of numpy and estimate the precision of answers, and if a ratio of the imaginary part to the real is less than it you can just drop it.

Also, interesting fact from tests: if roots are different, errors are very little, ~1e-16 or even 0, but when roots are the same, errors are almost never zero and have magnitude of 1e-8~1e-9.