Find a polynomial such that this proposed root finding algorithm fails.

Is this polynomial root finding algorithm below known, and under what conditions for the choice of polynomial coefficients does it find at least one root?

Description of the algorithm:

1. Consider the polynomial:

$$a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0 = 0$$

1. Replace $x^n$ with $y^n$ so that the expression becomes:

$$a_n y^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0 = 0$$

1. Solve for $y$ and you will get:

$$y_1 = \frac{(-1)^{1}(-1)^{\frac{1}{n}}(a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$y_2 = \frac{(-1)^{2}(-1)^{\frac{2}{n}}(a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$y_3 = \frac{(-1)^{3}(-1)^{\frac{3}{n}}(a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$\vdots$$

$$y_n = \frac{(-1)^{n}(-1)^{\frac{n}{n}}(a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

1. Replace $y_1...y_n$ with $x_1...x_n$ and replace all $x$ in the $n$-th equation with $x_n$:

$$x_1 = \frac{(-1)^{1}(-1)^{\frac{1}{n}}(a_{n-1} x_1^{n-1} + \cdots + a_2 x_1^2 + a_1 x_1 + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$x_2 = \frac{(-1)^{2}(-1)^{\frac{2}{n}}(a_{n-1} x_2^{n-1} + \cdots + a_2 x_2^2 + a_1 x_2 + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$x_3 = \frac{(-1)^{3}(-1)^{\frac{3}{n}}(a_{n-1} x_3^{n-1} + \cdots + a_2 x_3^2 + a_1 x_3 + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

$$\vdots$$

$$x_n = \frac{(-1)^{n}(-1)^{\frac{n}{n}}(a_{n-1} x_n^{n-1} + \cdots + a_2 x_n^2 + a_1 x_n + a_0)^{\frac{1}{n}}}{(a_n)^{\frac{1}{n}}}$$

1. Use $x_n=0$ as the seed point and iterate each equation in 4. infinitely many times. In the program below they are iterated 4000 times.

It then appears that when iterating the equations in 4. at least one of them will converge to a zero.

(*start*)
(*Mathematica 8*)
Clear[x, polynomial, X, y, polynomialdegree];
polynomialdegree = 7;
ycoefficient = RandomInteger[{-4, 4}]; polynomial =
Sum[RandomInteger[{-4, 4}]*x^n, {n, 0, polynomialdegree - 1}] +
If[ycoefficient == 0, -1, ycoefficient]*y^polynomialdegree;
Print["The polynomial is: ", polynomial];
nn = 4000;
Do[Clear[x, X, y]; X = y /. Solve[polynomial == 0, y][[i]];
x = 0;
Table[x = N[Round[X, 10^-30], 30], {n, 1, nn}];
y = x;
Print["x = ", x, " polynomial value at x is: ", polynomial];, {i, 1,
polynomialdegree}]
(*end*)


What I am really asking is to find a polynomial such that the proposed root finding algorithm fails.

Edit 10.9.2016:

This version of the algorithm works for complex numbers:

(*start*)
(*Mathematica 8*)
Clear[x, polynomial, X, y, polynomialdegree];
polynomialdegree = 5;
ycoefficient = -Rationalize[RandomComplex[{-10 - 10*I, 10 + 10*I}],
0]; polynomial =
Sum[Rationalize[RandomComplex[{-10 - 10*I, 10 + 10*I}], 0]*x^n, {n,
0, polynomialdegree - 1}] +
If[ycoefficient == 0, -1, ycoefficient]*y^polynomialdegree;
Print["The polynomial is: ", polynomial];
nn = 8000;
Do[Clear[x, X, y];
X = y /. Solve[polynomial == 0, y, WorkingPrecision -> 100][[i]];
x = 10;
Table[x = N[Round[X, 10^-30], 30], {n, 1, nn}];
y = x;
Print["x", i, " = ", x, " polynomial value at x", i, " is: ",
polynomial];, {i, 1, polynomialdegree}]
(*end*)

• If I'm understanding your construction correctly, we write $a_nx^n+f(x)$ for the polynomial, $x_0$ is a seed point and $x_i=\sqrt[n]{\frac{f(x_{i-1})}{a_n}}$. There is no guarantee that this will converge to a root and I can come up with examples where it does not. – Michael Burr May 22 '16 at 18:55
• @MichaelBurr Yes that is right, I added the seed point now to the description. The alternatives for $x_i$ are of course as many as the degree of the polynomial. Ok I see, then the method is not certain to find a root. – Mats Granvik May 22 '16 at 18:59
• If you could give an example polynomial when it cannot find at least one root, I would be happy. – Mats Granvik May 22 '16 at 19:01
• What do you do about multiple roots in the $y^n$ step? Consider trying to solve a quadratic equation both of whose roots are negative. – Rahul May 22 '16 at 19:23
• My formula above is in error, I think that it should be $x_i=\sqrt[n]{-\frac{f(x_{i-1})}{a_n}}$. – Michael Burr May 22 '16 at 19:25

I appears that the algorithm does not fail for the Wilkinson polynomial: $$w(x) = \prod_{i=1}^{20} (x - i) = (x-1)(x-2) \cdots (x-20)$$

but already for a degree $10$ Wilkinson like polynomial the computational computational cost becomes too high in order for one to be able to verify that the algorithm converges.

The Mathematica program below shows that the algorithm is not able to easily find any roots of the polynomial $p(x)$: $$p(x) = \prod_{i=1}^{10} (x - i) = (x-1)(x-2) \cdots (x-10)$$

(*start*)
(*Mathematica 8*)
Clear[x, polynomial, X, y, polynomialdegree];
polynomialdegree = 10;
polynomial =
Expand[Product[x - n, {n, 1, polynomialdegree}]] -
x^polynomialdegree + y^polynomialdegree;
Print["The polynomial is: ", polynomial];
Print["Number of iterations used for each iteration path:"]
nn = 400000
Do[Clear[x, X, y];
X = y /. Solve[polynomial == 0, y, WorkingPrecision -> 100][[i]];
x = 0;
Table[x = N[Round[X, 10^-30], 30], {n, 1, nn}];
y = x;
Print["x", i, " = ", x, " polynomial value at x", i, " is: ",
polynomial];, {i, 1, polynomialdegree}]
(*end*)


In the program the seed point is $x=0$. It appears that whenever the seed point is not a root of the polynomial then the algorithm becomes very slow.

Setting the seed point to and integer between $1$ and the degree of the Wilkinson like polynomial then the algorithm succeeds in a trivial way which of course does not count as success for the algorithm.

400000 iterations where used in the program and it appears that the algorithm finds an approximate root of the degree $10$ Wilkinson polynomial as:

$$x_2=5.99999552884940939567335401803+\frac{2.76805655896054526929533 i}{10^6}$$

So far I have not been able to find a polynomial for which the algorithm fails.

However if the convergence is this slow for the Wilkinson polynomial then dividing out the roots as they are found is not a good idea as pointed out here at Mathoverflow.