The Weierstrass-Durand-Kerner method
This is a method for finding all the roots of a polynomial simultaneously, by applying a sort of Newton-Raphson method. It gets its name from everybody associated wth it: Weierstrass published a version of the algorithm in 1891, and then rediscovered (independently) by Durand in 1960 and Kerner in 1966, as you can see at its Wikipedia page.
It's easier to describe by an example:
\[ x^4 - 26x^2 - 75x -56 = 0 \]
which has roots \((5\pm\sqrt{57})/2, (-5\pm\sqrt{3} i)/2\). Numerically, the roots are
\[ 6.27491721763537,\; -1.27491721763537,\; -2.5 + 0.866025403784439i,\; -2.5 - 0.866025403784439i \]
For a fourth degree equation such as above, if one approximation to the roots is given by \(x_0, x_1,x_2,x_3\), then the next approximation is given by
\[\begin{align*} y_0 &= x_0 - \frac{f(x_0)}{(x_0-x_1)(x_0-x_2)(x_0-x_3)}\\ y_1 &= x_1 - \frac{f(x_1)}{(x_1-x_2)(x_1-x_3)(x_1-x_0)}\\ y_2 &= x_2 - \frac{f(x_2)}{(x_2-x_3)(x_2-x_0)(x_2-x_1)}\\ y_3 &= x_3 - \frac{f(x_3)}{(x_3-x_0)(x_3-x_1)(x_3-x_2)} \end{align*}\]
Suppose however we create the temporary polynomial
\[ t(x) = (x-x_0)(x-x_1)(x-x_2)(x-x_3). \]
Since
\[\begin{align*} t'(x) &= (x-x_1)(x-x_2)(x-x_3) + (x-x_0)(x-x_2)(x-x_3) \\ & \qquad + (x-x_0)(x-x_1)(x-x_3) + (x-x_0)(x-x_1)(x-x_2) \end{align*}\]
we can write the above expressions for \(y_k\) more simply as
\[ y_k = x_k - \frac{f(x_k)}{t'(x_k)}. \]
which shows the connection with the Newton-Raphson method. Here's how it might be done in PARI/GP - more specifically, in the scripting language gp:
durandkerner(f,eps,N = 200) =
{
local(pd,xs,ys,dsa,dss,count,dt);
pd = poldegree(f);
xs = vector(pd,k,(0.6 + 0.8*I)^k);
ys = vector(pd);
ds = vector(pd);
dsa = 1.0;
dss = Vec([1.0]);
count = 1;
while(dsa > eps && count < N,
dt = deriv(prod(k = 1,pd,x - xs[k]),x);
for(k = 1,pd,
ys[k] = xs[k] - subst(f,'x,xs[k])/subst(dt,'x,xs[k]);
ds[k] = abs(ys[k] - xs[k]);
xs[k] = ys[k]
);
dsa = vecsum(ds)/pd;
dss = concat(dss,[dsa]);
count = count + 1;
);
return([count,dss,xs]);
}
We've chosen as our beginning values \((0.6 + 0.8i)^k\) - one of the requirements of the method is that all computations must be in complex numbers.
Here goes:
\p20
f = x^4 - 26*x^2 - 75*x - 56
s = durandkerner(f,0.1^20);
print("number of iterations: ",s[1]);
print()
vp = vecextract(s[2],"-10..-1"); \\ extracts the last ten differences
{for(k = 1,10,
printf("%.10g\n",vp[k]))} \\ and prints them
print()
{ for(k = 1,poldegree(f), \\ prints out the solutions
printf("%.15g\n",s[3][k]) ) }
which produces the output:
number of iterations: 17
1.123079176
0.6846519849
0.3550410771
0.1757081770
0.04978125947
0.003740061706
2.016207001 e-5
5.729914134 e-10
4.478174636 e-19
2.615741385 e-37
6.27491721763537-1.45113252725757 e-89I
-2.50000000000000+0.866025403784439I
-1.27491721763537+3.98066040407694 e-74I
-2.50000000000000-0.866025403784439I
As you see, the average absolute difference between successive iterates is about \(3.7\times 10^{-37}\), obtained in 17 iterations. And the values are certainly correct. You'll also notice that the imaginary parts of the first and third solutions are vanishingly small; in effect, zero.
For a bit of fun, let's try with a high precision:
\p1000
f = x^4 - 26*x^2 - 75*x - 56
s = durandkerner(f,0.1^1000);
the result of which (using the same printing script as above) is
number of iterations: 22
0.003740061706
2.016207001 e-5
5.729914134 e-10
4.478174636 e-19
2.615741385 e-37
8.409722008 e-74
8.032655807 e-147
6.583139811 e-293
3.794211385 e-585
9.907223519 e-1170
6.27491721763537+2.96119399627067 e-2361I
-2.50000000000000+0.866025403784439I
-1.27491721763537+1.01326544935096 e-2339I
-2.50000000000000-0.866025403784439I
It's only taken 5 more steps to get over 1000 place precision! This is because, once it hits its stride, so to speak, this method converges quadratically, which you can see in that the accuracy of the iteration bascally doubles each step.
We can also experiment with a large polynomial:
deg = 50;
f = x^deg + sum(k=0,deg-1,(random(20)-10)*x^k);
s = durandkerner(f,0.1^1000);
print("number of iterations: ",s[1]);
print()
vp = vecextract(s[2],"-10..-1"); \\ extracts the last ten differences
{for(k = 1,10,
printf("%.10g\n",vp[k]))}
which produces the output
number of iterations: 56
4.520878142 e-6
7.817132307 e-9
9.072720843 e-14
1.471330591 e-23
3.873798136 e-43
2.685287774 e-82
1.290323278 e-160
2.979297878 e-317
1.588344551 e-630
4.514465082 e-1257
(If you try this you'll get different numbers, seeing as you'll be sarting with a different random polynomial.) We won't print out all the solutions, as there are too many, but we can print out the absolute values \(\|f(x_k)\|\) of the function at each solution:
{
for(k=1,deg,
printf("%.10g\n",abs(subst(f,'x,s[3][k])))
)
}
which provides a list of 50 valaus beginning (in my case) with:
2.395686462 e-2332
9.801085659 e-2386
2.362011441 e-2291
4.227497886 e-2284
4.335737935 e-2330
These are vanishingly small; in effect zero.
Playing with the algorithm
To try this algorithms, open up the PARI/GP site, and go to "Try GP in your browser", or just go here. Then you can cut and paste the Durand-Kerner function into a cell and play with it.
For ease of playing, here are a couple of "helper" functions for printing the
output. The function printcomplex simply adds more space around the + or - in
a complex number, and the function printoutput does just that.
printcomplex(z) =
{
local(rz,iz);
rz = real(z); iz = imag(z);
if(imag(z) < 0, printf("%.10g - %.10gI\n",rz,-iz),
printf("%.10g + %.10gI\n",rz,iz));
}
printoutput(s) =
{
print("number of iterations: ",s[1]);
print();
if(length(s[2]) < 11, foreach(s[2],z,printf("%.10g\n",z)),
foreach(vecextract(s[2],"-10..-1"),z,printf("%.10g\n",z)));
print();
foreach(s[3],z,printcomplex(z));
}
So, open up GP "in your browser", and in a cell enter the durandkerner,
printcomplex, and printoutput functions. You can enlarge the cell from the
bottom right. Press the "Evaluate with PARI" function, or simply Shift+Enter.
Open up a new cell, and enter something like
\p50
f = x^4 + x^3 - 19
sol = durandkerner(f,0.1^50);
printoutput(sol)
Enjoy!