La formule de Cardano ne fonctionne pas avec numpy?

Question

--- en utilisant python 3 ---

Après les équations here, j'ai essayé de trouver toutes les racines réelles d'un troisième ordre polynomiale arbitraire. Malheureusement, ma mise en œuvre ne donne pas le résultat correct et je ne trouve pas l'erreur. Peut-être que vous êtes en mesure de le repérer en un clin d'œil et dites-moi.

(. Comme vous le remarquez, que les racines de la courbe verte sont mal)

Avec mes meilleures salutations

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 

L'exemple où vous voyez cela ne fonctionne pas:

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() 

Sortie:

[[ 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]] 

Answer 1

Voici mon coup à la solution. Votre code échoue pour le cas où R + np.sqrt(D) ou R - np.sqrt(D) est négatif. La raison en est this post. Fondamentalement, si vous faites a**(1/3) où a est négatif, numpy renvoie un nombre complexe. Cependant, nous voulons que S et T soient réels puisque la racine cubique d'un nombre réel négatif est simplement un nombre réel négatif (ignorons théorème de De Moivre pour le moment et concentrons-nous sur le code et non sur les maths). La façon de contourner le problème est de vérifier si S est réel, de le convertir en réel et de passer le S à la fonction from scipy.special import cbrt. De même pour T. code Exemple:

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() 

AVERTISSEMENT: La racine de sortie donne un avertissement, que vous pouvez probablement ignorer. La sortie est correcte. Cependant, le tracé montre une racine supplémentaire pour certaines raisons. Cela est probablement dû à votre code de traçage. Les racines imprimées semblent bien cependant.