La fórmula cúbica da un resultado erróneo (triplemente comprobado)

Question

Me gustaría resolver $ax^3 + bx^2 + cx + d = 0$ utilizando la fórmula cúbica.

He codificado tres versiones de esta fórmula, descritas en tres fuentes: MathWorld , EqWorld ,
y en el libro "El intento inalcanzable de evitar el Casus Irreducibilis para las ecuaciones cúbicas".

Aunque obtengo resultados idénticos en todas las versiones, estos resultados son incorrectos.

Por ejemplo, para $a=1$ , $b=2$ , $c=3$ , $d=4$ ,

Encuentro raíces incorrectas:
$x_1 = -0.1747 - 0.8521i$ ,
$x_2 = 0.4270 + 1.1995i$ ,
$x_3 = -2.2523 - 0.3474i$ .

Las raíces correctas son:
$x_1 = -1.6506$ ,
$x_2 = -0.1747 + 1.5469i$ ,
$x_3 = -0.1747 - 1.5469i$

En caso de que esté interesado, el código real está abajo. Gracias por su ayuda.

%% Wolfram version

Q = (3*c - b^2) / 9;
R = (9*b*c - 27*d - 2*b^3) / 54;

D = Q^3 + R^2;
S = (R + sqrt(D))^(1/3);
T = (R - sqrt(D))^(1/3);

x1 = - b/3 + (S + T);
x2 = - b/3 - (S + T) / 2 + sqrt(-3) * (S - T) / 2;
x3 = - b/3 - (S + T) / 2 - sqrt(-3) * (S - T) / 2;

%% Book version

omega1 = - 1/2 + sqrt(-3)/2;
omega2 = - 1/2 - sqrt(-3)/2;

p = (3*a*c - b^2) / (3*a^2);
q = (2*b^3 - 9*a*b*c + 27*(a^3)*d) / (27*a^3);

r = sqrt(q^2/4 + p^3/27);
s = (-q/2 + r)^(1/3);
t = (-q/2 - r)^(1/3);

x1 =        s +        t - b/(3*a);
x2 = omega1*s + omega2*t - b/(3*a);
x3 = omega2*s + omega1*t - b/(3*a);

%% Eqworld version

p = - 1/3  * (b/a)^2 + (c/a);
q =   2/27 * (b/a)^3 - (b*c)/(3*a^2) + d/a;

D = (p/3)^3 + (q/2)^2;
A = (-q/2 + sqrt(D))^(1/3);
B = (-q/2 - sqrt(D))^(1/3);

y1 = A + B;
y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

x1 = y1 - b / (3*a);
x2 = y2 - b / (3*a);
x3 = y3 - b / (3*a);

Answer 1

El problema ocurre cuando se toma una raíz cúbica: en la versión de Wolfram, cuando se calcula $$T = \sqrt[3]{R - \sqrt{D}} = \sqrt[3]{-\frac{35}{27}-\sqrt{\frac{50}{27}}}$$ y en los lugares correspondientes de las otras versiones de la fórmula cúbica.

La expresión dentro de la raíz cúbica es negativa (aproximadamente $-2.66$ ), y la fórmula cúbica funciona si se toma la raíz cúbica real (aproximadamente $-1.39$ ). Pero escribir $(R - \sqrt D)^{1/3}$ en muchos sistemas de álgebra computacional, en cambio, calcula el principal raíz cúbica: la raíz compleja con mayor parte real. Esta será la raíz cúbica real para un número positivo, pero aquí nos da aproximadamente $0.69 - 1.2 i$ y ese es el origen de su error.

No puedo identificar el lenguaje que estás usando al verlo, así que no sé cuál es la mejor manera de evitar este problema. En Mathematica, hay una función CubeRoot que siempre tomará la raíz real de un número real.

Siempre se puede hacer con algunos condicionales: definir $T$ para ser

• $(R - \sqrt{D})^{1/3}$ , si $R - \sqrt D \ge 0$ y
• $-(\sqrt D - R)^{1/3}$ , si $R - \sqrt D < 0$ .

(Y haz algo similar en todos los demás lugares en los que tomes una raíz cúbica).

---

Lo anterior funcionará para coeficientes reales; para coeficientes complejos (o incluso cuando $D$ es negativo), queremos ser más cuidadosos, porque entonces podría no existir una raíz real.

La idea es que aquí nos encontramos con problemas, no porque no estuviéramos tomando la raíz real necesariamente, sino porque tomamos dos raíces cúbicas que "no coincidían entre sí". Podemos reescribir la fórmula cúbica de varias maneras para que sólo tomemos un raíz cúbica, y expresar la otra raíz cúbica que tendríamos que tomar en términos de la primera. Entonces el problema no se produce.

Esto nos da una mejor manera de resolver el problema en código: después de computar $S = (R + \sqrt D)^{1/3}$ , set $T = -Q/S$ . La justificación es la siguiente: ya que $S^3 = R + \sqrt D$ y $T^3 = R - \sqrt D$ tenemos $S^3 T^3 = R^2 - D = -Q^3$ . Ingenuamente, podemos cancelar las raíces cúbicas y suponer $ST = -Q$ y esto funciona.

Ver también esta pregunta y su respuesta.

Answer 2

Esto no se ve bien.

y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

¿Quieres decir cbrt(-3)?