Durand-Kerner method for solving polynomial equations.
Python, 67 lines
Download
Copy to clipboard1 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))
|
Tags: math, mathematics
A better way to calculate polynomial values is this:
This is correct version of the Horner scheme for the given algorithm:
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:
But if I test it with the above code:
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!
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.