Finding all roots for a polynomial with either real or complex coefficients.
Let the polynomial be given by:
f(x)=anxn+an-1xn-1....a1x+a0
Where a
n is either a real or complex number. We try to find all roots to this polynomial by solving the equation for x.
anxn+an-1xn-1....a1x+a0=0
There exist many methods available to solve this problem. Among them are Newton’s, Bairstow’s, Halley’s, Laguerre’s, Graeffe’s, Alberth, Durand-Kerner, Jenkins-Traub methods. If you are interested in seeing all these methods in actions go to
http:www.hvks.com/Numerical/winsolve.php
For this Polynomial root finder we use one of the Newton’s variants original developed by Kaj Madsen in the early seventies [K.Madsen: "A root finding algorithm based on Newton Method" Bit 13 1973 page 71-75]. The reason is that it is fast, reliable and can also handle the traditional problems a Newton method can ran into when searching for roots.
Newton has quadratic convergence meaning that for each iteration the numbers of significant digits in the root doubles. What is remarkable by Madsen implementation is that it maintains the quadratic convergence even for multiple roots which otherwise will converge to a root on a much slower pace.
The method works by finding one root at a time, then divide the root up into the polynomial. The resulting polynomial is one degree lower than the original polynomial and you then repeat the process until all roots have been found. For quadratic polynomial or lesser degree we however solve the polynomial directly.
To start a search for a root we go though the following steps.
- The initial start value for the root search. Instead of using a fixed start guess we start with a guess that is less than the modulus of any root of f and as a general strategy try to find the roots in increasing order of modulus to ensure a stable deflation process. [J.H.Wilkinson: "Rounding erros in algebraic processes" Prentice-Hall 1963]
- The iterations phase. Madsen divide the iteration up into two stages. Stage 1 is the phase where we need to get into a position where we know the Newton will converge to a root. This is the most complicate stage and we would need to safeguard the search by being able to successful handle ‘saddle point’ or handle unreasonable steps side etc. Stage 2 is entered when we are sufficient close to the root that we can rely on the standard Newton step to do the job.
- Termination of the iterations. This is another important part if our stopping criteria is too aggressive we will never get there on the other hand if we are to relax the root will not be as accurate as we could possible obtain. For polynomial with real coefficients we use the stopping criteria proposed by D. Adams [D.Adams "A Stopping criterion for polynomial Root finding" Comm ACM 10 1967 pp 655-658]. For complex coefficients polynomial we use an upper bound of the errors in calculating f(x) for the polynomial.
By checking the verbose mode you can get a taste of how this method is working.