作为我正在编写的程序的一部分,我需要精确地求解三次方程(而不是使用数值根查找器):
a*x**3 + b*x**2 + c*x + d = 0.
我正在尝试从这里使用等式。但是,请考虑以下代码(这是Python,但这是非常通用的代码):
a = 1.0 b = 0.0 c = 0.2 - 1.0 d = -0.7 * 0.2 q = (3*a*c - b**2) / (9 * a**2) r = (9*a*b*c - 27*a**2*d - 2*b**3) / (54*a**3) print "q = ",q print "r = ",r delta = q**3 + r**2 print "delta = ",delta # here delta is less than zero so we use the second set of equations from the article: rho = (-q**3)**0.5 # For x1 the imaginary part is unimportant since it cancels out s_real = rho**(1./3.) t_real = rho**(1./3.) print "s [real] = ",s_real print "t [real] = ",t_real x1 = s_real + t_real - b / (3. * a) print "x1 = ", x1 print "should be zero: ",a*x1**3+b*x1**2+c*x1+d
但是输出是:
q = -0.266666666667 r = 0.07 delta = -0.014062962963 s [real] = 0.516397779494 t [real] = 0.516397779494 x1 = 1.03279555899 should be zero: 0.135412149064
因此输出不为零,因此x1实际上不是解决方案。维基百科文章中有错误吗?
ps:我知道numpy.roots可以解决这类方程式,但是我需要对数百万个方程式执行此操作,因此我需要实现它才能处理系数数组。
维基百科的符号(rho^(1/3), theta/3)并不意味着这rho^(1/3)是真实的部分,theta/3而是虚构的部分。而是在极坐标中。因此,如果您想要真正的零件,则可以选择rho^(1/3) * cos(theta/3)。
(rho^(1/3), theta/3)
rho^(1/3)
theta/3
rho^(1/3) * cos(theta/3)
我对您的代码进行了这些更改,并且对我有用:
theta = arccos(r/rho) s_real = rho**(1./3.) * cos( theta/3) t_real = rho**(1./3.) * cos(-theta/3)
(当然,s_real = t_real这里cos是偶数。)
s_real = t_real
cos