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) 10 years, 2 months ago  # | flag

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 10 years ago  # | flag

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 8 years, 11 months ago  # | flag

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) 8 years, 11 months ago  # | flag

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)
Python recipes (4591)
FB36's recipes (148)

Required Modules

  • (none specified)

Other Information and Tasks