El algoritmo de Jenkins-Traub para polinomios ceros es un método rápido de búsqueda de raíces polinomiales iterativo convergente globalmente publicado en 1970 por Michael A. Jenkins y Joseph F. Traub . Dieron dos variantes, una para polinomios generales con coeficientes complejos, comúnmente conocida como el algoritmo "CPOLY", y una variante más complicada para el caso especial de polinomios con coeficientes reales, comúnmente conocida como el algoritmo "RPOLY". Este último es "prácticamente un estándar en buscadores de raíces polinomiales de caja negra". [1]
Este artículo describe la variante compleja. Dado un polinomio P ,
con coeficientes complejos calcula aproximaciones a los n cerosde P ( z ), uno a la vez en orden de magnitud aproximadamente creciente. Después de calcular cada raíz, su factor lineal se elimina del polinomio. El uso de esta deflación garantiza que cada raíz se calcule solo una vez y que se encuentren todas las raíces.
La variante real sigue el mismo patrón, pero calcula dos raíces a la vez, ya sea dos raíces reales o un par de raíces complejas conjugadas. Al evitar la aritmética compleja, la variante real puede ser más rápida (por un factor de 4) que la variante compleja. El algoritmo de Jenkins-Traub ha estimulado una considerable investigación sobre teoría y software para métodos de este tipo.
Descripción general
El algoritmo de Jenkins-Traub calcula todas las raíces de un polinomio con coeficientes complejos. El algoritmo comienza verificando el polinomio en busca de raíces muy grandes o muy pequeñas. Si es necesario, los coeficientes se reescalan mediante un reajuste de la variable. En el algoritmo, las raíces adecuadas se encuentran una por una y generalmente en tamaño creciente. Después de encontrar cada raíz, el polinomio se desinfla dividiendo el factor lineal correspondiente. De hecho, la factorización del polinomio en el factor lineal y el polinomio desinflado restante ya es el resultado del procedimiento de búsqueda de raíces. El procedimiento de búsqueda de raíces tiene tres etapas que corresponden a diferentes variantes de la iteración de potencia inversa . Vea Jenkins y Traub . [2] También se puede encontrar una descripción en Ralston y Rabinowitz [3] p. 383. El algoritmo es similar en espíritu al algoritmo de dos etapas estudiado por Traub. [4]
Procedimiento de búsqueda de raíces
Comenzando con el polinomio actual P ( X ) de grado n , se calcula la raíz más pequeña de P (x) . Para ello, se construye una secuencia de los denominados polinomios H. Estos polinomios son todos de grado n - 1 y se supone que convergen al factor de P ( X ) que contiene todas las raíces restantes. La secuencia de polinomios H ocurre en dos variantes, una variante no normalizada que permite conocimientos teóricos fáciles y una variante normalizada de polinomios que mantienen los coeficientes en un rango numéricamente sensible.
La construcción de los polinomios H depende de una secuencia de números complejos llamados turnos. Estos mismos cambios dependen, al menos en la tercera etapa, de los polinomios H previos . Los polinomios H se definen como la solución a la recursividad implícita
- y
Una solución directa a esta ecuación implícita es
donde la división polinomial es exacta.
Algorítmicamente, se usaría, por ejemplo, el esquema de Horner o la regla de Ruffini para evaluar los polinomios eny obtener los cocientes al mismo tiempo. Con los cocientes resultantes p ( X ) yh ( X ) como resultados intermedios, el siguiente polinomio H se obtiene como
Dado que el coeficiente de grado más alto se obtiene de P (X) , el coeficiente principal de es . Si esto se divide, el polinomio H normalizado es
Etapa uno: proceso sin cambios
Para colocar . Por lo general, se elige M = 5 para polinomios de grados moderados hasta n = 50. Esta etapa no es necesaria solo por consideraciones teóricas, pero es útil en la práctica. Destaca en los polinomios H el cofactor (del factor lineal) de la raíz más pequeña.
Etapa dos: proceso de turno fijo
El desplazamiento para esta etapa se determina como un punto cercano a la raíz más pequeña del polinomio. Está ubicado casi al azar en el círculo con el radio de la raíz interior, que a su vez se estima como la solución positiva de la ecuación.
Dado que el lado izquierdo es una función convexa y aumenta monótonamente desde cero hasta el infinito, esta ecuación es fácil de resolver, por ejemplo, mediante el método de Newton .
Ahora elige en el círculo de este radio. La secuencia de polinomios, , se genera con el valor de turno fijo . Durante esta iteración, la aproximación actual para la raíz
se rastrea. La segunda etapa se termina con éxito si las condiciones
- y
se cumplen simultáneamente. Si no hubo éxito después de varias iteraciones, se intentará con un punto aleatorio diferente en el círculo. Por lo general, se usan 9 iteraciones para polinomios de grado moderado, con una estrategia de duplicación para el caso de fallas múltiples.
Etapa tres: proceso de cambio variable
La ahora se generan usando los cambios variables que son generados por
siendo la última raíz estimada de la segunda etapa y
- dónde es el polinomio H normalizado , es decir dividido por su coeficiente principal.
Si el tamaño del paso en la etapa tres no cae lo suficientemente rápido a cero, entonces la etapa dos se reinicia usando un punto aleatorio diferente. Si esto no tiene éxito después de un pequeño número de reinicios, el número de pasos en la etapa dos se duplica.
Convergencia
Se puede demostrar que, siempre que L se elige suficientemente grande, es λ siempre converge a una raíz de P .
El algoritmo converge para cualquier distribución de raíces, pero puede fallar al encontrar todas las raíces del polinomio. Además, la convergencia es ligeramente más rápida que la convergencia cuadrática de la iteración de Newton-Raphson, sin embargo, utiliza al menos el doble de operaciones por paso.
¿Qué le da al algoritmo su poder?
Comparar con la iteración de Newton-Raphson
La iteración usa la P dada y. En contraste, la tercera etapa de Jenkins-Traub
es precisamente una iteración de Newton-Raphson realizada sobre ciertas funciones racionales . Más precisamente, Newton-Raphson se realiza en una secuencia de funciones racionales
Para suficientemente largo,
está tan cerca como se desea de un polinomio de primer grado
dónde es uno de los ceros de . Aunque la Etapa 3 es precisamente una iteración de Newton-Raphson, no se realiza la diferenciación.
Análisis de los polinomios H
Dejar ser las raíces de P ( X ). Los llamados factores de Lagrange de P (X) son los cofactores de estas raíces,
Si todas las raíces son diferentes, entonces los factores de Lagrange forman una base del espacio de polinomios de grado como máximo n - 1. Mediante el análisis del procedimiento de recursividad se encuentra que los polinomios H tienen la representación de coordenadas
Cada factor de Lagrange tiene un coeficiente principal 1, de modo que el coeficiente principal de los polinomios H es la suma de los coeficientes. Por tanto, los polinomios H normalizados son
Órdenes de convergencia
Si la condición se mantiene para casi todas las iteraciones, los polinomios H normalizados convergerán al menos geométricamente hacia .
Con la condición de que
uno obtiene las estimaciones asintóticas para
- Nivel 1:
- para la etapa 2, si s está lo suficientemente cerca de:
- y
- y para la etapa 3:
- y
- dando lugar a un orden de convergencia superior al cuadrático , dónde es la proporción áurea .
Interpretación como iteración de potencia inversa
Todas las etapas del algoritmo complejo de Jenkins-Traub pueden representarse como el problema de álgebra lineal para determinar los valores propios de una matriz especial. Esta matriz es la representación de coordenadas de un mapa lineal en el espacio n -dimensional de polinomios de grado n - 1 o menos. La idea principal de este mapa es interpretar la factorización
con una raíz y el factor restante de grado n - 1 como la ecuación del vector propio para la multiplicación con la variable X , seguido del cálculo del resto con el divisor P ( X ),
Esto mapea polinomios de grado como máximo n - 1 a polinomios de grado como máximo n - 1. Los valores propios de este mapa son las raíces de P ( X ), ya que la ecuación del vector propio dice
lo que implica que , es decir, es un factor lineal de P ( X ). En la base monomial el mapa linealestá representado por una matriz compañera del polinomio P , como
la matriz de coeficientes resultante es
A esta matriz se le aplica la iteración de potencia inversa en las tres variantes de no desplazamiento, desplazamiento constante y desplazamiento de Rayleigh generalizado en las tres etapas del algoritmo. Es más eficiente realizar las operaciones de álgebra lineal en aritmética polinomial y no mediante operaciones matriciales, sin embargo, las propiedades de la iteración de potencia inversa siguen siendo las mismas.
Coeficientes reales
El algoritmo de Jenkins-Traub descrito anteriormente funciona para polinomios con coeficientes complejos. Los mismos autores también crearon un algoritmo de tres etapas para polinomios con coeficientes reales. Consulte Jenkins y Traub Un algoritmo de tres etapas para polinomios reales mediante iteración cuadrática . [5] El algoritmo encuentra un factor lineal o cuadrático que funciona completamente en aritmética real. Si los algoritmos complejos y reales se aplican al mismo polinomio real, el algoritmo real es aproximadamente cuatro veces más rápido. El algoritmo real siempre converge y la tasa de convergencia es mayor que la de segundo orden.
Una conexión con el algoritmo QR desplazado
Existe una conexión sorprendente con el algoritmo QR desplazado para calcular los valores propios de la matriz. Consulte Dekker y Traub El algoritmo QR modificado para matrices hermitianas . [6] Nuevamente, los cambios pueden verse como una iteración de Newton-Raphson en una secuencia de funciones racionales que convergen en un polinomio de primer grado.
Software y pruebas
El software para el algoritmo Jenkins-Traub se publicó como Jenkins y Traub Algorithm 419: Zeros of a Complex Polynomial . [7] El software para el algoritmo real fue publicado como Jenkins Algorithm 493: Zeros of a Real Polynomial . [8]
Los métodos han sido probados extensamente por muchas personas. Como se predijo, disfrutan de una convergencia más rápida que la cuadrática para todas las distribuciones de ceros.
Sin embargo, hay polinomios que pueden causar pérdida de precisión, como se ilustra en el siguiente ejemplo. El polinomio tiene todos sus ceros en dos semicírculos de diferentes radios. Wilkinson recomienda que es deseable para una deflación estable que se calculen primero los ceros más pequeños. Los cambios de la segunda etapa se eligen de modo que los ceros en el semicírculo más pequeño se encuentren primero. Después de la deflación, se sabe que el polinomio con los ceros en el semicírculo está mal acondicionado si el grado es grande; ver Wilkinson, [9] p. 64. El polinomio original era de grado 60 y sufría una grave inestabilidad deflacionista.
Referencias
- ^ Prensa, WH, Teukolsky, SA, Vetterling, WT y Flannery, BP (2007), Recetas numéricas: El arte de la informática científica, 3ª ed., Cambridge University Press, página 470.
- ^ Jenkins, MA y Traub, JF (1970), una iteración de cambio de variables de tres etapas para ceros polinomiales y su relación con la iteración de Rayleigh generalizada , número. Matemáticas. 14, 252-263.
- ^ Ralston, A. y Rabinowitz, P. (1978), Un primer curso de análisis numérico, 2ª ed., McGraw-Hill, Nueva York.
- ^ Traub, JF (1966), una clase de funciones de iteración globalmente convergentes para la solución de ecuaciones polinomiales , matemáticas. Comp., 20 (93), 113-138.
- ^ Jenkins, MA y Traub, JF (1970), Un algoritmo de tres etapas para polinomios reales mediante iteración cuadrática , SIAM J. Numer. Anal., 7 (4), 545-566.
- ^ Dekker, TJ y Traub, JF (1971), El algoritmo QR desplazado para matrices hermitianas , Lin. Álgebra Appl., 4 (2), 137-154.
- ^ Jenkins, MA y Traub, JF (1972), Algoritmo 419: Ceros de un polinomio complejo , Comm. ACM, 15, 97–99.
- ^ Jenkins, MA (1975), Algoritmo 493: Ceros de un polinomio real , ACM TOMS, 1, 178-189.
- ^ Wilkinson, JH (1963), Redondeo de errores en procesos algebraicos, Prentice Hall, Englewood Cliffs, Nueva Jersey
enlaces externos
- Una aplicación de Windows descargable gratuita que utiliza el método Jenkins-Traub para polinomios con coeficientes reales y complejos.
- RPoly ++ Una implementación C ++ optimizada para SSE del algoritmo RPOLY.