未能实施Cardano方法。复数

的立方根为了提高三次方程np.roots的表现，我尝试实施Cardan(o) Method：未能实施Cardano方法。复数

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

它运作良好时，根实：

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))

但在失败复杂的情况：

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的评价d v由立方根。 Theory说uv=-p/3，但这里uv=pJ/3：(u,v)是不是一个好根。

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

编辑

@Sally文章后，我可以准确的问题。好的一对并不总是(u,v)，它可以是(u,vJ)或(uJ,v)。所以问题可能是：

是否有一个简单的方法来捕捉好对？

而终极：目前，通过编写此代码与Numba，它比np.roots快20倍。

有没有更好的算法来计算三次方程的三个根？

来源

2016-03-05 B. M.

你在寻找类似De Moivre定理的东西吗？ – AMACB

不是真的。 'sin'和'cos'很耗时。我看着一个快速的方法来知道'（v，Jv，J？v）'必须与'u'相关联以找到好的根。 –

他们不应该是，除非你使用他们*很多*：http://stackoverflow.com/a/3836083/5277935 – AMACB

您正确识别该问题：有立方根的三个可能的值在复平面中，产生9个可能的对((-q+r)/2)**(1/3)和((-q-r)/2)**(1/3)。在这9个中，只有3对导致了正确的根：即u * v = -p/3。一个简单的解决方法是用v=-p/(3*u)替换v的公式。这可能也是一个加速：分裂应该比立方根快。

然而u可能等于或接近于零，在这种情况下，划分变得可疑。事实上，在你的第一个例子中，它使精度稍微更差。下面是一个数字稳健的方法：在return语句之前插入这两行。

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绝对值。也许可以通过存储三个候选人来略微提高性能，这样胜者就不必重新计算。

来源

2016-03-06 04:52:04

感谢您的确认和建议。我想有一个更直接和对称的方法，但我目前没有看到它。我在这个目标中编辑我的帖子。 –

最终版本在这里。 http://stackoverflow.com/questions/35795663/fastest-way-to-find-the-smallest-positive-real-root-of-quartic-polynomial-4-degr/35829660#35829660。我更喜欢避免用“u”来限制舍入问题，但精度不是很好。 –

由于作为平方根之一的r的符号是空闲的（分别为：的u和v的角色是可以互换的u+v和u^3,v^3为二次多项式的根）设置

u3 = (abs(q+r)>abs(q-r))? -(q+r) : -(q-r) 

u = u3**(1/3) 
v = -p/(3*u)

这可以确保除数总是尽可能大，减少了错误的商和最小化零（零）附近可能成为问题的案例数量。

来源

2016-03-06 12:58:38 LutzL

谢谢。然而它不会在'q = p = 0'上失败，这并不罕见。 –

当然，你必须事先检查，抓住并处理这些特殊情况。人们已经尝试过卡尔达诺的健壮实现，因为我在某处读过，它变得相当长。（= q）**（1/3）* J^k'，'q == 0 => z = 0，+ - （ - p）**（1/2）'等 – LutzL

未能实施Cardano方法。复数

回答

相关问题