Databases
StackOverflow2013
Postsdbo
PODurand-kerner implementation doesn't work
 

Note that there are some explanatory texts on larger screens.

plurals
  1. PODurand-kerner implementation doesn't work
    primarykey
    data
    text
    <p>What's wrong with this implementation of the Durand-Kerner algorithm (<a href="http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method" rel="nofollow">here</a>) ?</p> <pre><code>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]) &lt; epsilon for i, n in enumerate(new)): return new roots = new </code></pre> <p>When I try it, I have to stop it with <code>KeyboardInterrupt</code> because it doesn't stop!<br> <code>poly</code> is a Polynomial instance of the <code>pypol</code> library.</p> <p>Thank you in advance, rubik</p> <p><strong>EDIT</strong>: Using a numpy polynomial it takes 9 iterations:</p> <pre><code>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)] </code></pre> <p>Using a pypol polynomial it never finishes (it is probably a bug in pypol):</p> <pre><code>In [3]: roots.d2(poly1d([1, -3, 3, -5])) ^C--------------------------------------------------------------------------- KeyboardInterrupt </code></pre> <p>but I can't find the bug!! </p> <p><strong>EDIT2</strong>: Comparing the <code>__call__</code> method with Martin's Poly:</p> <pre><code>&gt;&gt;&gt; p = Poly(-5, 3, -3, 1) &gt;&gt;&gt; from pypol import poly1d &gt;&gt;&gt; p2 = poly1d([1, -3, 3, -5]) &gt;&gt;&gt; for i in xrange(-100000, 100000): assert p(i) == p2(i) &gt;&gt;&gt; &gt;&gt;&gt; for i in xrange(-10000, 10000): assert p(complex(1, i)) == p2(complex(1, i)) &gt;&gt;&gt; for i in xrange(-10000, 10000): assert p(complex(i, i)) == p2(complex(i, i)) &gt;&gt;&gt; </code></pre> <p><strong>EDIT3</strong>: pypol works fine if the roots aren't complex numbers:</p> <pre><code>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)] </code></pre> <p>So it doesn't works only when the roots are complex numbers!</p> <p><strong>EDIT4</strong>: I wrote a slightly different implementation for numpy polynomials and saw that after one iteration the roots (of the Wikipedia polynomial) are different:</p> <pre><code>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 </code></pre> <p><strong>EDIT5</strong>: Hey! If I change the line: <code>if all(n == roots[i] ... )</code> into <code>if all(str(n) == str(roots[i]) ... )</code> it finishes and returns the right roots!!!</p> <pre><code>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)] </code></pre> <p>But the question is: why does it work with a different complex numbers comparation??</p> <p><strong>UPDATE</strong><br> Now it works, and I've done some tests:</p> <pre><code>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)) </code></pre>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. USrubik
    1. USrubik
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    2. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. PODurand-kerner implementation doesn't work
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. PODurand-kerner implementation doesn't work
      1. USaaronasterling
      1. VTFavorite
    3. VO
      singulars
      1. PODurand-kerner implementation doesn't work
      1. USrubik
      1. VTFavorite
    1. COCommenting because it's only a guess (or a copy/paste issue), but should the `if all...` line be indented, and thus in the `for i,r...` loop? Currently, `i` will always be the last value in `enumerate(roots)`.
      singulars
      1. PODurand-kerner implementation doesn't work
      1. USBenn
    2. COIt is not a copy/past issue, it is out of the for loop because I have to check _all_ the elements, not only the current one. But I can try `if r == roots[i] or abs(r - roots[i]) < epsilon: return new` too.
      singulars
      1. PODurand-kerner implementation doesn't work
      1. USrubik
    3. COIf I change that it returns only one incorrect root...
      singulars
      1. PODurand-kerner implementation doesn't work
      1. USrubik
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload