Utt*_*wut 5 python algorithm math numpy scipy
[我想要的]是找到四次函数的唯一一个最小的正实根a x ^ 4 + b x ^ 3 + c x ^ 2 + d x + e
[现有方法] 我的方程是用于碰撞预测中,最大程度是四次函数为F(X) = 一个的x ^ 4 + b的x ^ 3 + c ^ X ^ 2 + d X + ë 和A,B,C,d ,e coef可以是正/负/零(实浮点值).所以我的函数f(x)可以是四次,三次或二次,取决于a,b,c,d,e输入系数.
目前我使用numpy来查找root,如下所示.
导入numpy
root_output = numpy.roots([a,b,c,d,e])
来自numpy模块的" root_output "可以是所有可能的实/复根,具体取决于输入系数.所以我必须逐个查看" root_output ",并检查哪个根是最小的实数正值(root> 0？)
[问题] 我的程序需要多次执行numpy.roots([a,b,c,d,e]),因此执行numpy.roots很多次对我的项目来说太慢了.每次执行numpy.roots时,(a,b,c,d,e)值总是会改变
我的尝试是在Raspberry Pi2上运行代码.以下是处理时间的示例.
您能否指导我如何在最快的解决方案中找到最小的正实根？使用scipy.optimize或实现一些算法来加速寻找根或任何建议将是伟大的.
谢谢.
[解]
二次函数只需要真正的正根(请注意除零)
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 []
Quartic Function for quartic function,你可以使用纯python/numba版本作为@BM的下面答案.我还从@ BM的代码中添加了另一个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)
[问题法拉利的方法]我们面临的问题时四次方程式的系数为[0.00614656,-0.0933333333333,0.527664995846,-1.31617928376,1.21906444869]从numpy.roots和法拉利方法的输出是完全不同的(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))
另一个答案：
\n\n使用分析方法（Ferrari、Cardan）来完成此操作，并使用即时编译（Numba）来加速代码：
\n\n首先让我们看看改进：
\n\nIn [2]: P=poly1d([1,2,3,4],True)\n\nIn [3]: roots(P)\nOut[3]: array([ 4., 3., 2., 1.])\n\nIn [4]: %timeit roots(P)\n1000 loops, best of 3: 465 \xc2\xb5s per loop\n\nIn [5]: ferrari(*P.coeffs)\nOut[5]: ((1+0j), (2-0j), (3+0j), (4-0j))\n\nIn [5]: %timeit ferrari(*P.coeffs) #pure python without jit\n10000 loops, best of 3: 116 \xc2\xb5s per loop \nIn [6]: %timeit ferrari(*P.coeffs) # with numba.jit\n100000 loops, best of 3: 13 \xc2\xb5s per loop\n然后是丑陋的代码：
\n\n对于订单 4：
\n\n@jit(nopython=True)\ndef ferrari(a,b,c,d,e):\n "resolution of P=ax^4+bx^3+cx^2+dx+e=0"\n "CN all coeffs real."\n "First shift : x= z-b/4/a => P=z^4+pz^2+qz+r"\n z0=b/4/a\n a2,b2,c2,d2 = a*a,b*b,c*c,d*d \n p = -3*b2/(8*a2)+c/a\n q = b*b2/8/a/a2 - 1/2*b*c/a2 + d/a\n r = -3/256*b2*b2/a2/a2 +c*b2/a2/a/16-b*d/a2/4+e/a\n "Second find X so P2=AX^3+BX^2+C^X+D=0"\n A=8\n B=-4*p\n C=-8*r\n D=4*r*p-q*q\n y0,y1,y2=cardan(A,B,C,D)\n if abs(y1.imag)<abs(y0.imag): y0=y1 \n if abs(y2.imag)<abs(y0.imag): y0=y2 \n a0=(-p+2*y0.real)**.5\n if a0==0 : b0=y0**2-r\n else : b0=-q/2/a0\n r0,r1=roots2(1,a0,y0+b0)\n r2,r3=roots2(1,-a0,y0-b0)\n return (r0-z0,r1-z0,r2-z0,r3-z0) \n对于订单 3：
\n\nJ=exp(2j*pi/3)\nJc=1/J\n\n@jit(nopython=True) \ndef cardan(a,b,c,d):\n u=empty(2,complex128)\n z0=b/3/a\n a2,b2 = a*a,b*b \n p=-b2/3/a2 +c/a\n q=(b/27*(2*b2/a2-9*c/a)+d)/a\n D=-4*p*p*p-27*q*q\n r=sqrt(-D/27+0j) \n u=((-q-r)/2)**0.33333333333333333333333\n v=((-q+r)/2)**0.33333333333333333333333\n w=u*v\n w0=abs(w+p/3)\n w1=abs(w*J+p/3)\n w2=abs(w*Jc+p/3)\n if w0<w1: \n if w2<w0 : v*=Jc\n elif w2<w1 : v*=Jc\n else: v*=J \n return u+v-z0, u*J+v*Jc-z0,u*Jc+v*J-z0\n对于订单 2：
\n\n@jit(nopython=True)\ndef roots2(a,b,c):\n bp=b/2 \n delta=bp*bp-a*c\n u1=(-bp-delta**.5)/a\n u2=-u1-b/a\n return u1,u2 \n可能需要进一步测试，但有效。
\n| 归档时间： |
|
| 查看次数： |
2377 次 |
| 最近记录： |