2

この 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の多項式インスタンスです。pypol

事前にありがとう、ルービック

編集:numpy多項式を使用すると、9回の反復が必要です:

In [1]: import numpy as np

In [2]: roots.d1(np.poly1d([1, -3, 3, -5]))
3
[(1.3607734793516519+2.0222302921553128j), (-1.3982133295376746-0.69356635962504309j), (3.0374398501860234-1.3286639325302696j)]
[(0.98096328371966801+1.3474626910848715j), (-0.3352519326012724-0.64406860772816388j), (2.3542886488816044-0.70339408335670761j)]
[(0.31718054925650596+0.93649454851955749j), (0.49001572078718736-0.9661410790307261j), (2.1928037299563066+0.029646530511168612j)]
[(0.20901563897345796+1.5727420147652911j), (0.041206038662691125-1.5275192097633465j), (2.7497783223638508-0.045222805001944255j)]
[(0.21297050700971876+1.3948274731404162j), (0.18467846583682396-1.3845653821841168j), (2.6023510271534573-0.010262090956299326j)]
[(0.20653075193800668+1.374878742771485j), (0.20600107336130213-1.3746529207714699j), (2.5874681747006911-0.00022582200001499547j)]
[(0.20629950692533283+1.3747296033941407j), (0.20629947661265013-1.374729584400741j), (2.5874010164620169-1.899339978055233e-08j)]
[(0.20629947401589896+1.3747296369986031j), (0.20629947401590082-1.3747296369986042j), (2.5874010519682002+9.1830687539942581e-16j)]
[(0.20629947401590029+1.3747296369986026j), (0.20629947401590026-1.3747296369986026j), (2.5874010519681994+1.1832913578315177e-30j)]
Out[2]: 
[(0.20629947401590029+1.3747296369986026j),
 (0.20629947401590029-1.3747296369986026j),
 (2.5874010519681994+0j)]

pypol 多項式を使用しても終了しません (おそらく pypol のバグです):

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 多項式のわずかに異なる実装を作成し、1 回の反復後に (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))
4

2 に答える 2

1

あなたのアルゴリズムはうまく見えます。ウィキペディアの例ではうまくいきます

import operator
class Poly:
    def __init__(self, *koeff):
        self.koeff = koeff
        self.degree = len(koeff)-1

    def __call__(self, val):
        res = 0
        x = 1
        for k in self.koeff:
            res += x*k
            x *= val
        return res

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

print durand_kerner(Poly(-5,3,-3,1))

与える

[(0.20629947401590026+1.3747296369986026j), 
 (0.20629947401590026-1.3747296369986026j), 
 (2.5874010519681994+8.6361685550944446e-78j)]
于 2010-10-30T08:48:06.567 に答える
1

「EDIT 5」について：これは、 str() が数値を完全な精度でフォーマットしないために発生します。

>>> print str((2.5874010519681994+8.636168555094445e-78j))
(2.58740105197+8.63616855509e-78j)
>>> print repr((2.5874010519681994+8.636168555094445e-78j))
(2.5874010519681994+8.636168555094445e-78j)
>>>

だから、それをしないでください。

いずれにせよ、コード内の等価テスト:

if all(n == roots[i] or abs(n - roots[i]) < epsilon for i, n in enumerate(new)):

冗長です。の場合n == roots[i]、abs(n - roots[i])ゼロになるので、単に行うことができます

if all(abs(n - roots[i]) < epsilon for i, n in enumerate(new)):

epsilonそして、デフォルトがどうあるべきかを理解するために少し努力します。コメントで指摘したように、解決はX**3 == 1収束しますが、デフォルトのイプシロンが小さすぎて、永遠にループするのを止めることができません。1.12e-16デフォルトのイプシロンの方が良いようです。

娯楽のために、不適切な Poly([-1, 3, -3, 1]) を試してください ... 3 つのルートはすべて (1+0j) に等しくなります ... 600 回以上の反復が必要で、最後の 10 回のエラーまたはそのため、左側のフィールドのはるか遠くから合理的なソリューションに到達するまで、反復は驚くほど飛び跳ねます。

于 2010-10-30T23:07:51.813 に答える