Algoritmo numéricamente estable para resolver la ecuación cuadrática cuando $a$ es muy pequeño o $0$

Question

Resolver $a x^2 + bx +c=0$ para $x$ da

$$x = \frac{-b \pm \sqrt{b^2-4ac}}{2a} \text{, para } a \ne 0$$

Pero para $a = 0$ obtenemos

$$x=-\frac{c}{b}$$

¿Cómo implementar un algoritmo numéricamente estable para calcular $x$ a partir de $a,b,c$ que permita que $a$ esté cerca de cero o sea cero?

Answer 1

Las respuestas anteriores identifican correctamente que existen dos fórmulas cuadráticas, y que cada una tiene un rango diferente de estabilidad numérica; sin embargo, se pierden las otras sutiles inestabilidades.

Para recapitular: la fórmula cuadrática estándar, \begin{align*} x=\frac{-b\pm\sqrt{b^2-4\,a\,c}}{2\,a}, \end{align*} no es única. Se puede generar una segunda fórmula multiplicando el numerador y el denominador por el "conjugado" del numerador, es decir, \begin{align*} x&=\frac{-b\pm\sqrt{b^2-4\,a\,c}}{2\,a}\times\frac{-b\mp\sqrt{b^2-4\,a\,c}}{-b\mp\sqrt{b^2-4\,a\,c}}= \frac{2\,c}{-b\mp\sqrt{b^2-4\,a\,c}}, \end{align*} donde $\mp=-(\pm)$. Estas dos expresiones son analíticamente equivalentes; sin embargo, cada una es numéricamente inestable para ciertos valores de los coeficientes.

La inestabilidad discutida en las otras publicaciones es si el término, $4\,a\,c$, es pequeño en comparación con $b^2$ y el signo del término radical y $b$ son los mismos, entonces ocurre cancelación catastrófica.

Afortunadamente, las dos fórmulas cuadráticas tienen signos opuestos en el término radical para las mismas raíces, por lo que es posible evitar la cancelación catastrófica seleccionando la forma estable. Es decir,

\begin{align*} x_1 &= \frac{-b-{\rm sign}(b)\sqrt{b^2-4\,a\,c}}{2a}& x_2 = \frac{c}{a\,x_1} \end{align*}

Para mayor claridad, existen al menos tres inestabilidades numéricas más en esta fórmula:

1. Cuando $a=0: así la ecuación es lineal y tiene a lo sumo una raíz dada por $-c/b$. Se debe tener en cuenta que si $b=0$ también, la ecuación es en realidad $c=0$ y deja $x$ indefinido/sin restricciones.
2. Cuando $a\neq0$ y $c=0: así una raíz está en $x=0$. En este caso, la segunda forma de la ecuación cuadrática dará un NaN ($0/0$) para la segunda raíz. Por lo tanto, es mejor determinar la segunda raíz factorizando la raíz en cero para dar $-b/a$.
3. Desbordamiento: la representación de punto flotante puede desbordar/sub-desbordar durante el cálculo y esto es más crítico para el coeficiente $b$. Definiendo $\mathcal{R}_{max}$ como el número de punto flotante representable más grande ($\approx1.8\times10^{308}$ para doble precisión) entonces $b^2$ desbordará a infinito en $|b|>\sqrt{\mathcal{R}_{max}}$ de acuerdo con la especificación IEEE 754 (que es "solo" $\approx1.304\times10^{154}$). Puede que no definas esto como una inestabilidad, más bien como un error por estar fuera de rango, ¡pero creo que el rango superior reducido de $b$ es sorprendente y realmente me tomó por sorpresa!
4. ??: Estoy bastante seguro de que me he perdido algunas más, la aritmética de punto flotante puede ser muy sutil.

He escrito una biblioteca de álgebra en tiempo de compilación y puedes echarle un vistazo a mi implementación de la fórmula cuadrática. A continuación, se muestra una implementación "más estable" de la fórmula cúbica, descaradamente robada de D. Herbison-Evans. Para polinomios de orden superior me baso en técnicas de bisección basadas en el teorema de Budan, ¡que no es tan difícil de hacer estable!

Answer 2

A veces la Fórmula de Citardauq $$x=\frac{2c}{-b\mp \sqrt{b^2-4ac}}$$ puede ser útil. Se elige el signo que no da lugar a cancelaciones en el denominador.

Answer 3

Tenga en cuenta que solo para uno de $x_1=\frac{-b-\sqrt{b^2-4ac}}{2a}$ y $x_2=\frac{-b+\sqrt{b^2-4ac}}{2a}$ se realiza la "resta cerrada", por lo que puede calcular uno de ellos y el segundo usando la fórmula de Viète: $x_1 x_2=\frac{c}{a}$.