Формула Кардано не работает с numpy?

--- используя питон 3 ---

Следуя уравнениям здесь, я попытался найти все действительные корни произвольного многочлена третьего порядка . К сожалению, моя реализация не дает правильного результата, и я не могу найти ошибку. Может быть, вы сможете заметить это в мгновение ока и сказать мне.

(Как вы заметили, неверны только корни зеленой кривой.)

С наилучшими пожеланиями

import numpy as np
def find_cubic_roots(a,b,c,d):
    # with ax³ + bx² + cx + d = 0
    a,b,c,d = a+0j, b+0j, c+0j, d+0j
    all_ = (a != np.pi)

    Q = (3*a*c - b**2)/ (9*a**2)
    R = (9*a*b*c - 27*a**2*d - 2*b**3) / (54 * a**3)
    D = Q**3 + R**2
    S = (R + np.sqrt(D))**(1/3)
    T = (R - np.sqrt(D))**(1/3)

    result = np.zeros(tuple(list(a.shape) + [3])) + 0j
    result[all_,0] = - b / (3*a) + (S+T)
    result[all_,1] = - b / (3*a)  - (S+T) / 2 + 0.5j * np.sqrt(3) * (S - T)
    result[all_,2] = - b / (3*a)  - (S+T) / 2 -  0.5j * np.sqrt(3) * (S - T)

    return result

Пример, где вы видите, что это не работает:

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
a = np.array([2.5])
b = np.array([-5])
c = np.array([0])

x = np.linspace(-2,3,100)
for i, d in enumerate([-8,0,8]):
    d = np.array(d)
    roots = find_cubic_roots(a,b,c,d)
    ax.plot(x, a*x**3 + b*x**2 + c*x + d, label = "a = %.3f, b = %.3f, c = %.3f, d = %.3f"%(a,b,c,d), color = colors[i])
    print(roots)
    ax.plot(x, x*0)
    ax.scatter(roots,roots*0,  s = 80)
ax.legend(loc = 0)
ax.set_xlim(-2,3)
plt.show()

Простой пример

Выход:

[[ 2.50852567+0.j        -0.25426283+1.1004545j -0.25426283-1.1004545j]]
[[ 2.+0.j  0.+0.j  0.-0.j]]
[[ 1.51400399+1.46763129j  1.02750817-1.1867528j  -0.54151216-0.28087849j]]

python numpy calculus
person varantir    schedule 13.09.2016    source источник
comment
это питон 3 или питон 2? потому что в питоне 2 1/3 = 0.   -  person Jean-François Fabre    schedule 13.09.2016
comment
python3, я напишу это выше.   -  person varantir    schedule 13.09.2016

Ответы (1)


arrow_upward
1
arrow_downward

Вот мой удар по решению. Ваш код не работает в случае, когда R + np.sqrt(D) или R - np.sqrt(D) имеет отрицательное значение. Причина в этом сообщении. В основном, если вы делаете a**(1/3), где a отрицательно, numpy возвращает комплексное число. Однако на самом деле мы хотим, чтобы S и T были реальными, поскольку кубический корень из отрицательного действительного числа — это просто отрицательное действительное число (давайте пока проигнорируем теорему Де Муавра и сосредоточимся на коде, а не на математике). ). Чтобы обойти это, нужно проверить, является ли S реальным, привести его к реальному и передать S функции from scipy.special import cbrt. Аналогично для T. Пример кода:

import numpy as np
import pdb
import math
from scipy.special import cbrt
def find_cubic_roots(a,b,c,d, bp = False):

    a,b,c,d = a+0j, b+0j, c+0j, d+0j
    all_ = (a != np.pi)

    Q = (3*a*c - b**2)/ (9*a**2)
    R = (9*a*b*c - 27*a**2*d - 2*b**3) / (54 * a**3)
    D = Q**3 + R**2
    S = 0 #NEW CALCULATION FOR S STARTS HERE
    if np.isreal(R + np.sqrt(D)):
        S = cbrt(np.real(R + np.sqrt(D)))
    else:
        S = (R + np.sqrt(D))**(1/3)
    T = 0 #NEW CALCULATION FOR T STARTS HERE
    if np.isreal(R - np.sqrt(D)):
        T = cbrt(np.real(R - np.sqrt(D)))
    else:
        T = (R - np.sqrt(D))**(1/3)

    result = np.zeros(tuple(list(a.shape) + [3])) + 0j
    result[all_,0] = - b / (3*a) + (S+T)
    result[all_,1] = - b / (3*a)  - (S+T) / 2 + 0.5j * np.sqrt(3) * (S - T)
    result[all_,2] = - b / (3*a)  - (S+T) / 2 -  0.5j * np.sqrt(3) * (S - T)
    #if bp:
        #pdb.set_trace()
    return result


import matplotlib.pyplot as plt
fig, ax = plt.subplots()
a = np.array([2.5])
b = np.array([-5])
c = np.array([0])
x = np.linspace(-2,3,100)
for i, d in enumerate([-8,0,8]):
    d = np.array(d)
    if d == 8:
        roots = find_cubic_roots(a,b,c,d, True)
    else:
        roots = find_cubic_roots(a,b,c,d)

    ax.plot(x, a*x**3 + b*x**2 + c*x + d, label = "a = %.3f, b = %.3f, c = %.3f, d = %.3f"%(a,b,c,d))
    print(roots)
    ax.plot(x, x*0)
    ax.scatter(roots,roots*0,  s = 80)
ax.legend(loc = 0)
ax.set_xlim(-2,3)
plt.show()

ОТКАЗ ОТ ОТВЕТСТВЕННОСТИ: выходной корень выдает некоторое предупреждение, которое вы вероятно можете игнорировать. Выход правильный. Однако на графике почему-то виден лишний корень. Вероятно, это связано с вашим кодом построения графика. Хотя отпечатанные корни выглядят нормально.

person TuanDT    schedule 13.09.2016
comment
@varantir с удовольствием. - person TuanDT; 14.09.2016