Newton Maehly

Roots of a polynomial

Basically

The basic idea to find a root of a function is to start at a particular point and run along the graph of the function till it crosses the cero line. In a computer program this can be done by starting at point P₁ and laying a tangent onto the graph. This tangent should cross the cero line somewhere at x₂ and there where it crosses it is already a bit closer to the cero crossing of the graph than where we started. From x₂ we go to P₂ and do the same again: Lay another tangent onto the graph and go to its cero crossing x₃ and do the same again till we are close enough at the cero crossing. Close enough means the difference between one such iteration and the next is small enough.

In mathematics that means with

We get

and so on...

The method of Newton to find the roots of a polynomial basically does this by running a kind of polynomial evaluation by Horner for f(x) and f’(x) in a combination.

For the part f(x) it starts by dividing the polynomial

by (x - p) whit p <> 0 and bringing it to the form

Expanding this and doing a parameter comparison leads to

a₀ = b₀, a₁ = b₁ – p*b₀, …, a_n = b_n – p*b_n-1

This solved for b gives an algorithm to calculate the b parameters in a loop:

b₀ = a₀, b_j = a_j + p*b_j-1 for j = 1,2,...,n

This might not look too interesting, but if we set x = p and look further above, we see that the polynomial becomes equal bn in this case. That means if x = p we just have to run through the loop and get the value of the polynomial at the end in b_n. O.k. that’s still not that interesting :-)

If we have a look about the differentiation of f(x), things might become somewhat more interesting.

For the part f’(x) we first substitute

And get

with R_n = b_n

partially differentiated

And with x = p

To evaluate this we divide it by x – p again and get

With

This gives the algorithm to calculate the c parameters from the b parameters from above in a loop:

c₀ = b₀, c_j = b_j + pc_j-1 for j = 1,2,...,n-1

and at the end of the loop we get P’(x) = c_n-1.

Now the loop for the b parameters from above and the loop for the c parameters can be put together.

The first step is

b₀ = a₀ and c₀ = b₀

With this we start a loop looking like

a₁ = b₁ – p*b₀ and c₁ = b₁ + pc₀

a₂ = b₂ – p*b₁ and c₂ = b₂ + pc₁

…

a_i= b_i – p*b_i-1 and c_i = b_i + pc_i-1

For the c part we have to skip the last cycle of the loop. In C# that looks like

bn = a[0];
cn = bn;
for (i = 1; i < order; i++)
{
      bn = xk * bn + a[i];
      if (i < order - 1)
            cn = xk * cn + bn;
}

Whit a[0] as the first supporting point of the polynomial which is defined in the float array[1..order] a.

At the end of the first loop we get the next point x_k by subtracting bn/cn from the starting point. At the end of the following loops we get the next x_k by subtracting bn/cn from the last x_k.

In my sample for a polynomial of the 5. degree with a₀ = 1.0, a₁ = -4.0, a₂ = 0.0, a₃ = 10.0, a₄ = -1.0 and a₅ = -6.0 I start at 5.0 and get x_k = 4.29411 after one cycle, x_k = 3.765 after the second and x_k = 3.00… as one first root of the polynomial.

Now we got one root, but there are more than just one. To find the others it would be a possible solution, to divide the starting polynomial by (x-x_k) and start the root finding algorithm again and do this for each root we find. That’s called deflation.
As xk is not the exact root of the polynomial, the division by (x-x_k) brings a inaccuracy into the polynomial that increases for each root we find. That’s not the best solution.

Therefore Maehly modified the algorithm a bit.

To understand his idea we have a look at a polynomial that has 3 roots. With the roots z₁, z₂, z₃ this looks like

P(x) = (x – z₁) * (x – z₂) * (x – z₃)

And the differentiation of it is

P’(x) = (x – z₂) * (x – z₃) + (x – z₁) * (x – z₃) + (x – z₁) * (x – z₂)

The Division P‘(x)/ P(x) becomes

And separated

That means the formula from above

Can be written a

And we can carry out an implicit deflation for each root by modifying the formula for x_k a

For n found roots.

That means in the algorithm we have to distinguish between “now root found jet” and “found at least one root” and use the unmodified formula for the first case and the modified for the second case. That looks like

if (n == 0)
{
      if (cn != 0) // no root found jet => go the the next P
      {
            dx = bn / cn;
            xk = xk - dx;
      }
      else
            if(Math.Abs(bn) < 1E-8) // that's already a root point
                  dx = 0.0;
}
else
{
      if (bn != 0) // one root found => keep it and subtract it from the polynomial
      {
            xz = 0;
            for (j = 0; j < n; j++)
                  xz = xz + 1.0 / (xk - x[j]);
            dx = 1.0 / (cn / bn - xz);
            xk = xk - dx;
      }
      else
            xk = xk + 0.5; // we started right at one root point just move a bit
}

With n holding the number of found roots and x[j] the found roots themselves.

Now we just have to take care that our algorithm converges and does not run mad. This can be done by checking for delta(x_k) <= delta(x_k-1) and delta(x_k) <> 0 and x_k does not become huge.

Another point is the case that c_n or b_n is 0:

If we don’t have at least one root found and c_n and b_n both are 0, then our starting point is a first root already and we can take this one right away.

If we have one or more roots already found and b_n is 0, we don’t use that as the next x_k will be calculated including all the found roots. In this case we just have to move a little bit by increasing (or decreasing) x_k to get away from this point. I do this by adding 0.5 to x_k arbitrary.

To avoid running into roots we already have found, I move the starting point after each found root to the mean of the last starting point and the newly found root.

In this manner the root finding algorithm by Newton-Meahly works quite fine. It found all the roots of different sample polynomials of the 5. and the 4. degree :-)

I implemented a little online solver in JavaqScript. That's here

C# Demo Project Newton Maehly

Newton_Maehly.zip