Encontrar todas las raíces de un polinomio usando el método de Newton-Raphson.

Question

¿Existe una formulación general para encontrar todas las raíces de un polinomio, especialmente las complejas, utilizando el Método de Newton-Raphson?

Answer 1

Sí, ¡existe tal método! Consulta el llamado "Cómo encontrar todas las raíces de polinomios complejos por el método de Newton", de Hubbard, Schliecher y Sutherland.

No solo hay un método sino que, debido a la estabilidad de los puntos fijos bajo la iteración de la función del método de Newton, hay un método muy bueno. De hecho, los autores señalan que fueron llevados a este tema en su estudio de medidas invariantes del mapeo de Henon y, en ese contexto, tuvieron una necesidad genuina de encontrar raíces de polinomios cuyo grado era de varios miles. Evidentemente, su método tuvo éxito donde las herramientas de software establecidas fallaron.

La idea básica es encontrar una colección de semillas iniciales distribuidas de tal manera que se garantice que, para cada raíz, hay al menos una de las semillas que converge a esa raíz. Este conjunto es bastante grande pero puedes detenerte cuando hayas encontrado todas las raíces. La multiplicidad de la raíz se puede determinar por la velocidad de convergencia.

He incluido una implementación en Mathematica a continuación llamada NewtonSolve. Aquí tienes un ejemplo:

NewtonSolve[(z^3 - 1) (z - I)^3, z, 
  ErrorTolerance -> 10^-20] // Chop
(* Salida: {{2.34482*10^-8 + 1. I, 1}, {0. + 1. I, 3}, 
   {-0.5 + 0.866025 I, 1}, {-0.5 - 0.866025 I, 1}}
*)

Nota que el comando devuelve una lista de pares {raíz, multiplicidad} y que se identificó correctamente que $i$ tiene una multiplicidad de 3. Por supuesto, hay polinomios que desafiarán a cualquier buscador de raíces y no afirmo que este código sea de calidad de producción, pero sí ilustra las ideas.

Aquí tienes una imagen que ilustra aún más lo que está sucediendo con este ejemplo:

Los puntos verdes grandes son exactamente las raíces del polinomio, cada una dentro de su cuenca de atracción sombreada. La observación fundamental en el documento es que cada una de estas cuencas de atracción se extiende hasta el infinito. Como resultado, podemos elegir un anillo suficientemente grande de semillas iniciales para que cada cuenca de atracción contenga una de las semillas. El anillo para este polinomio en particular (de hecho, cualquier polinomio de grado seis con raíces en el disco unitario) se muestra como el anillo de puntos negros.

La función resuelve un ejemplo supuestamente difícil propuesto en los comentarios anteriores, con la multiplicidad correcta, en menos de un segundo:

NewtonSolve[(z - 1)^50, z] // Timing
(* Salida: {0.293804, {{1. + 0. I, 50}}} *)

---

Aquí está el código para NewtonSolve:

Options[NewtonSolve] = {
   NormalizationFactor -> Automatic,
   ErrorTolerance -> $MachineEpsilon
   };
NewtonSolve[poly_, var_,
   opts___] := Module[
    {p, z, nwt, cp, cpp,
     coeffs, deg, maxRootRadius,
     mult, sameRootQ, newRootQ,
     maxIters, rootsSoFar, initPoints,
     numInitPoints, rmPair, normFactor,
     tol},

    normFactor = NormalizationFactor /. {opts} /. 
      Options[NewtonSolve];
    tol = ErrorTolerance /. {opts} /.
      Options[NewtonSolve];

    coeffs = CoefficientList[poly, var];
    deg = Length[coeffs] - 1;
    If[normFactor === Automatic,
     maxRootRadius = Max[Abs[Drop[coeffs, -1]/Last[coeffs]]] + 1,
     maxRootRadius = normFactor];
    p[z_] = poly /. var -> maxRootRadius*z;

    cp = Compile[{{z, _Complex}}, p[z]];
    cpp = Compile[{{z, _Complex}}, Evaluate[p'[z]]];
    nwt = Compile[{{z, _Complex}},
      Evaluate[z - p[z]/p'[z]]];
    mult[z0_Complex] := If[nwt[z0] == z0, 1,
      mult[FixedPointList[nwt, z0, 1000]]];
    mult[orbit_List] := Module[{f, div, diffs},
      f[x_] := If[x != 1.0, 1/(x - 1)];
      div[{a_, b_}] := If[b != 0, a/b];
      diffs = Abs[(#[[2]] - #[[1]]) & /@ 
         Partition[Drop[orbit, -1], 2, 1]];
      statisticalMode[Round[f /@ (div /@ 
             Partition[diffs, 2, 1])] + 1  /. Round[_] -> 0]];
    sameRootQ[z1_, z2_, 1] := If[Abs[cp[z1]] >= Abs[cp[z2]],
      Abs[z1 - z2] <= Abs[3 cp[z1]/cpp[z1]],
      Abs[z2 - z1] <= Abs[3 cp[z2]/cpp[z2]]];
    sameRootQ[z1_, z2_, mm_ /; mm > 1] := Abs[z1 - z2] < 100 tol;
    newRootQ[z_, mm_] := Module[{justRootsSoFar, flag},
      flag = True;
      justRootsSoFar = First /@ Flatten[rootsSoFar];
      Scan[If[sameRootQ[#, z, mm], Return[flag = False]] &,
       justRootsSoFar];
      flag];
    maxIters = 
     Ceiling[(Log[1 + Sqrt[2]] - Log[tol])/(
      Log[deg] - Log[deg - 1])];
    rootsSoFar = {};
    numRootsSoFar = 0;
    initPoints = N[S[deg]];
    numInitPoints = Length[initPoints];

    i = 0;
    While[i++ < numInitPoints  && numRootsSoFar < deg,
     orbit = NestWhileList[nwt, initPoints[[i]],
       (Abs[cp[#1]] >= tol && #1 != #2) &, 2, maxIters];
     nextRoot = Last[orbit];
     If[(Abs[cp[nextRoot]] < tol || orbit[[-2]] == orbit[[-1]]),
      m = mult[nextRoot];
      If[m > 1,
       nextRoot = FixedPoint[nwt, nextRoot, maxIters]];
      If[newRootQ[nextRoot, m],
       numRootsSoFar = numRootsSoFar + m;  
       rootsSoFar = {rootsSoFar, rmPair[nextRoot, m]}]]
     ];
    Flatten[rootsSoFar] /. rmPair[r_, mm_] -> 
      {r*maxRootRadius, mm}
    ] /; (PolynomialQ[poly, var] && 
     NumericQ[poly /. var :> Random[]]);

statisticalMode[list_] :=
      Module[{c, mc, ms},
            c = {Length[#], First[#]} & /@ Split[Sort[list]];
            If[Length[c] === 1, Return[c[[1, 2]]]];
            mc = Max[First[Transpose[c]]];
            If[mc === 1, Return[{}]];
            ms = Cases[c, {mc , val_} -> val];
            If[Length[ms] == 1, ms[[1]], ms];
    ms[[1]]] /; VectorQ[list] && (Length[list] > 0);
statisticalMode[{}] = 1;

ring[r_, d_] := With[{n = Ceiling[8.32547 d Log[d]]},
   Flatten[Transpose[Partition[Table[r Exp[2 Pi I j/n], 
       {j, 0, n - 1}], Ceiling[n/d], Ceiling[n/d], 1, {{}}]]]];
S[d_] := With[{R = 1 + Sqrt[2], s = Ceiling[0.26632 Log[d]]},
   Flatten[Table[ring[R (1 - 1/d)^((2 \[Nu] - 1)/(4 s)), d],
     {\[Nu], 1, s}]]];

Answer 2

Tal como se ha señalado, el método de Newton-Raphson, al igual que la mayoría de otros métodos similares, convergen a una sola raíz aislada. Sin embargo, existen métodos para encontrar todas las raíces de un polinomio.

Uno que conozco, es crear la matriz compañera, y luego encontrar (numéricamente) los eigenvalores de esa matriz, que son precisamente las raíces de tu polinomio.

Si utilizas Matlab, esto es precisamente lo que hace la rutina 'roots'. Por experiencia personal con esa rutina, corre muy rápido tanto para polinomios de grados bajos como altos.

Answer 3

Si bien en general la respuesta es no, existe el método de Weierstrass, el cual se basa en el método de Newton multidimensional e intenta encontrar todas las raíces complejas. Al igual que el método de Newton, no se garantiza que tenga éxito y solo funciona bien cuando las raíces son simples.

El método de Weierstrass probablemente es lo más cercano a lo que estás buscando.