Durand-Kerner算法(here)的这种实现有什么问题?在
def durand_kerner(poly, start=complex(.4, .9), epsilon=10**-16):#float('-inf')):
roots = []
for e in xrange(poly.degree):
roots.append(start ** e)
while True:
new = []
for i, r in enumerate(roots):
new_r = r - (poly(r))/(reduce(operator.mul, [(r - r_1) for j, r_1 in enumerate(roots) if i != j]))
new.append(new_r)
if all(n == roots[i] or abs(n - roots[i]) < epsilon for i, n in enumerate(new)):
return new
roots = new
当我尝试它时,我必须用keyboardInterrupt
来停止它,因为它不会停止!poly
是{
提前谢谢你, 魔方
编辑:使用numpy多项式需要9次迭代:
^{pr2}$使用pypol多项式它永远不会结束(这可能是pypol中的一个bug):
In [3]: roots.d2(poly1d([1, -3, 3, -5]))
^C---------------------------------------------------------------------------
KeyboardInterrupt
但我找不到虫子!!在
EDIT2:比较__call__
方法与Martin's Poly:
>>> p = Poly(-5, 3, -3, 1)
>>> from pypol import poly1d
>>> p2 = poly1d([1, -3, 3, -5])
>>> for i in xrange(-100000, 100000):
assert p(i) == p2(i)
>>>
>>> for i in xrange(-10000, 10000):
assert p(complex(1, i)) == p2(complex(1, i))
>>> for i in xrange(-10000, 10000):
assert p(complex(i, i)) == p2(complex(i, i))
>>>
EDIT3:如果根不是复数,pypol可以正常工作:
In [1]: p = pypol.funcs.from_roots([4, -2, 443, -11212])
In [2]: durand_kerner(p)
Out[2]: [(4+0j), (443+0j), (-2+0j), (-11212+0j)]
所以它不只是在根是复数的时候才起作用!在
EDIT4:我为numpy多项式编写了一个稍有不同的实现,并且发现在一次迭代之后,(Wikipedia多项式的)根是不同的:
In [4]: d1(numpyp.poly1d([1, -3, 3, -5]))
Out[4]:
[(0.98096328371966801+1.3474626910848715j),
(-0.3352519326012724-0.64406860772816388j),
(2.3542886488816044-0.70339408335670761j)]
In [5]: d2(pypol.poly1d([1, -3, 3, -5]))
Out[5]:
[(0.9809632837196679+1.3474626910848717j),
(-0.33525193260127306-0.64406860772816377j),
(2.3542886488816048-0.70339408335670772j)] ## here
EDIT5:嘿!如果我将if all(n == roots[i] ... )
改为if all(str(n) == str(roots[i]) ... )
它就结束并返回正确的根!!!在
In [9]: p = pypol.poly1d([1, -3, 3, -5])
In [10]: roots.durand_kerner(p)
Out[10]:
[(0.20629947401590029+1.3747296369986026j),
(0.20629947401590013-1.3747296369986026j),
(2.5874010519681994+0j)]
但问题是:为什么它能用不同的复数比较呢??在
更新
现在它起作用了,我做了一些测试:
In [1]: p = pypol.poly1d([1, -3, 3, -1])
In [2]: p
Out[2]: + x^3 - 3x^2 + 3x - 1
In [3]: pypol.roots.cubic(p)
Out[3]: (1.0, 1.0, 1.0)
In [4]: durand_kerner(p)
Out[4]:
((1+0j),
(1.0000002484566535-2.708692281244913e-17j),
(0.99999975147728026+2.9792265510301965e-17j))
In [5]: q = x ** 3 - 1
In [6]: q
Out[6]: + x^3 - 1
In [7]: pypol.roots.cubic(q)
Out[7]: (1.0, (-0.5+0.8660254037844386j), (-0.5-0.8660254037844386j))
In [8]: durand_kerner(q)
Out[8]: ((1+0j), (-0.5-0.8660254037844386j), (-0.5+0.8660254037844386j))
你的算法看起来很好,对我来说在维基百科中的例子是可行的
给予
^{pr2}$关于您的“edit5”:发生这种情况是因为str()没有将数字格式化为完全精度。在
所以不要那样做。在
在任何情况下,代码中的相等性测试:
^{pr2}$是多余的;如果}将为零,所以您可以
n == roots[i]
,那么{再花点功夫找出}会收敛,但您的默认epsilon太小,无法永远停止循环。
epsilon
的默认值是什么;正如我在一篇评论中指出的,求解{1.12e-16
看起来是默认epsilon更好的选择。在为了好玩,试试不适合的Poly([-1,3,-3,1])。。。三个根都等于(1+0j)。。。它需要600多次迭代,最后10次左右迭代中的错误惊人地跳来跳去,直到它在左视场中从很远的地方得到一个合理的解决方案。在
相关问题 更多 >
编程相关推荐