# How to find more than one root of a polynomial?

This program finds the first root of the function f, defined in the code. There are 5 roots of this function. (x=1,2,3,4,5) I wish to find all of the roots in this program and print them to the screen. In the main function, you can see that I have attempted to write a while loop to do this(commented out in the code below). But this doesn't work, it just keeps printing the same result to the screen infinitely. Any help would be appreciated.

// This program finds the roots of a polynomial using the half interval search method

#include <iostream> //preprocessor directive needed in order to use std::cout
//and std::cin
#include <iomanip> //preprocessor directive needed in order to use
//a manipulator which uses an argument

using namespace std; //avoids having to uses std:: with cout and cin

//declare the variables
const double a = 0.0, b = 6.0; //limits for f()
double px = 1e-6;              //define the precision
double x = a, dx = 1.0;        //start value and initial interval
double root;

//prototypes for f()
double f (double x);
double halfIntervalSearch (double a, double b, double dx, double px);

int main (int argc, char* argv[])
{
//Inform the user what the program does
cout << "This program finds the roots of\n"
"f(x) = x^5 - 15x^4 + 85x^3 - 225x^2 + 274x -120 = 0\n"
"which is bracketed by [" << a << ", " << b <<
"] using the half interval search method." << endl;

cout << "Searching...\n"; //Lets user know what is happening

root = halfIntervalSearch(a,b,dx,px); //pass the variable

//print out results
//while (x <= b)
{
//x = root + 2*px; //reset value of x
cout << "\nRoot at x = " << setw(10) << setprecision(7) << fixed
<< root << " of value " << f(root) << endl;
}
//allow user to see the results before ending the program
cout << "\n\nPress Enter to end the program";
cin.get();
}

////////////////////////////////////////////////////////////////////////////////
//define f()
double f (double x)
{
return (x*x*x*x*x)-(15*x*x*x*x)+(85*x*x*x)-(225*x*x)+(274*x)-120;
}
///////////////////////////////////////////////////////////////////////////////
//define halfIntervalSearch()
double halfIntervalSearch (double a, double b, double dx, double px)
{
while ((x <= b) && (dx > px))
{
while ( f(x)* f(x+dx) > 0) //while there are no roots
{
//keep track while program finds roots
cout << "x = " << setw(18) << x << " | dx = " << setw(15)
<< right << dx << endl;
x += dx; //step forward
}

dx = dx/2; //half the interval
}
return x;
}

• There are libraries for this in many languages; you could use one of them or check what algorithms they use. One approach is to form the companion matrix and apply an eigenvalue solver. – David Ketcheson Apr 5 '13 at 2:32
• I'm reading your code, and the while loop on condition x <= b doesn't change x or b (the assignment to the former is commented out). What do you really want to happen? – hardmath Apr 5 '13 at 3:47
• @David Ketcheson : I looked up the algorithm that involves computing the eigenvalues of the companion matrix: A = diag(ones(n-1,1),-1); A(1,:) = -c(2:n+1)./c(1); eig(A) this is the matlab verson. (which I can't read) Can you tell me what a companion matrix is, and how to find it using c++ ? Thanks – flamingohats Apr 5 '13 at 11:29
• I just want to point out that companion matrix is usually ill-conditioned and that for high order polynomials (>10), you may run into numerical issues. If you need high accuracy you could use this method, then use the eigenvalues as starting points for something such as Newton's method since you know you are very close to the roots. – Godric Seer Apr 5 '13 at 13:26
• See links to the companion matrix approach/efficiency in my earlier SciComp Answer. In a broad sense a companion matrix is one whose characteristic polynomial is the given polynomial. Often one writes an explicit companion matrix down in upper Hessenberg form making direct use of the polynomial's coefficients, but it's possible to produce a tridiagonal symmetric "companion matrix" from Sturm sequences if the polynomial has only real roots, and from there standard eigenvalue methods will do. – hardmath Apr 5 '13 at 15:09

If you can stand to use complex arithmetic, simultaneous iteration methods might be preferable for computing all the roots of your polynomial. The simplest simultaneous iteration method, the (Weierstrass-)Durand-Kerner method, is effectively equivalent to applying Newton-Raphson to the Vieta relations relating the coefficients and roots of a polynomial, treated as simultaneous equations. (See also this math.SE answer.) Thus, this method is also quadratically convergent, if the polynomial has no repeated roots.

Customarily, the method is started with initial approximations lying around a circle in the complex plane, with the center and radius determined from the coefficients (which is why I mentioned the need for complex arithmetic). If you know in advance that the roots are real, and you already have good approximations for all the roots, one could also use Durand-Kerner for polishing purposes.

A related simultaneous iteration method is Ehrlich-Aberth(-Maehly), mentioned briefly by Gert in his answer. The method has cubic convergence if the roots of your polynomial are all simple, but the iteration formula is a bit more complicated, involving the evaluation of the polynomial's derivative (which of course can be done with e.g. Horner). Picking between these two methods would require some testing on your part.

As pointed out by David Ketcheson, one method is to use the companion matrix and find its eigenvalues (this is what Matlab does for the roots function).

However, if you want to code everything by yourself, you can try to use the Sturm sequences, which you use as a first step to find an interval with only one zero. Then, you can apply one standard methods for finding the root within that interval (Newton, bisection,...).

That is easier to code, but probably less robust than the first approach.

Let's lay out some specs for what the halfIntervalSearch function is supposed to do, and go from there to changes necessary in your code.

I'd assume from the name that it is supposed to find (and return) a root of a (hardcoded) function f using the method of interval bisection. The mathematics behind this is that a continuous function which changes sign in an interval must have at least one root there.

So a typical implementation would be to call the bisection method with such an interval to have the interval refined by successive evaluations at midpoints, with the result that each evaluation results either in finding a (sufficiently approximate) root or in reducing the width of the interval by half (by retaining that half of the earlier interval over which the function is known to change signs).

Your code for halfIntervalSearch does something like this, but I suspect it suffers from being overly indirect in its approach. Of the four arguments passed to this function, the first of them is not used at all. More to the point, instead of a single looping construct which preserves the change-of-sign property of an interval, your implementation using a nested pair of looping constructs whose inner loop is theoretically capable of "stepping forward" over any particular root (possibly all of them): x += dx. I don't think this is a sound approach.

Now the search for roots of polynomials is difficult in complete generality, so we should not set the bar too high. However, going back to the mathematical motivation (Intermediate Value Theorem), we would be better served to implement a routine that iteratively refines an interval with a change-of-signs until a (sufficiently approximate) root is found.

The higher level issue is then how to call/make-use-of such a routine. This is typically done by making a coarsely stepped search of the "big interval" of interest (which in your program is [a,b]) to detect changes-of-sign. Once a (modest size) subinterval where the function changes sign is found, the bisection method routine would be called on just that subinterval (guaranteeing us at least to get one root from that subinterval).

It is a difficult problem to locate all roots of a general polynomial. If a root has even multiplicity, the polynomial will not change signs (except that it becomes zero exactly at that root). More generally you may have a pair of roots that are so close together, the coarsely stepped search passes over both of them in one step, thereby missing these.

Still I recommend you focus on the development of your halfIntervalSearch function in a way that conserves the interval having a change-of-signs, and once that is in place we can perhaps revisit the somewhat vexing topics of roots with multiplicity > 1, etc.

• This is a very clear guide, thank you. I am still working on improving this function so that it works for at least a few polynomials of relatively low powers (perhaps to powers of 6), for different step sizes dx, and for any interval (defined by the user) – flamingohats Apr 9 '13 at 22:01

Another approach would be the deflation approach. In this approach you apply typically an iterative method like Newton's to the original polynomial to find a first root. Next, you apply again Newton's algorithm but this time to the deflated polynomial, i.e. the original polynomial divided by $x-x_1$ where $x_1$ is the first root. By using Horner's rule intelligently you can stably evaluate the deflated polynomial without calculating the coefficients explicitly. You can repeat this until you have fully deflated the polynomial. Typically one would then do one or two extra Newton steps on the original polynomial using the roots just calculated as initial values. I believe this approach is called Maehly's algorithm.

• One should of course be careful with deflation; the wrong order of deflation might have you end with not too accurate roots. See e.g. Peters/Wilkinson, which advocates the use of composite deflation, in an attempt to have the best of both the forward and backward forms. – J. M. May 11 '13 at 18:22