I have an univariate polynomial of some degree – how do I numerically find all of its real roots?
I never thought I would ask this question – everyone knows how to find polynomial roots, right..? Please, can someone explain to me how this works? This is what I think I understand so far:
Newton method – it can find only one root at a time (if it converges) so in order to list all the roots I need to somehow guess where to start from. How to guess?
Bisection – also finds only one root and it needs two endpoints with opposite signs. How to guess those endpoints? What if the function simply does not have opposite signs, for example $f(x) = x*x$?
Do I need to know at least the interval I want to search in? Say, that I'm only interested in roots from -1000 to 1000? Or is it possible to list all the roots without any initial guesses?
Then I found about http://en.wikipedia.org/wiki/Sturm%27s_theorem but I don't really understand it yet – is it somehow helpful in this context?
Is the approach basically combination of the above? Count the signs with Descartes (what for?), find the intervals with Sturm, check if there are opposite signs, if no – find the extreme, if yes – bisect?
And I also read somewhere that it is possible to covert the polynomial into a system of linear equations (matrix?) and solve that? Is this possible?
(Note: I search this forum and found this but I don't really understand any of it – "Homotopy Continuation"?)
Best Answer
You can use Sturm's Theorem to find initial guesses for Newton's method, along the lines of Johannes's comment.
The matrix formulation you're referring to is the companion matrix of the polynomial. Roots of the polynomial become eigenvalues of the companion matrix, so you can then use any eigenvalue algorithm to find the roots.
However, in my experience, if you need to find all of the roots of an arbitrary polynomial, the Jenkins-Traub algorithm is easily the best, both in terms of efficiency and robustness, and it's what commercial software like Mathematica uses under the hood when you ask it to numerical solve for polynomial roots. The one downside is that it is far from trivial to implement yourself; fortunately it's not hard to find codes on the Internet, in both Fortran and C.