Welcome, guest | Sign In | My Account | Store | Cart

Durand-Kerner method for solving polynomial equations.

Python, 67 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67``` ```# Durand-Kerner method for solving polynomial equations # http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method # http://en.wikipedia.org/wiki/Horner%27s_method # http://en.wikipedia.org/wiki/Properties_of_polynomial_roots # FB - 20130428 import random import math def bound(coefficients): coefficients.reverse() n = len(coefficients) - 1 b = 0.0 for i in range(n): b += abs(coefficients[i] / coefficients[i + 1]) coefficients.reverse() return b def polynomial(z, coefficients): # Horner method t = complex(0, 0) for c in reversed(coefficients): t = t * z + c return t eps = 1e-7 # max error allowed def DurandKerner(coefficients): n = len(coefficients) - 1 roots = [complex(0, 0)] * n bnd = bound(coefficients) retry = True while retry: retry = False # set initial roots as random points within bounding circle for k in range(n): r = bnd * random.random() theta = 2.0 * math.pi * random.random() roots[k] = complex(r * math.cos(theta), r * math.sin(theta)) itCtr = 0 rootsNew = roots[:] flag = True while flag: flag = False for k in range(n): temp = complex(1.0, 0.0) for j in range(n): if j != k: temp *= roots[k] - roots[j] rootsNew[k] = roots[k] - polynomial(roots[k], coefficients) / temp if abs(roots[k] - rootsNew[k]) > eps: # print abs(roots[k] - rootsNew[k]) flag = True if math.isnan(rootsNew[k].real) or math.isnan(rootsNew[k].imag): flag = False retry = True print 'retrying...' break roots = rootsNew[:] itCtr += 1 print "iteration count: " + str(itCtr) return roots # example # x**3-3*x**2+3*x-5=0 coefficients = [complex(-5, 0), complex(3, 0), complex(-3, 0), complex(1, 0)] print "coefficients: " + str(coefficients) print "roots: " + str(DurandKerner(coefficients)) ```

4 comments

FB36 (author) 12 years, 7 months ago

A better way to calculate polynomial values is this:

``````def polynomial(z, coefficients): # Horner scheme
t = complex(0, 0)
for c in coefficients:
t = t * z + c
return t
``````
Petr Kobalicek 12 years, 5 months ago

This is correct version of the Horner scheme for the given algorithm:

``````def polynomial(z, coefficients): # Horner scheme
k = len(coefficients) - 1;
t = coefficients[k]

while k > 0:
k -= 1
t = t * z + coefficients[k]

return t
``````
Shriramana Sharma 11 years, 3 months ago

Hello -- I'm trying to use this code, and testing to find the roots of 2x6+3x5-x4-5x3+8x**2-x+10 :

The correct roots (obtained using SymPy's nroots) are:

``````[-1.66994457011139 - 1.04653556965579*I,
-1.66994457011139 + 1.04653556965579*I,
-0.150681569320409 - 0.883232554208288*I,
-0.150681569320409 + 0.883232554208288*I,
1.0706261394318 - 0.676257117998739*I,
1.0706261394318 + 0.676257117998739*I]
``````

But if I test it with the above code:

``````>>> c=complex
>>> DurandKerner([c(2),c(3),c(-1),c(-5),c(8),c(-1),c(10)])
[(nan+nanj), (nan+nanj), (nan+nanj), (nan+nanj), (nan+nanj), (nan+nanj)]
``````

And with the replaced polynomial function given by Petr, each time it returns a different set of random complex numbers (with one NaN). Please tell me what I am doing wrong? Thank you!

FB36 (author) 11 years, 3 months ago

I tested it and got the same result.

I remember when I had written this code I also had got bad results in some test cases. And I had done some research in the Internet and found Durand-Kerner Method itself is unstable. It does not always converge.

I do not know if there is a known way to force the algorithm towards convergence if it diverges.

Of course I cannot rule out somehow there maybe a mistake in the code but right now I believe there is not. Libraries like SymPy, NumPy obviously using other methods.

One way to resolve this would be to find other Durand-Kerner implementations in the Internet and compare the results.

 Created by FB36 on Mon, 5 Sep 2011 (MIT)

Required Modules

• (none specified)