Encontrar todas las raíces de un polinomio

Question

¿Es posible, para un polinomio arbitrario en una variable con coeficientes enteros, determinar las raíces del polinomio en el campo complejo con una precisión arbitraria? Cuando estaba investigando esto, encontré algunos artículos sobre la continuación de homotopía que parecen resolver este problema (para las soluciones reales al menos), ¿es eso correcto? ¿O hay restricciones para que la continuación de homotopía funcione? ¿Es necesario que la región de la solución esté acotada?

Answer 1

Este argumento es problemático; véase el comentario de Andrej Bauer más abajo.

---

Claro. No tengo ni idea de cómo es un algoritmo eficiente, pero como sólo has preguntado si es posible te ofreceré uno terrible.

Lema: Dejemos que $f(z) = z^n + a_{n-1} z^{n-1} + \cdots + a_0$ sea un polinomio complejo y que $R = \max(1, |a_{n-1}| + \cdots + |a_0|)$ . Entonces todas las raíces de $f$ se encuentran en el círculo de radio $R$ centrado en el origen.

Prueba. Si $|z| > R$ entonces $|z|^n > R |z|^{n-1} \ge |a_{n-1} z^{n-1}| + \cdots + |a_0|$ por lo que por la desigualdad del triángulo no hay tal $z$ es una raíz.

Ahora subdivide el disco de radio $R$ en, digamos, una malla de cuadrados de lado $\varepsilon > 0$ y evaluar el polinomio en todos los puntos de la malla. A medida que el tamaño de la malla tiende a cero, encontrarás puntos que se aproximan a los ceros con una precisión arbitraria.

También hay muchos algoritmos especializados para encontrar raíces de polinomios en el artículo de Wikipedia .

Answer 2

Editar - Grandes noticias: A continuación doy un algoritmo y muestro que funciona, en un sentido específico. Reconozco que no funciona perfectamente, en cierto sentido. ¡Tenemos una solución para eso - un poco de pre-cálculo y funciona perfectamente!

El algoritmo siguiente produce $a_1,\dots,a_k$ de manera que cada $a_j$ está dentro de $\epsilon$ de un cero de $p$ y cada cero de $p$ está dentro de $\epsilon$ de algunos $a_j$ . No es perfecto, el problema es que puede ocurrir que varios $a_j$ se agrupan en torno al mismo cero; en cierto sentido, hemos localizado todos los ceros dentro de $\epsilon$ pero no lo sabemos que $a_j$ están cerca de más de un cero. Si sabíamos que $|z_i-z_j|\ge 5\epsilon$ por cada dos ceros $z_i\ne z_j$ podríamos arreglar eso. (Porque entonces $D(a_j,\epsilon)$ contiene exactamente un cero $z_j$ de $p$ y si $|a_i-a_j|<2\epsilon$ entonces $z_i=z_j$ . Así que $(a_1,\dots,a_k)$ nos da $(z_1,\dots,z_n)$ : Sabemos que hay un cero $z_1\in D(a_1,\epsilon)$ . Llame a $a_1$ la aproximación a $z_1$ . Tire el $a_i$ con $|a_i-a_1|<2\epsilon$ y continuar...)

Así que podemos arreglar eso:

Lema 1. Dado $p$ sin raíces múltiples podemos calcular $\epsilon>0$ tal que dos raíces cualesquiera de $p$ están separados por al menos $\epsilon$ .

En primer lugar

Lema 2. Dado $p$ sin raíces múltiples podemos calcular $\delta>0$ tal que $p(z)=0$ implica $|p'(z)|\ge\delta$ .

Prueba: El algoritmo euclidiano nos da $r_1,r_2\in\Bbb Q[z]$ tal que $$1=r_1(z)p(z)+r_2(z)p'(z).$$ Calcular $R$ así que $p(z)=0$ implica $|z|\le R$ y, a continuación, calcular $\delta$ para que $|z|\le R$ implica $|r_2(z)|\le1/\delta$ .

Prueba del lema 2: Digamos que todos los ceros de $p$ están contenidas en $|z|\le R$ y decir $|p''(z)|\le c$ en $|z|\le R$ . Sea $\delta$ sea como en el lema 1, y supongamos que $p(a)=0$ . Entonces el teorema de Taylor muestra que $$|p(z)|\ge\delta|z-a|-\frac12 c|z-a|^2\quad(|z|\le R),$$ por lo que si $0<|z-a|<2\delta/c$ entonces $|p(z)|>0$ .

Así que se comienza el algoritmo sustituyendo primero $p$ por $p/gcd(p,p')$ para asegurarse de que no hay ceros múltiples y luego disminuir $\epsilon$ si es necesario para asegurar que todos los ceros están separados por $3\epsilon$ .

En este momento este post empieza más o menos en la mitad - quería poner la gran noticia en la parte superior (el Lemma 1 es algo que me ha estado molestando). Comienza aquí:

Comienza: De hecho, el algoritmo en la respuesta aceptada hace funcionan, al menos con un poco de elaboración en cuanto a cómo decidimos rechazar o aceptar un determinado cuadrado durante la subdivisión.

Editar: Por alguna razón pensé que la respuesta aceptada utilizaba una búsqueda binaria, como la que se muestra a continuación. Mirando de nuevo veo que no lo hace, sino que comienza con una cuadrícula de cuadrados de tamaño $\epsilon$ . Así que no importa lo que he dicho sobre la respuesta aceptada. Parece probable que la búsqueda binaria a continuación sea mucho más eficiente, ya que un cuadrado de tamaño mucho mayor que $\epsilon$ a menudo será rechazado en un solo paso. También hay una prueba más abajo de que la búsqueda binaria funciona.

En cuanto a eficiencia : Hice una pequeña implementación en Python (ver aplicación abajo) y buscó las raíces de $z^2+1$ con varios $\epsilon$ . Los resultados aparecieron "instantáneamente". No me molesté en averiguar cómo cronometrarlo realmente, pero conté el número de iteraciones, es decir, el número de cuadrados "probados". Obtuve esto:

eps = 0.1: count = 516
eps = 0.01: count = 868
eps = 0.001: count = 1156
eps = 0.0001: count = 1428
eps = 1e-05: count = 1812
eps = 1e-06: count = 2100
eps = 1e-07: count = 2388
eps = 1e-08: count = 2772
eps = 1e-09: count = 3060
eps = 1e-10: count = 3348

Nota que es obsoleta. Con la versión mejorada de reject( $Q$ ) a continuación obtenemos count = 1108 para $\epsilon = 10^{-10}$ .

La dependencia del tiempo de ejecución de $\epsilon$ aparece mucho mejor de lo que esperaba - ciertamente aniquila totalmente el $c/\epsilon^2$ en la respuesta aceptada.

En retrospectiva, esto tiene sentido. Digamos que ejecutamos el algoritmo con $\epsilon=\epsilon_0$ y luego otra vez con $\epsilon=\epsilon_0/10$ . La definición de reject( $Q$ ) es independiente de $\epsilon$ Por lo tanto, cada cuadrado rechazado la primera vez también será rechazado la segunda vez; el único trabajo extra en la segunda ejecución es la subdivisión de los cuadrados que fueron aceptados la primera vez, y el número de cuadrados aceptados es muy pequeño. (Con $z^2+1$ había ocho plazas aceptadas para cada $\epsilon$ que es lo que esperaba, o al menos esperaba). Hmm, según ese argumento debería ser algo así como $\log(1/\epsilon)$ lo que parece ser más o menos consistente con el experimento anterior.

Funciona" en este sentido : Dado $\epsilon>0$ podemos encontrar un número finito de puntos $a_1,\dots a_k$ de manera que cada $a_j$ está dentro de $\epsilon$ de un cero y cada cero está dentro de $\epsilon$ de algunos $a_j$ .

Es decir, podemos aproximar el conjunto cero en la métrica de Hausdorff. Eso parece que debería ser suficiente. No funciona a la perfección, en el sentido de que podría ocurrir que dos ceros estén cerca del mismo $a_j$ . Está claro que las raíces múltiples van a causar problemas con este tipo de búsqueda; parece claro que dos raíces cercanas también van a dificultar las cosas. Podemos eliminar las raíces múltiples, pero no hay mucho que podamos hacer con dos raíces cercanas...

Se dijo que el algoritmo aceptado no funciona porque $|p(z)|$ pequeño no implica que $z$ está cerca de un cero de $p$ . Pero hace implican que bajo algunas condiciones, para un valor adecuado de "pequeño":

Lema 3. Supongamos que $p$ está completo, $p(a)\ne0$ y $p'(a)\ne 0$ . Sea $r=2|p(a)|/|p'(a)|$ y elija $\delta_2>0$ tal que $|p''(z)|\le \delta_2$ para $|z-a|\le r$ . Si $\frac12|p(a)|\delta_2<|p'(a)|^2$ entonces $p$ tiene exactamente un cero en $|z-a|<r$ .

Prueba de ello: Sea $g(z)=p(a)+p'(a)(z-a)$ y observe que $g$ tiene exactamente un cero en $|z-a|<r$ . Si $|z-a|=r$ entonces $|g(z)|\ge r|p'(a)|-|p(a)|=|p(a)|$ mientras que el teorema de Taylor muestra que $|p(z)-g(z)|\le\frac12 r^2\delta_2$ . La hipótesis implica que $\frac12 r^2\delta_2<|p(a)|$ por lo que el teorema de Rouche dice que $p$ tiene el mismo número de ceros que $g$ en $|z-a|<r$ .

Ahora, el algoritmo: Dejemos que $\epsilon>0$ . Sustitución de $p$ por $p/gcd(p,p')$ podemos suponer que cada cero de $p$ es simple.

Dados los coeficientes de $p$ ciertamente podemos encontrar un cuadrado (cerrado) $Q_0$ que contiene todos los ceros de $p$ y entonces podemos encontrar $\delta_1$ y $\delta_2$ tal que $|p'(z)|\le\delta_1$ y $|p''(z)|\le\delta_2$ siempre que $d(z,Q_0)\le\epsilon$ .

Construimos un árbol de cuadrados, dividiendo $Q_0$ en cuatro subcuadros iguales y luego dividir cada uno de ellos en cuatro subcuadros, etc. Para cada cuadrado $Q$ obtenida de esta manera realizamos dos pruebas:

def rechazar( $Q$ ):

Diga $a$ es el centro de $Q$ y que $r$ sea la distancia de $a$ a una esquina de $Q$ . Si $|p(a)|>r\delta_1$ declaramos $Q$ sea un nodo terminal y no realice ninguna otra subdivisión; a continuación diremos que $Q$ fue "rechazado". Tenga en cuenta que si $Q$ es rechazado está claro que $Q$ no puede contener un cero de $p$ .

Editar: Eso es estúpido, usar un límite superior global en $|p'|$ para un límite superior en $Q$ . De hecho $\min(\delta_1,|p'(a)|+r\delta_2)$ es un límite superior para $|p'|$ en $Q$ . Así que cambiar lo anterior: Rechazar $Q$ si $|p(a)|>r\min(\delta_1,|p'(a)|+r\delta_2)$ . (Esto hizo un enorme mejora en la búsqueda de los ceros de $\sin(z)$ en $|z|<4$ ya que $\cos$ suele ser mucho menor que su máximo global).

De hecho, llegamos a las plazas pequeñas muy rápidamente, y una vez $Q$ es pequeño, rara vez o nunca ocurrirá que $\min(\delta_1,|p'(a)|+r\delta_2)=\delta_1$ . Entonces, ¿por qué no olvidar $\delta_1$ y simplificar las cosas diciendo que rechazamos $Q$ si $|p(a)|>r(|p'(a)|+r\delta_2)$ ? Eso podría funcionar, pero estropea la prueba de abajo de que el algoritmo termina.

def aceptar( $Q$ ):

Diga $a$ es el centro de $Q$ y que $r=2|p(a)|/|p'(a)|\in[0,\infty]$ . Si $diam(Q)<\epsilon$ , $\frac12|p(a)|\delta_2<|p'(a)|^2$ y $r<\epsilon$ declaramos $Q$ para ser un nodo terminal y no realizar ninguna otra subdivisión; diremos a continuación que $Q$ fue "aceptado". Tenga en cuenta que si $Q$ se acepta entonces $p$ tiene exactamente un cero en $|z-a|<r$ por el lema. (A menos que $p(a)=0$ , en cuyo caso $r=0$ ; en cualquier caso $p$ tiene un cero en $|z-a|\le r$ .)

Nota puede y probablemente sucederá que rechace( $Q$ )=aceptar( $Q$ )=Verdadero. No me di cuenta hasta que ejecuté algo de código que esto significa que es importante probar reject( $Q$ ) primero. Hacerlo al revés aumenta la longitud de la lista $a_1,\dots,a_k$ volvió, con más $a_j$ agrupación alrededor de cada cero.

Si el algoritmo termina, hemos dividido $Q_0$ en un número finito de cuadrados, cada uno de los cuales fue aceptado o rechazado, y hemos terminado.

Por qué hemos terminado , suponiendo que el algoritmo termine: Digamos que $a_1,\dots,a_k$ son los centros de los cuadrados aceptados. Para cada $j$ , $p$ tiene exactamente una raíz en $|z-a_j|\le r_j$ ; ya que $r_j<\epsilon$ esto dice cada $a_j$ está dentro de $\epsilon$ de una raíz. A la inversa, como los cuadrados rechazados no contienen ninguna raíz, toda raíz se encuentra en algún cuadrado aceptado $Q$ ; ya que $diam(Q)<\epsilon$ esto demuestra que cada raíz está dentro de $\epsilon$ de algunos $a_j$ .

El algoritmo termina.

Prueba: Por el lema de Konig basta con demostrar que nuestro árbol no tiene infinitas ramas. Digamos que $Q_1,Q_2,\dots$ es una secuencia de cuadrados tal que $Q_{k+1}$ es uno de los cuatro cuartos de $Q_k$ necesitamos demostrar que algunos $Q_k$ se acepta o se rechaza.

Diga $a_k$ es el centro de $Q_k$ . Diga $a_k\to z$ .

Supongamos que $p(z)\ne0$ y que $r_k$ sea la distancia desde $a_k$ a una esquina de $Q_k$ . Entonces ciertamente existe $k$ tal que $|p(a_k)|>r_k\delta_1$ ya que $r_k\to0$ pero $p(a_k)\not\to0$ . Así que $Q_k$ es rechazado.

Supongamos que $p(z)=0$ . Entonces $p'(z)\ne0$ . Así que $p(a_k)\to0$ y $p'(a_k)\not\to0$ y si se observan las condiciones de la definición de "aceptado", esto deja claro que algunos $Q_k$ fue aceptada.

La puesta en práctica: También puedo incluir el código que mencioné antes, por si alguien quiere probarlo. Nota que esto no es exactamente lo que se describe arriba: search(Q_0,f,df,eps,d1,d2) funciona para cualquier función holmórfica $f$ , devolviendo todos los ceros en $Q_0$ . En particular, no trata de calcular un $Q_0$ que contiene todos los ceros, ya que $f$ no es un polinomio; tampoco calcula $\delta_1$ y $\delta_2$ hay que calcular los límites superiores adecuados d1 y d2 para $|f'|$ y $|f''|$ a ti mismo. También simplemente no es correcto, ya que utiliza cálculos de punto flotante para comprobar si $p(a)=0$ . No pretende ser una versión de producción, sólo una prueba rápida y sucia de si la cosa realmente funciona. Lo hace, mucho mejor de lo que esperaba... (Si obtienes un SyntaxError o IndentError he metido la pata añadiendo cuatro espacios al inicio de cada línea).

En fin:

from math import sqrt

class square():
  def __init__(Q,  a, r):
    """If a= x+iy and r>0 then square(a,r) is supposed to represent the
square [x-r,x+r] x [y-r,y+r]."""
    Q.a = a
    Q.r = float(r)

  def quarters(Q): 
    vs = [0.5+0.5j, 0.5-0.5j, -0.5+0.5j, -0.5-0.5j]
    return [square(Q.a + v*Q.r, Q.r/2)for v in vs]

def search(Q_0,f,df,eps,d1,d2):
  """f should be analytic  in the eps-nbd of Q_0. NOTE proof that the algorithm
terminates requires that all zeroes of f are simple. df=f', d1 d2 upper bounds
for |f'| and |f''| in eps-nbd of Q_0. Returns list approximating all zeroes  in
Q_0, in the sense that everything in returned list is within eps of a zero
and every zero is within eps of something in list"""
  s = sqrt(2)

  def reject(Q):
    #return abs(f(Q.a)) > s*d1*Q.r
    #No! There's a better upper bound for the derivative in Q
     #if the derivative at a happens to be small:
    return abs(f(Q.a)) > s*min(d1,abs(df(Q.a))+d2*s*Q.r)*Q.r

  def accept(Q):
    z = abs(f(Q.a))
    dz = abs(df(Q.a))
    return (2*s*Q.r<eps) and (2*z<eps*abs(dz)) and (z*d2<2*dz*dz)

  todo = [Q_0]
  accepted = []

  while len(todo)>0:
    newtodo = []
    for Q in todo:
      for q in Q.quarters():
        if not reject(q):
          if accept(q):
            accepted.append(q)
          else:
            newtodo.append(q)
    todo = newtodo
  return [Q.a for Q in accepted]

Answer 3

El artículo de la wikipedia http://en.wikipedia.org/wiki/Root-finding_algorithm ofrece enlaces a muchos métodos diferentes para encontrar raíces de polinomios. (Comience en la sección titulada "Encontrar raíces de polinomios".) Muchos de los métodos son incomparables, en el sentido de que funcionan más rápido o más lento que otros dependiendo del polinomio específico.

Answer 4

El método de continuación de homotopía es bueno para encontrar todas las soluciones COMPLEJAS con una precisión arbitraria, y está implementado en el paquete de Geometría Algebraica Numérica en Macaulay 2, por ejemplo. El método es más general. Puede resolver un sistema de ecuaciones polinómicas en muchas variables. De hecho, es un problema más difícil encontrar todas las soluciones REALES SIN encontrar todas las soluciones complejas.

Por lo que tengo entendido, la región de la solución no necesita estar acotada para que la continuación de homotopía funcione. También puedes "proyectar" tu problema si es necesario, para no tener que preocuparte de que los caminos de homotopía se vayan al infinito. Algunos métodos asumen que las soluciones son todas simples, pero hay formas de evitarlo. Una de ellas es el método de "deflación".

Answer 5

Uno de los semi-recomendados para encontrar raíces en el plano complejo es el método de Laguerre, que por alguna razón no está incluido en el artículo de Wikipedia sobre búsqueda de raíces.

http://en.wikipedia.org/wiki/Laguerre's_method

La razón por la que sé esto es una conferencia en un coloquio, hace mucho tiempo, de Steven Smale sobre la complejidad del método de Newton, durante la cual William Kahan se levantó y sostuvo por qué el método de Newton no valía nada y el de Laguerre era mucho mejor.

No puedo decir si insiste en encontrar todas las raíces con gran precisión. Tal vez se podría dividir por $(x - r_k)^{n_k}$ cada vez que una raíz $r_k$ con multiplicidad $n_k$ y buscar las raíces del nuevo polinomio, utilizando esos resultados como valores semilla para encontrar las raíces exactas utilizando el polinomio original.