Encontrar raíces polinómicas es un problema de larga data que ha sido objeto de mucha investigación a lo largo de la historia. Prueba de ello es que hasta el siglo XIX, álgebra significaba esencialmente teoría de ecuaciones polinómicas .
Encontrar la raíz de un polinomio lineal (grado uno) es fácil y solo necesita una división: la ecuación general tiene solución Para polinomios cuadráticos (grado dos), la fórmula cuadrática produce una solución, pero su evaluación numérica puede requerir cierto cuidado para asegurar la estabilidad numérica . Para los grados tres y cuatro, hay soluciones en forma cerrada en términos de radicales , que generalmente no son convenientes para la evaluación numérica, por ser demasiado complicadas e involucrar el cálculo de varias raíces n -ésimas cuyo cálculo no es más fácil que el cálculo directo de las raíces del polinomio (por ejemplo, la expresión de las raíces reales de un polinomio cúbico puede involucrar raíces cúbicas no reales ). Para polinomios de grado cinco o superior, el teorema de Abel-Ruffini afirma que, en general, no hay expresión radical de las raíces.
Así pues, salvo en el caso de grados muy bajos, la búsqueda de raíces de polinomios consiste en hallar aproximaciones de las raíces. Por el teorema fundamental del álgebra , un polinomio de grado n tiene exactamente n raíces reales o complejas contando las multiplicidades .
De ello se deduce que el problema de encontrar raíces para polinomios puede dividirse en tres subproblemas diferentes;
Para encontrar una raíz, el método de Newton y otros métodos iterativos generales generalmente funcionan bien.
Para encontrar todas las raíces, posiblemente el método más confiable es el algoritmo Francis QR que calcula los valores propios de la matriz compañera correspondiente al polinomio, implementado como el método estándar [1] en MATLAB .
El método más antiguo para encontrar todas las raíces es empezar por encontrar una sola raíz. Cuando se ha encontrado una raíz r , se puede eliminar del polinomio dividiendo el binomio x – r . El polinomio resultante contiene las raíces restantes, que se pueden encontrar iterando sobre este proceso. Sin embargo, excepto para grados bajos, esto no funciona bien debido a la inestabilidad numérica : el polinomio de Wilkinson muestra que una modificación muy pequeña de un coeficiente puede cambiar drásticamente no solo el valor de las raíces, sino también su naturaleza (real o compleja). Además, incluso con una buena aproximación, cuando se evalúa un polinomio en una raíz aproximada, se puede obtener un resultado que está lejos de estar cerca de cero. Por ejemplo, si un polinomio de grado 20 (el grado del polinomio de Wilkinson) tiene una raíz cercana a 10, la derivada del polinomio en la raíz puede ser del orden de esto implica que un error de en el valor de la raíz puede producir un valor del polinomio en la raíz aproximada que sea del orden de
Para evitar estos problemas, se han elaborado métodos que calculan todas las raíces simultáneamente, con cualquier precisión deseada. Actualmente, el método más eficiente es el método de Aberth . Hay disponible una implementación gratuita con el nombre de MPSolve . Se trata de una implementación de referencia que puede encontrar rutinariamente las raíces de polinomios de grado mayor que 1000, con más de 1000 dígitos decimales significativos.
Los métodos para calcular todas las raíces se pueden utilizar para calcular raíces reales. Sin embargo, puede resultar difícil decidir si una raíz con una pequeña parte imaginaria es real o no. Además, como el número de raíces reales es, en promedio, proporcional al logaritmo del grado, [2] es un desperdicio de recursos informáticos calcular las raíces no reales cuando uno está interesado en las raíces reales.
El método más antiguo para calcular el número de raíces reales y el número de raíces en un intervalo resulta del teorema de Sturm , pero los métodos basados en la regla de los signos de Descartes y sus extensiones ( los teoremas de Budan y Vincent ) son generalmente más eficientes. Para encontrar raíces, todos proceden reduciendo el tamaño de los intervalos en los que se buscan raíces hasta obtener intervalos que contengan cero o una raíz. Luego, los intervalos que contienen una raíz pueden reducirse aún más para obtener una convergencia cuadrática del método de Newton a las raíces aisladas. Los principales sistemas de álgebra computacional ( Maple , Mathematica , SageMath , PARI/GP ) tienen cada uno una variante de este método como algoritmo predeterminado para las raíces reales de un polinomio.
La clase de métodos se basa en convertir el problema de encontrar raíces polinómicas al problema de encontrar valores propios de la matriz compañera del polinomio, [1] en principio, se puede utilizar cualquier algoritmo de valores propios para encontrar las raíces del polinomio. Sin embargo, por razones de eficiencia se prefieren métodos que emplean la estructura de la matriz, es decir, se pueden implementar en forma libre de matriz. Entre estos métodos se encuentran el método de potencia , cuya aplicación a la transposición de la matriz compañera es el método clásico de Bernoulli para encontrar la raíz de mayor módulo. El método de potencia inversa con desplazamientos, que encuentra primero alguna raíz más pequeña, es lo que impulsa la variante compleja ( cpoly ) del algoritmo de Jenkins-Traub y le da su estabilidad numérica. Además, tiene una convergencia rápida con el orden (donde es la proporción áurea ) incluso en presencia de raíces agrupadas. Esta rápida convergencia tiene un costo de tres evaluaciones polinomiales por paso, lo que resulta en un residuo de O (| f ( x )| 2+3 φ ) , que es una convergencia más lenta que con tres pasos del método de Newton.
El método más utilizado para calcular una raíz es el método de Newton , que consiste en las iteraciones del cálculo de
Partiendo de un valor bien elegido
Si f es un polinomio, el cálculo es más rápido cuando se utiliza el método de Horner o la evaluación con preprocesamiento para calcular el polinomio y su derivada en cada iteración.
Aunque la convergencia es generalmente cuadrática , puede converger mucho más lentamente o incluso no converger en absoluto. En particular, si el polinomio no tiene raíz real y es real, entonces el método de Newton no puede converger. Sin embargo, si el polinomio tiene una raíz real, que es mayor que la raíz real mayor de su derivada, entonces el método de Newton converge cuadráticamente a esta raíz más grande si es mayor que esta raíz más grande (hay formas fáciles de calcular un límite superior de las raíces, consulte Propiedades de las raíces de polinomios ). Este es el punto de partida del método de Horner para calcular las raíces.
Cuando se ha encontrado una raíz r , se puede utilizar la división euclidiana para eliminar el factor x – r del polinomio. Calcular una raíz del cociente resultante y repetir el proceso proporciona, en principio, una forma de calcular todas las raíces. Sin embargo, este esquema iterativo es numéricamente inestable; los errores de aproximación se acumulan durante las factorizaciones sucesivas, de modo que las últimas raíces se determinan con un polinomio que se desvía ampliamente de un factor del polinomio original. Para reducir este error, se puede, para cada raíz que se encuentre, reiniciar el método de Newton con el polinomio original y esta raíz aproximada como valor inicial.
Sin embargo, no hay garantía de que esto permita encontrar todas las raíces. De hecho, el problema de encontrar las raíces de un polinomio a partir de sus coeficientes puede ser muy mal condicionado . Esto se ilustra con el polinomio de Wilkinson : las raíces de este polinomio de grado 20 son los 20 primeros números enteros positivos; cambiar el último bit de la representación de 32 bits de uno de sus coeficientes (igual a -210) produce un polinomio con solo 10 raíces reales y 10 raíces complejas con partes imaginarias mayores que 0,6.
El método de Halley y el método de Laguerre están estrechamente relacionados con el método de Newton . Ambos utilizan el polinomio y sus dos primeras derivaciones para un proceso iterativo que tiene una convergencia cúbica . Combinando dos pasos consecutivos de estos métodos en una única prueba, se obtiene una tasa de convergencia de 9, a costa de 6 evaluaciones de polinomios (con la regla de Horner). Por otro lado, la combinación de tres pasos del método de Newton da una tasa de convergencia de 8 a costa del mismo número de evaluaciones de polinomios. Esto da una ligera ventaja a estos métodos (menos clara para el método de Laguerre, ya que se debe calcular una raíz cuadrada en cada paso).
Al aplicar estos métodos a polinomios con coeficientes reales y puntos de partida reales, el método de Newton y el de Halley se mantienen dentro de la línea de números reales. Hay que elegir puntos de partida complejos para encontrar raíces complejas. Por el contrario, el método de Laguerre con una raíz cuadrada en su evaluación abandonará el eje real por sí solo.
Si el polinomio dado solo tiene coeficientes reales, se puede evitar realizar cálculos con números complejos. Para ello, hay que encontrar factores cuadráticos para pares de raíces complejas conjugadas. La aplicación del método multidimensional de Newton a esta tarea da como resultado el método de Bairstow .
La variante real del algoritmo Jenkins-Traub es una mejora de este método.
El método simple de Durand-Kerner y el ligeramente más complicado de Aberth encuentran simultáneamente todas las raíces utilizando únicamente aritmética de números complejos . Los algoritmos acelerados para la evaluación de múltiples puntos y la interpolación similares a la transformada rápida de Fourier pueden ayudar a acelerarlos para grados grandes del polinomio. Es aconsejable elegir un conjunto de puntos iniciales asimétricos, pero distribuidos uniformemente. La implementación de este método en el software gratuito MPSolve es una referencia por su eficiencia y precisión.
Otro método con este estilo es el método de Dandelin-Gräffe (a veces también atribuido a Lobachevsky ), que utiliza transformaciones polinómicas para elevar al cuadrado repetida e implícitamente las raíces. Esto magnifica enormemente las varianzas en las raíces. Aplicando las fórmulas de Viète , se obtienen aproximaciones fáciles para el módulo de las raíces y, con un poco más de esfuerzo, para las raíces mismas.
Existen varias pruebas rápidas que indican si un segmento de la línea real o una región del plano complejo no contiene raíces. Al acotar el módulo de las raíces y subdividir recursivamente la región inicial indicada por estos límites, se pueden aislar pequeñas regiones que pueden contener raíces y luego aplicar otros métodos para localizarlas con exactitud.
Todos estos métodos implican encontrar los coeficientes de versiones desplazadas y escaladas del polinomio. Para grados grandes, los métodos acelerados basados en FFT resultan viables.
Para raíces reales, consulte las siguientes secciones.
El algoritmo de Lehmer-Schur utiliza la prueba de Schur-Cohn para círculos; una variante, el algoritmo de bisección global de Wilf, utiliza un cálculo de número sinuoso para regiones rectangulares en el plano complejo.
El método de división del círculo utiliza transformaciones polinómicas basadas en FFT para encontrar factores de alto grado correspondientes a grupos de raíces. La precisión de la factorización se maximiza utilizando una iteración de tipo Newton. Este método es útil para encontrar las raíces de polinomios de alto grado con precisión arbitraria; tiene una complejidad casi óptima en este contexto. [ cita requerida ]
Encontrar las raíces reales de un polinomio con coeficientes reales es un problema que ha recibido mucha atención desde principios del siglo XIX y que sigue siendo un campo activo de investigación. La mayoría de los algoritmos de búsqueda de raíces pueden encontrar algunas raíces reales, pero no pueden certificar haber encontrado todas las raíces. Los métodos para encontrar todas las raíces complejas, como el método de Aberth , pueden proporcionar las raíces reales. Sin embargo, debido a la inestabilidad numérica de los polinomios (véase el polinomio de Wilkinson ), pueden necesitar una aritmética de precisión arbitraria para decidir qué raíces son reales. Además, calculan todas las raíces complejas cuando solo unas pocas son reales.
De ello se deduce que la forma estándar de calcular raíces reales es calcular primero intervalos disjuntos, llamados intervalos de aislamiento , de modo que cada uno contenga exactamente una raíz real y juntos contengan todas las raíces. Este cálculo se llama aislamiento de raíces reales . Al tener un intervalo de aislamiento, se pueden utilizar métodos numéricos rápidos, como el método de Newton , para mejorar la precisión del resultado.
El algoritmo más antiguo y completo para el aislamiento de raíces reales resulta del teorema de Sturm . Sin embargo, parece ser mucho menos eficiente que los métodos basados en la regla de los signos de Descartes y el teorema de Vincent . Estos métodos se dividen en dos clases principales, una que utiliza fracciones continuas y la otra que utiliza bisección. Ambos métodos han mejorado drásticamente desde principios del siglo XXI. Con estas mejoras alcanzan una complejidad computacional similar a la de los mejores algoritmos para calcular todas las raíces (incluso cuando todas las raíces son reales).
Estos algoritmos se han implementado y están disponibles en Mathematica (método de fracciones continuas) y Maple (método de bisección). Ambas implementaciones pueden encontrar rutinariamente las raíces reales de polinomios de grado superior a 1000.
Las raíces múltiples son altamente sensibles, se sabe que están mal condicionadas e imprecisas en el cálculo numérico en general. Un método de Zhonggang Zeng (2004), implementado como un paquete MATLAB , calcula raíces múltiples y multiplicidades correspondientes de un polinomio con precisión incluso si los coeficientes son inexactos. [3] [4] [5]
El método se puede resumir en dos pasos. Sea el polinomio dado. El primer paso determina la estructura de multiplicidad aplicando factorización libre de cuadrados con un algoritmo numérico de máximo común divisor. [5] Esto permite escribir como
donde son las multiplicidades de las raíces distintas. Esta ecuación es un sistema sobredeterminado para tener variables en ecuaciones que coinciden con coeficientes (el coeficiente principal no es una variable). La solución de mínimos cuadrados ya no está mal condicionada en la mayoría de los casos. El segundo paso aplica el algoritmo de Gauss-Newton para resolver el sistema sobredeterminado para las raíces distintas.
La sensibilidad de raíces múltiples se puede regularizar debido a una propiedad geométrica de raíces múltiples descubierta por William Kahan (1972) y el modelo de sistema sobredeterminado mantiene las multiplicidades .
Para los polinomios cuyos coeficientes se dan exactamente como números enteros o racionales , existe un método eficiente para factorizarlos en factores que solo tienen raíces simples y cuyos coeficientes también se dan exactamente. Este método, llamado factorización libre de cuadrados , se basa en que las raíces múltiples de un polinomio son las raíces del máximo común divisor del polinomio y su derivada.
La factorización libre de cuadrados de un polinomio p es una factorización donde cada uno es 1 o un polinomio sin raíces múltiples, y dos diferentes no tienen ninguna raíz común.
Un método eficiente para calcular esta factorización es el algoritmo de Yun .
Zeng, Zhonggang (2004). "Algoritmo 835: MultRoot – Un paquete de Matlab para calcular raíces y multiplicidades polinómicas". ACM Transactions on Mathematical Software . 30 : 218-236. doi :10.1145/992200.992209. S2CID 18188044.
