Método de Durand-Kerner

En análisis numérico , el método de Weierstrass o método de Durand-Kerner , descubierto por Karl Weierstrass en 1891 y redescubierto independientemente por Durand en 1960 y Kerner en 1966, es un algoritmo de búsqueda de raíces para resolver ecuaciones polinómicas . ^{[ 1 ]} En otras palabras, el método se puede utilizar para resolver numéricamente la ecuación f ( x ) = 0, donde f es un polinomio dado, que se puede tomar escalado de manera que el coeficiente principal sea 1. 

Esta explicación considera ecuaciones de cuarto grado . Se puede generalizar fácilmente a otros grados.

Sea f el polinomio definido por

${\displaystyle f(x)=x^{4}+ax^{3}+bx^{2}+cx+d}$

para todo x . Los números conocidos a , b , c , d son los coeficientes .

Sean los números (potencialmente complejos) P , Q , R , S las raíces de este polinomio f . Entonces

${\displaystyle f(x)=(xP)(xQ)(xR)(xS)}$

para todo x . Se puede aislar el valor P de esta ecuación:

${\displaystyle P=x-{\frac {f(x)}{(xQ)(xR)(xS)}}.}$

Entonces, si se usa como una iteración de punto fijo

${\displaystyle x_{1}:=x_{0}-{\frac {f(x_{0})}{(x_{0}-Q)(x_{0}-R)(x_{0}-S)}},}$

Es fuertemente estable en el sentido de que cada punto inicial x ₀ ≠ Q , R , S proporciona después de una iteración la raíz P = x ₁ . Además, si se reemplazan los ceros Q , R y S por aproximaciones q ≈ Q , r ≈ R , s ≈ S , tales que q , r , s no son iguales a P , entonces P sigue siendo un punto fijo de la iteración de punto fijo perturbada.

${\displaystyle x_{k+1}:=x_{k}-{\frac {f(x_{k})}{(x_{k}-q)(x_{k}-r)(x_{k}-s)}},}$

desde

${\displaystyle P-{\frac {f(P)}{(Pq)(Pr)(Ps)}}=P-0=P.}$

Nótese que el denominador sigue siendo distinto de cero. Esta iteración de punto fijo es una aplicación de contracción para x alrededor de P.

La clave del método ahora consiste en combinar la iteración de punto fijo para P con iteraciones similares para Q , R y S en una iteración simultánea para todas las raíces.

Inicializar p , q , r , s :

p ₀  := (0.4 + 0.9 i ) ⁰ ,

q ₀  := (0.4 + 0.9 i ) ¹ ,

r ₀  := (0.4 + 0.9 i ) ² ,

s ₀  := (0.4 + 0.9 i ) ³ .

No hay nada especial en elegir 0,4  +  0,9 i excepto que no es ni un número real ni una raíz de la unidad .

Realiza las sustituciones para n = 1, 2, 3, ...:

${\displaystyle p_{n}=p_{n-1}-{\frac {f(p_{n-1})}{(p_{n-1}-q_{n-1})(p_{n-1}-r_{n-1})(p_{n-1}-s_{n-1})}},}$

${\displaystyle q_{n}=q_{n-1}-{\frac {f(q_{n-1})}{(q_{n-1}-p_{n})(q_{n-1}-r_{n-1})(q_{n-1}-s_{n-1})}},}$

${\displaystyle r_{n}=r_{n-1}-{\frac {f(r_{n-1})}{(r_{n-1}-p_{n})(r_{n-1}-q_{n})(r_{n-1}-s_{n-1})}},}$

${\displaystyle s_{n}=s_{n-1}-{\frac {f(s_{n-1})}{(s_{n-1}-p_{n})(s_{n-1}-q_{n})(s_{n-1}-r_{n})}}.}$

Repita el proceso hasta que los números p , q , r , s dejen de variar con respecto a la precisión deseada. Entonces, tendrán los valores P , Q , R , S en algún orden y con la precisión elegida. De esta forma, el problema queda resuelto.

Tenga en cuenta que debe utilizarse aritmética de números complejos y que las raíces se encuentran simultáneamente, en lugar de una a una.

Este procedimiento iterativo, al igual que el método de Gauss-Seidel para ecuaciones lineales, calcula un número a la vez basándose en los números ya calculados. Una variante de este procedimiento, como el método de Jacobi , calcula un vector de aproximaciones de raíces a la vez. Ambas variantes son algoritmos eficaces para la búsqueda de raíces.

También se podrían elegir los valores iniciales para p , q , r , s mediante algún otro procedimiento, incluso aleatoriamente, pero de tal manera que

• están dentro de un círculo no demasiado grande que también contiene las raíces de f ( x ), por ejemplo, el círculo alrededor del origen con radio${\displaystyle 1+\max {\big (}|a|,|b|,|c|,|d|{\big )}}$(donde 1, a , b , c , d son los coeficientes de f ( x ))

y eso

• no están demasiado cerca el uno del otro,

lo cual puede convertirse en una preocupación cada vez mayor a medida que aumenta el grado del polinomio.

Si los coeficientes son reales y el polinomio tiene grado impar, entonces debe tener al menos una raíz real. Para hallarla, use un valor real de p₀ como estimación inicial y haga que _q₀ y r₀ , etc., _{sean pares} complejos conjugados . Entonces _, la iteración conservará estas propiedades; es decir, pₙ siempre será real, y qₙ y rₙ _, etc., siempre serán conjugados. De esta manera, pₙ convergerá a una raíz real P. Alternativamente, _{haga que todas las estimaciones iniciales sean} reales; permanecerán así.

Este ejemplo proviene de la referencia Jacoby (1992). La ecuación resuelta es x³ ⁻ 3x² + 3x − 5 = 0. Las primeras 4 iteraciones mueven p , q , r aparentemente de forma caótica, pero luego las raíces se ubican con una ^precisión de 1 decimal. Después de la iteración número 5, tenemos 4 decimales correctos, y la iteración número 6 confirma que las raíces calculadas son fijas. Este comportamiento general es característico del método. Nótese también que, en este ejemplo, las raíces se utilizan tan pronto como se calculan en cada iteración. En otras palabras, el cálculo de cada segunda columna utiliza el valor de las columnas calculadas anteriormente.

Nótese que la ecuación tiene una raíz real y un par de raíces complejas conjugadas, y que la suma de las raíces es  3.

Para cada n -tupla de números complejos, existe exactamente un polinomio mónico de grado n que los tiene como raíces (manteniendo las multiplicidades). Este polinomio se obtiene multiplicando todos los factores lineales correspondientes, es decir:

${\displaystyle g_{\vec {z}}(X)=(X-z_{1})\cdots (X-z_{n}).}$

Este polinomio tiene coeficientes que dependen de los ceros prescritos,

${\displaystyle g_{\vec {z}}(X)=X^{n}+g_{n-1}({\vec {z}})X^{n-1}+\cdots +g_{0}({\vec {z}}).}$

Esos coeficientes son, salvo un signo, los polinomios simétricos elementales.${\displaystyle \alpha _{1}({\vec {z}}),\dots ,\alpha _{n}({\vec {z}})}$de grados 1,...,n .

Para hallar todas las raíces de un polinomio dado${\displaystyle f(X)=X^{n}+c_{n-1}X^{n-1}+\cdots +c_{0}}$con vector de coeficientes${\displaystyle (c_{n-1},\dots ,c_{0})}$Simultáneamente, ahora es lo mismo que encontrar un vector solución para el sistema de Vieta.

${\displaystyle {\begin{matrix}c_{0}&=&g_{0}({\vec {z}})&=&(-1)^{n}\alpha _{n}({\vec {z}})&=&(-1)^{n}z_{1}\cdots z_{n}\\c_{1}&=&g_{1}({\vec {z}})&=&(-1)^{n-1}\alpha _{n-1}({\vec {z}})\\&\vdots &\\c_{n-1}&=&g_{n-1}({\vec {z}})&=&-\alpha _{1}({\vec {z}})&=&-(z_{1}+z_{2}+\cdots +z_{n}).\end{matrix}}}$

El método de Durand-Kerner se obtiene como el método de Newton multidimensional aplicado a este sistema. Es algebraicamente más cómodo tratar esas identidades de coeficientes como la identidad de los polinomios correspondientes,${\displaystyle g_{\vec {z}}(X)=f(X)}$En el método de Newton se busca, dado algún vector inicial${\displaystyle {\vec {z}}}$, para un vector de incremento${\displaystyle {\vec {w}}}$de tal manera que${\displaystyle g_{{\vec {z}}+{\vec {w}}}(X)=f(X)}$se satisface hasta los términos de segundo orden y superiores en el incremento. Para ello se resuelve la identidad.

${\displaystyle f(X)-g_{\vec {z}}(X)=\sum _{k=1}^{n}{\frac {\partial g_{\vec {z}}(X)}{\partial z_{k}}}w_{k}=-\sum _{k=1}^{n}w_{k}\prod _{j\neq k}(X-z_{j}).}$

Si los números${\displaystyle z_{1},\dots ,z_{n}}$Si son diferentes por pares, entonces los polinomios en los términos del lado derecho forman una base del espacio n -dimensional.${\displaystyle \mathbb {C} [X]_{n-1}}$de polinomios con grado máximo n − 1. Por lo tanto, una solución  ${\displaystyle {\vec {w}}}$En este caso existe la ecuación del incremento. Las coordenadas del incremento${\displaystyle {\vec {w}}}$se obtienen simplemente evaluando la ecuación de incremento

${\displaystyle -\sum _{k=1}^{n}w_{k}\prod _{j\neq k}(X-z_{j})=f(X)-\prod _{j=1}^{n}(X-z_{j})}$

en los puntos${\displaystyle X=z_{k}}$, lo que resulta en

${\displaystyle -w_{k}\prod _{j\neq k}(z_{k}-z_{j})=-w_{k}g_{\vec {z}}'(z_{k})=f(z_{k})}$, eso es${\displaystyle w_{k}=-{\frac {f(z_{k})}{\prod _{j\neq k}(z_{k}-z_{j})}}.}$

En el anillo cociente (álgebra) de clases de residuos módulo ƒ ( X ), la multiplicación por X define un endomorfismo cuyos ceros de ƒ ( X ) son autovalores con las multiplicidades correspondientes. Al elegir una base, el operador de multiplicación se representa mediante su matriz de coeficientes A , la matriz compañera de ƒ ( X ) para dicha base.

Dado que todo polinomio puede reducirse módulo ƒ ( X ) a un polinomio de grado n − 1 o inferior, el espacio de clases de residuos puede identificarse con el espacio de polinomios de grado acotado por n − 1. Una base específica del problema puede tomarse de la interpolación de Lagrange como el conjunto de n polinomios.    

${\displaystyle b_{k}(X)=\prod _{1\leq j\leq n,\;j\neq k}(X-z_{j}),\quad k=1,\dots ,n,}$

dónde${\displaystyle z_{1},\dots ,z_{n}\in \mathbb {C} }$son números complejos distintos entre sí. Nótese que las funciones núcleo para la interpolación de Lagrange son${\displaystyle L_{k}(X)={\frac {b_{k}(X)}{b_{k}(z_{k})}}}$.

Para el operador de multiplicación aplicado a los polinomios base se obtiene a partir de la interpolación de Lagrange

dónde${\displaystyle w_{j}=-{\frac {f(z_{j})}{b_{j}(z_{j})}}}$Son de nuevo las actualizaciones de Weierstrass.

Por lo tanto, la matriz compañera de ƒ ( X ) es

${\displaystyle A=\mathrm {diag} (z_{1},\dots ,z_{n})+{\begin{pmatrix}1\\\vdots \\1\end{pmatrix}}\cdot \left(w_{1},\dots ,w_{n}\right).}$

Del caso de la matriz transpuesta del teorema del círculo de Gershgorin se deduce que todos los autovalores de A , es decir, todas las raíces de ƒ ( X ), están contenidos en la unión de los discos.${\displaystyle D(a_{k,k},r_{k})}$con un radio${\displaystyle r_{k}=\sum _{j\neq k}{\big |}a_{j,k}{\big |}}$.

Aquí uno tiene${\displaystyle a_{k,k}=z_{k}+w_{k}}$, por lo que los centros son las siguientes iteraciones de la iteración de Weierstrass y los radios${\displaystyle r_{k}=(n-1)\left|w_{k}\right|}$que son múltiplos de las actualizaciones de Weierstrass. Si las raíces de ƒ ( X ) están todas bien aisladas (en relación con la precisión computacional) y los puntos${\displaystyle z_{1},\dots ,z_{n}\in \mathbb {C} }$Si las aproximaciones a estas raíces son suficientemente cercanas, entonces todos los discos serán disjuntos, de modo que cada uno contendrá exactamente un cero. Los puntos medios de los círculos serán mejores aproximaciones de los ceros.

Cada matriz conjugada${\displaystyle TAT^{-1}}$de A es también una matriz compañera de ƒ ( X ). Elegir T como matriz diagonal deja la estructura de A invariante. La raíz cercana a${\displaystyle z_{k}}$está contenido en cualquier círculo aislado con centro${\displaystyle z_{k}}$independientemente de T. Elegir la matriz diagonal óptima T para cada índice da como resultado mejores estimaciones (ver referencia Petkovic et al. 1995).

La conexión entre la expansión en serie de Taylor y el método de Newton sugiere que la distancia desde${\displaystyle z_{k}+w_{k}}$a la raíz correspondiente es del orden${\displaystyle O{\big (}|w_{k}|^{2}{\big )}}$, si la raíz está bien aislada de raíces cercanas y la aproximación está suficientemente cerca de la raíz. Así, una vez que la aproximación está cerca, el método de Newton converge cuadráticamente ; es decir, el error se eleva al cuadrado con cada paso (lo que reducirá enormemente el error una vez que sea menor que 1). En el caso del método de Durand-Kerner, la convergencia es cuadrática si el vector${\displaystyle {\vec {z}}=(z_{1},\dots ,z_{n})}$está cerca de alguna permutación del vector de las raíces de f .

Para la conclusión de convergencia lineal existe un resultado más específico (véase la referencia Petkovic et al. 1995). Si el vector inicial${\displaystyle {\vec {z}}}$y su vector de actualizaciones de Weierstrass${\displaystyle {\vec {w}}=(w_{1},\dots ,w_{n})}$satisface la desigualdad

${\displaystyle \max _{1\leq k\leq n}|w_{k}|\leq {\frac {1}{5n}}\min _{1\leq j<k\leq n}|z_{k}-z_{j}|,}$

entonces esta desigualdad también se cumple para todas las iteraciones, todos los discos de inclusión${\displaystyle D{\big (}z_{k}+w_{k},(n-1)|w_{k}|{\big )}}$son disjuntos y se cumple la convergencia lineal con un factor de contracción de 1/2. Además, los discos de inclusión pueden elegirse en este caso como

${\displaystyle D\left(z_{k}+w_{k},{\tfrac {1}{4}}|w_{k}|\right),\quad k=1,\dots ,n,}$

cada una contiene exactamente un cero de f .

El método de Weierstrass/Durand-Kerner no es generalmente convergente: en otras palabras, no es cierto que para cada polinomio, el conjunto de vectores iniciales que finalmente converge a raíces sea abierto y denso. De hecho, existen conjuntos abiertos de polinomios cuyos conjuntos de vectores iniciales también son abiertos y convergen a ciclos periódicos distintos de las raíces (véase Reinke et al.).

• Weierstrass, Karl (1891). "Neuer Beweis des Satzes, dass jede ganze racionale Function einer Veränderlichen dargestellt werden kann als ein Product aus linearen Functionen derselben Veränderlichen" . Sitzungsberichte der königlich preussischen Akademie der Wissenschaften zu Berlin . Archivado desde el original el 2 de noviembre de 2013 . Consultado el 31 de octubre de 2013 .
• Durand, E. (1960). "Ecuaciones del tipo F ( x )  =  0: Racines d'un polinome". En Masón; et  al. (eds.). Soluciones Numériques des Equations Algébriques . vol.  1.
• Kerner, Immo O. (1966). "Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen". Matemática numérica . 8 (3): 290– 294. doi : 10.1007/BF02162564 . S2CID 115307022 . 
• Prešić, Marica (1980). "Un teorema de convergencia para un método de determinación simultánea de todos los ceros de un polinomio" (PDF) . Publications de l'Institut Mathématique . Nouvelle Série. 28 (42): 158– 168.
• Petkovic, MS, Carstensen, C. y Trajkovic, M. (1995). "Fórmula de Weierstrass y métodos para encontrar ceros". Numerische Mathematik . 69 (3): 353– 372. CiteSeerX 10.1.1.53.7516 . doi : 10.1007/s002110050097 . S2CID 18594004 .  {{cite journal}}: CS1 maint: varios nombres: lista de autores ( enlace )
• Bo Jacoby, Nulpunkter for polynomier , CAE-nyt (una publicación periódica del Dansk CAE Gruppe [Grupo CAE danés]), 1988.
• Agnethe Knudsen, Numeriske Metoder (notas de clase), Københavns Teknikum.
• Bo Jacoby, Numerisk løsning af ligninger , Bygningsstatiske meddelelser (Publicado por la Sociedad Danesa de Ciencias e Ingeniería Estructurales) volumen 63 no. 3–4, 1992, págs.  83–105.
• Gourdon, Xavier (1996). Combinación, algoritmos y geometría de polinomios . París: École Polytechnique. Archivado desde el original el 28 de octubre de 2006 . Consultado el 22 de agosto de 2006 .
• Victor Pan (mayo de 2002): Búsqueda de raíces de polinomios univariados con menor precisión computacional y mayores tasas de convergencia . Informe técnico, Universidad de la Ciudad de Nueva York.
• Neumaier, Arnold (2003). "Clústeres de ceros de polinomios que encierran" . Journal of Computational and Applied Mathematics . 156 (2): 389– 401. Bibcode : 2003JCoAM.156..389N . doi : 10.1016/S0377-0427(03)00380-7 .
• Jan Verschelde, El método de Weierstrass (también conocido como método Durand-Kerner) , 2003.
• Bernhard Reinke, Dierk Schleicher y Michael Stoll, " El buscador de raíces de Weierstrass no es generalmente convergente ", 2020
  • Bernhard Reinke, Dierk Schleicher y Michael Stoll: "El buscador de raíces de Weierstrass-Durand-Kerner no es generalmente convergente", Math. comp. vol.92 (2023), págs.839-866. DOI: https://doi.org/10.1090/mcom/3783 .

• Ada Generic_Roots usando el método Durand–Kerner (archivo) — unaimplementación de código abierto en Ada
• Raíces polinomiales : una implementación de código abierto en Java.
• Extracción de raíces de polinomios  : El método de Durand-Kerner — incluye unademostración en un applet de Java