无法实现Cardano方法.复数的立方根 [英] Fail to implement Cardano method. Cube root of a complex number

问题描述

为了提高三次方程的np.roots性能，我尝试实现 Cardan(o)方法:

In order to improve np.roots performance on cubic equation, I try to implement Cardan(o) Method :

def cardan(a,b,c,d):
    #"resolve P=ax^3+bx^2+cx+d=0" 
    #"x=z-b/a/3=z-z0 => P=z^3+pz+q"
    z0=b/3/a
    a2,b2 = a*a,b*b    
    p=-b2/3/a2 +c/a
    q=(b/27*(2*b2/a2-9*c/a)+d)/a
    D=-4*p*p*p-27*q*q+0j
    r=sqrt(-D/27)
    J=-0.5+0.86602540378443871j # exp(2i*pi/3)
    u=((-q+r)/2)**(1/3)
    v=((-q-r)/2)**(1/3)
    return u+v-z0,u*J+v/J-z0,u/J+v*J-z0

当根为真时，它会很好地工作:

It works well when roots are real:

In [2]: P=poly1d([1,2,3],True)
In [3]: roots(P)
Out[3]: array([ 3.,  2.,  1.])
In [4]: cardan(*P)
Out[4]: ((3+0j), (1+0j), (2+1.110e-16j))

但是在复杂情况下失败:

But fails in the complex case :

In [8]: P=poly1d([1,-1j,1j],True)
In [9]: P
Out[9]: poly1d([ 1., -1.,  1., -1.])
In [10]: roots(P)
Out[10]: array([  1.0000e+00+0.j,   7.771e-16+1.j,   7.771e-16-1.j])
In [11]: cardan(*P)
Out[11]: ((1.366+0.211j),(5.551e-17+0.577j),(-0.366-0.788j))

我想问题是立方根对u和v的求值. 理论说uv=-p/3，但在这里uv=pJ/3:(u,v)不是很好的一对.

I guess that the problem is the evaluation of u and v by cube roots . Theory say uv=-p/3, but here uv=pJ/3: (u,v) is not a good pair of roots.

在所有情况下获得正确配对的最佳方法是什么?

What is the best way to obtain a correct pair in all cases ?

编辑

在@Sally发布之后，我可以解决问题.好的对并不总是(u,v)，它可以是(u,vJ)或(uJ,v).所以问题可能是:

After @Sally post, I can precise the problem. The good pair is not always (u,v), it can be (u,vJ) or (uJ,v). So the question could be :

• 有没有一种简单的方法来捕捉好配对?

最终:通过使用Numba编译此代码，它的速度比np.roots快20倍.

And ultimate : for the moment, by compiling this code with Numba, it's 20x faster than np.roots .

• 是否有更好的算法来计算三次方程的三个根?

推荐答案

您正确地识别了问题:在复杂平面中存在三次可能的立方根值，从而产生了((-q+r)/2)**(1/3)和((-q-r)/2)**(1/3)的9对可能.在这9个中，只有3对引出正确的根:即u * v = -p/3的根.一个简单的解决方法是将v的公式替换为v=-p/(3*u).这可能也是一种提速:除法应该比求立方根更快.

You correctly identified the problem: there are three possible values of cubic root in a complex plane, resulting in 9 possible pairs of ((-q+r)/2)**(1/3) and ((-q-r)/2)**(1/3). Of these 9, only 3 pairs lead to the correct roots: namely, those for which u*v = -p/3. An easy fix is to replace the formula for v with v=-p/(3*u). This is probably also a speed-up: division should be faster than taking cubic root.

但是u可能等于或接近于零，在这种情况下，该除法变得可疑.确实，在您的第一个示例中，它已经使精度稍差了.这是一种数值上可靠的方法:在return语句之前插入这两行.

However u might be equal or close to zero, in which case the division becomes suspect. Indeed, already in your first example it makes the precision slightly worse. Here is a numerically robust approach: insert these two lines before the return statement.

choices = [abs(u*v*J**k+p/3) for k in range(3)]
v = v*J**choices.index(min(choices))

这将遍历v的三个候选，选择使u*v+p/3的绝对值最小的那个.也许可以通过存储三个候选者来稍微改善此处的性能，从而不必重新计算获胜者.

This loops over the three candidates for v, picking whichever minimizes the absolute value of u*v+p/3. Perhaps one can slightly improve the performance here by storing the three candidates so that the winner does not have to be recomputed.

这篇关于无法实现Cardano方法.复数的立方根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！