Raíz rápida y robusta de un polinomio cúbico con restricciones.

Question

Estoy buscando una rápida y eficaz método para encontrar una raíz de un polinomio cúbico

$x^3 + px^2 + qx + r$

Para hacer la búsqueda más robusto y más rápido, me gustaría aprovechar estas propiedades:

1. El polinomio tiene exactamente un real positivo de la raíz (otras dos raíces son complejas o real y negativo).
2. Sólo el valor positivo de la raíz es necesario.
3. Hay una buena estimación inicial sobre el valor de la raíz.

Hasta ahora mi enfoque fue aplicar directamente el método de Newton a la función con el valor inicial y que me daría un resultado decente en sólo un par de iteraciones.

Sin embargo, en algunos casos, la iteración sería la causa de la conjetura de la raíz para saltar más cerca de uno de los negativos de las raíces y el método incorrectamente empezar a converger hacia esa raíz en su lugar. Si bien es posible detectar esta situación, es difícil llevar la iteración de vuelta en la pista de la derecha.

Hay un artículo interesante acerca de la solución de cuárticas y cúbicas y también un ejemplo de aplicación, pero los métodos son muy genérico, demasiado lento para mis necesidades y tener la robustez de los problemas de su propio.

Sería interesante saber si hay una manera de hacer mi iterativo de búsqueda más robusto o, posiblemente, de que hay un más rápido método de análisis de tomar ventaja de las propiedades adicionales.

Answer 1

Creo que se puede evitar el método de Newton-Raphson en total, desde cúbico es solucionable en los radicales.

Aquí está la completa algoritmo (que trabaja bajo las restricciones descritas en el problema), en Mathematica:

Clear[PositiveCubicRoot];
PositiveCubicRoot[p_, q_, r_] := 
 Module[{po3 = p/3, a, b, det, abs, arg},
  b = ( po3^3 - po3 q/2 + r/2);
  a = (-po3^2 + q/3);
  det = a^3 + b^2;
  If[det >= 0,
   det = Power[Sqrt[det] - b, 1/3];
   -po3 - a/det + det
   ,
   (* evaluate real part, imaginary parts cancel anyway *)
   abs = Sqrt[-a^3];
   arg = ArcCos[-b/abs];
   abs = Power[abs, 1/3];
   abs = (abs - a/abs);
   arg = -po3 + abs*Cos[arg/3]
   ]
  ]

Entonces

In[222]:= PositiveCubicRoot[-2.52111798, -71.424692, -129.51520]

Out[222]= 10.499

Sin embargo, si el método de Newton-Raphson el método debe ser utilizado, una buena estimación inicial es imprescindible. División binaria, es un buen método para aislar la raíz en el caso en cuestión. Empezamos con un punto arbitrario $x$, elegí $x=1$, el doble, mientras que el polinomio en ese punto es negativo. A continuación, hacer la división binaria de un número determinado de veces (el código siguiente hace dos veces). En última instancia, polaco fuera con el de Newton-Raphson iteraciones:

In[283]:= 
NewtonRaphsonStartingPoint[{p_, q_, r_}] := Module[{x1=0, x2=1,f1,f2,xm,fm},
   f1 = r + x1 (q + (p + x1) x1);
   While[(f2 = r + x2 (q + (p + x2) x2)) <= 0, 
      x1 = x2; f1 = f2; x2 = 2 x2];
   Do[xm = (x1 + x2)/2; fm = r + xm (q + (p + xm) xm); 
    If[fm <= 0, f1 = fm; x1 = xm, f2 = fm; x2 = xm], {i, 2}];
   (f2 x2 - f1 x1)/(f2 - f1)
];
NewtonRaphsonIterate[{p_, q_, r_}, x0_Real] := 
 FixedPoint[
  Function[x, x - (r + x (q + (p + x) x))/(q + (2 p + 3 x) x)], x0]

In[285]:= 
NewtonRaphson[p_, q_, r_] := 
 NewtonRaphsonIterate[{p, q, r}, NewtonRaphsonStartingPoint[{p, q, r}]]

In[286]:= NewtonRaphson[-2.52111798, -71.424692, -129.51520]

Out[286]= 10.499

Answer 2

Si su función$f$ es convexa y aumenta, o cóncava y disminuye, en un intervalo$[a,b]$ que contiene una solución de$f(x) = 0$, el método de Newton que comienza en$b$ convergerá a solución. De manera similar, si es cóncavo y en aumento, o convexo y en disminución, y comienzas en$a$. Para un cúbico, es fácil identificar los puntos críticos y la inflexión y, por lo tanto, encontrar dicho intervalo.