Решение кубического уравнения ⇐ Python

[Anonymous]

В рамках программы, которую я пишу, мне нужно точно решить кубическое уравнение (вместо использования числового средства поиска корня):

Код: Выделить всё

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 решит такое уравнение, но мне нужно сделать это для миллионов уравнений и так Мне нужно реализовать это для работы с массивами коэффициентов.

Подробнее здесь: https://stackoverflow.com/questions/182 ... c-equation