[我想要做的]是找到四次函数 ax^4 + bx^3 + cx^2 + dx + e 的唯一最小正实根。
[现有方法]
我的等式用于碰撞预测，最高次数为四次函数，如 f(x) = ax^4 + bx^3 + cx^2 + dx + e，而 a,b,c,d,e 系数可以是正数、负数或零（实数浮点值）。因此，我的函数 f(x) 可以是四次、三次或二次，具体取决于输入系数 a、b、c、d、e。
目前，我使用 NumPy 查找根，如下所示。
NumPy模块中的“root_output”可以是根据输入系数可能的所有实/复根。因此，我必须逐个查看“root_output”，并检查哪个根是最小的正实数值（根>0？）。 [问题] 我的程序需要多次执行numpy.roots([a,b,c,d,e])，但多次执行numpy.roots对于我的项目来说太慢了。当执行numpy.roots时，(a,b,c ,d ,e)的值总是在每次更改。
我尝试在Raspberry Pi2上运行代码。以下是处理时间的示例：
谢谢。 [解决方案]
目前，我使用 NumPy 查找根，如下所示。
import numpy
root_output = numpy.roots([a, b, c ,d ,e])
NumPy模块中的“root_output”可以是根据输入系数可能的所有实/复根。因此，我必须逐个查看“root_output”，并检查哪个根是最小的正实数值（根>0？）。 [问题] 我的程序需要多次执行numpy.roots([a,b,c,d,e])，但多次执行numpy.roots对于我的项目来说太慢了。当执行numpy.roots时，(a,b,c ,d ,e)的值总是在每次更改。
我尝试在Raspberry Pi2上运行代码。以下是处理时间的示例：
- 在PC上运行多次numpy.roots： 1.2秒
- 在Raspberry Pi2上运行多次numpy.roots： 17秒
谢谢。 [解决方案]
- 二次函数只需要实际正根（请注意除零）
def SolvQuadratic(a, b ,c):
d = (b**2) - (4*a*c)
if d < 0:
return []
if d > 0:
square_root_d = math.sqrt(d)
t1 = (-b + square_root_d) / (2 * a)
t2 = (-b - square_root_d) / (2 * a)
if t1 > 0:
if t2 > 0:
if t1 < t2:
return [t1, t2]
return [t2, t1]
return [t1]
elif t2 > 0:
return [t2]
else:
return []
else:
t = -b / (2*a)
if t > 0:
return [t]
return []
- 四次函数 对于四次函数，您可以使用纯Python / Numba版本，如@B.M.的下面的答案所示。 我还添加了另一个来自@B.M代码的Cython版本。 您可以将以下代码用作.pyx文件，然后编译它以获得比纯Python快约2倍的速度（请注意舍入问题）。
import cmath
cdef extern from "complex.h":
double complex cexp(double complex)
cdef double complex J=cexp(2j*cmath.pi/3)
cdef double complex Jc=1/J
cdef Cardano(double a, double b, double c, double d):
cdef double z0
cdef double a2, b2
cdef double p ,q, D
cdef double complex r
cdef double complex u, v, w
cdef double w0, w1, w2
cdef double complex r1, r2, r3
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
r=cmath.sqrt(-D/27+0j)
u=((-q-r)/2)**0.33333333333333333333333
v=((-q+r)/2)**0.33333333333333333333333
w=u*v
w0=abs(w+p/3)
w1=abs(w*J+p/3)
w2=abs(w*Jc+p/3)
if w0<w1:
if w2<w0 : v = v*Jc
elif w2<w1 : v = v*Jc
else: v = v*J
r1 = u+v-z0
r2 = u*J+v*Jc-z0
r3 = u*Jc+v*J-z0
return r1, r2, r3
cdef Roots_2(double a, double complex b, double complex c):
cdef double complex bp
cdef double complex delta
cdef double complex r1, r2
bp=b/2
delta=bp*bp-a*c
r1=(-bp-delta**.5)/a
r2=-r1-b/a
return r1, r2
def SolveQuartic(double a, double b, double c, double d, double e):
"Ferrarai's Method"
"resolution of P=ax^4+bx^3+cx^2+dx+e=0, coeffs reals"
"First shift : x= z-b/4/a => P=z^4+pz^2+qz+r"
cdef double z0
cdef double a2, b2, c2, d2
cdef double p, q, r
cdef double A, B, C, D
cdef double complex y0, y1, y2
cdef double complex a0, b0
cdef double complex r0, r1, r2, r3
z0=b/4.0/a
a2,b2,c2,d2 = a*a,b*b,c*c,d*d
p = -3.0*b2/(8*a2)+c/a
q = b*b2/8.0/a/a2 - 1.0/2*b*c/a2 + d/a
r = -3.0/256*b2*b2/a2/a2 + c*b2/a2/a/16 - b*d/a2/4+e/a
"Second find y so P2=Ay^3+By^2+Cy+D=0"
A=8.0
B=-4*p
C=-8*r
D=4*r*p-q*q
y0,y1,y2=Cardano(A,B,C,D)
if abs(y1.imag)<abs(y0.imag): y0=y1
if abs(y2.imag)<abs(y0.imag): y0=y2
a0=(-p+2*y0)**.5
if a0==0 : b0=y0**2-r
else : b0=-q/2/a0
r0,r1=Roots_2(1,a0,y0+b0)
r2,r3=Roots_2(1,-a0,y0-b0)
return (r0-z0,r1-z0,r2-z0,r3-z0)
[Ferrari方法的问题] 当四次方程的系数为[0.00614656，-0.0933333333333，0.527664995846，-1.31617928376，1.21906444869]时，使用numpy.roots和Ferrari方法得到的输出完全不同（numpy.roots是正确的输出）。
import numpy as np
import cmath
J=cmath.exp(2j*cmath.pi/3)
Jc=1/J
def ferrari(a,b,c,d,e):
"Ferrarai's Method"
"resolution of P=ax^4+bx^3+cx^2+dx+e=0, coeffs reals"
"First shift : x= z-b/4/a => P=z^4+pz^2+qz+r"
z0=b/4/a
a2,b2,c2,d2 = a*a,b*b,c*c,d*d
p = -3*b2/(8*a2)+c/a
q = b*b2/8/a/a2 - 1/2*b*c/a2 + d/a
r = -3/256*b2*b2/a2/a2 +c*b2/a2/a/16-b*d/a2/4+e/a
"Second find y so P2=Ay^3+By^2+Cy+D=0"
A=8
B=-4*p
C=-8*r
D=4*r*p-q*q
y0,y1,y2=Cardano(A,B,C,D)
if abs(y1.imag)<abs(y0.imag): y0=y1
if abs(y2.imag)<abs(y0.imag): y0=y2
a0=(-p+2*y0)**.5
if a0==0 : b0=y0**2-r
else : b0=-q/2/a0
r0,r1=Roots_2(1,a0,y0+b0)
r2,r3=Roots_2(1,-a0,y0-b0)
return (r0-z0,r1-z0,r2-z0,r3-z0)
#~ @jit(nopython=True)
def Cardano(a,b,c,d):
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
r=cmath.sqrt(-D/27+0j)
u=((-q-r)/2)**0.33333333333333333333333
v=((-q+r)/2)**0.33333333333333333333333
w=u*v
w0=abs(w+p/3)
w1=abs(w*J+p/3)
w2=abs(w*Jc+p/3)
if w0<w1:
if w2<w0 : v*=Jc
elif w2<w1 : v*=Jc
else: v*=J
return u+v-z0, u*J+v*Jc-z0, u*Jc+v*J-z0
#~ @jit(nopython=True)
def Roots_2(a,b,c):
bp=b/2
delta=bp*bp-a*c
r1=(-bp-delta**.5)/a
r2=-r1-b/a
return r1,r2
coef = [0.00614656, -0.0933333333333, 0.527664995846, -1.31617928376, 1.21906444869]
print("Coefficient A, B, C, D, E", coef)
print("")
print("numpy roots: ", np.roots(coef))
print("")
print("ferrari python ", ferrari(*coef))