在python中找到四次多项式的最小正实根的最快方法 [英] fastest way to find the smallest positive real root of quartic polynomial 4 degree in python
问题描述
[我想要的] 是找到四次函数 a x ^ 4 + b x ^的唯一最小的正实根3 + c x ^ 2 + d x + e
[What I want] is to find the only one smallest positive real root of quartic function ax^4 + bx^3 + cx^2 + dx + e
[现有方法] 我的方程式是用于碰撞预测的，最大程度是四次函数为 f(x) = a x ^ 4 + b x ^ 3 + c x ^ 2 + d x + e 和 a，b，c，d，e 系数可以为正/负/零(实际浮点值).因此，根据输入系数a，b，c，d和e，我的函数 f(x)可以是四次，三次或二次.
[Existing Method] My equation is for collision prediction, the maximum degree is quartic function as f(x) = ax^4 + bx^3 + cx^2 + dx + e and a,b,c,d,e coef can be positive/negative/zero (real float value). So my function f(x) can be quartic, cubic, or quadratic depending on a, b, c ,d ,e input coefficient.
当前，我使用numpy来查找根，如下所示.
Currently I use numpy to find root as below.
导入numpy
root_output = numpy.roots([a，b，c，d，e])
root_output = numpy.roots([a, b, c ,d ,e])
numpy模块中的" root_output "可以是所有可能的实数/复数根，具体取决于输入系数.因此，我必须一一查看" root_output "，并检查哪个根是最小的实际正值(root> 0?)
The "root_output" from numpy module can be all possible real/complex roots depending on input coefficient. So I have to look at "root_output" one by one, and check which root is the smallest real positive value (root>0?)
[问题] 我的程序需要多次执行 numpy.roots([a，b，c，d，e])，因此对于我的项目而言，执行numpy.roots的许多次都太慢了.和(a，b，c，d，e)的值每次执行numpy.roots时都会更改
[The Problem] My program need to execute numpy.roots([a, b, c, d, e]) many times, so many times of executing numpy.roots is too slow for my project. and (a, b, c ,d ,e) value is always changed every time when executing numpy.roots
我的尝试是在Raspberry Pi2上运行代码.下面是处理时间的示例.
My attempt is to run the code on Raspberry Pi2. Below is example of processing time.
- 在PC上多次运行numpy.roots: 1.2秒
- 在Raspberry Pi2上多次运行numpy.roots: 17秒
- Running many many times of numpy.roots on PC: 1.2 seconds
- Running many many times of numpy.roots on Raspberry Pi2: 17 seconds
能否请您指导我如何在最快的解决方案中找到最小的正实根?使用scipy.optimize或实施某种算法来加快查找根源或您的任何建议的方法都会很棒.
Could you please guide me how to find the smallest positive real root in the fastest solution? Using scipy.optimize or implement some algorithm to speed up finding root or any advice from you will be great.
谢谢.
[解决方案]
-
二次函数仅需要实数正根(请注意除以零)
Quadratic function only need real positive roots (please aware of division by zero)
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倍的速度(请注意四舍五入的问题).
Quartic Function for quartic function, you can use pure python/numba version as the below answer from @B.M.. I also add another cython version from @B.M's code. You can use below code as .pyx file and then compile it to get about 2x faster than pure python (please aware of rounding issues).
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是正确的输出).
[Problem of Ferrari's method] We're facing the problem when the coefficients of quartic equation is [0.00614656, -0.0933333333333, 0.527664995846, -1.31617928376, 1.21906444869] the output from numpy.roots and ferrari methods is entirely different (numpy.roots is correct output).
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))
推荐答案
另一个答案:
使用分析方法(法拉利，
do it with analytic methods (Ferrari,Cardan), and speed the code with Just in Time compilation (Numba) :
让我们首先看看改进之处:
Let see the improvement first :
In [2]: P=poly1d([1,2,3,4],True)
In [3]: roots(P)
Out[3]: array([ 4., 3., 2., 1.])
In [4]: %timeit roots(P)
1000 loops, best of 3: 465 µs per loop
In [5]: ferrari(*P.coeffs)
Out[5]: ((1+0j), (2-0j), (3+0j), (4-0j))
In [5]: %timeit ferrari(*P.coeffs) #pure python without jit
10000 loops, best of 3: 116 µs per loop
In [6]: %timeit ferrari(*P.coeffs) # with numba.jit
100000 loops, best of 3: 13 µs per loop
然后是丑陋的代码:
对于订单4:
@jit(nopython=True)
def ferrari(a,b,c,d,e):
"resolution of P=ax^4+bx^3+cx^2+dx+e=0"
"CN all coeffs real."
"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 X so P2=AX^3+BX^2+C^X+D=0"
A=8
B=-4*p
C=-8*r
D=4*r*p-q*q
y0,y1,y2=cardan(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.real)**.5
if a0==0 : b0=y0**2-r
else : b0=-q/2/a0
r0,r1=roots2(1,a0,y0+b0)
r2,r3=roots2(1,-a0,y0-b0)
return (r0-z0,r1-z0,r2-z0,r3-z0)
对于订单3:
J=exp(2j*pi/3)
Jc=1/J
@jit(nopython=True)
def cardan(a,b,c,d):
u=empty(2,complex128)
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=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
对于订单2:
@jit(nopython=True)
def roots2(a,b,c):
bp=b/2
delta=bp*bp-a*c
u1=(-bp-delta**.5)/a
u2=-u1-b/a
return u1,u2
可能需要进一步测试，但是有效.
Probably needs to be test furthermore, but efficient.
这篇关于在python中找到四次多项式的最小正实根的最快方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
