Fractal de Newton

Fractal de Newton-Raphson Por Ariel S.

Para crear este fractal se realiza un proceso similar al de Julia. El color de cada punto del fractal dependerá de la velocidad con la que dicho punto converge a una de las raices de un polinomio dado, utilizando el método numérico de aproximación de racies de Newton-Raphson.

Si tomamos la función tendremos que dependiendo de la condición inicial $z_{0}=x_{0}+iy_{0}$ que tomemos al aplicar el método de Newton, en caso de converger lo hará a alguna de las raíces de Recordemos que éste es un método iterativo que funciona de la siguiente manera a partir del $z_{0}$ dado generando la sucesión definida por

MATH

Lo que hacemos entonces es tomar un punto fijo $z_{0}$ $\in \U{2102}$ y generar la sucesión descripta. Dado un cierto grado de tolerancia para el error de la aproximacion lo que nos importará a la hora de graficar es la cantidad de iteraciones que se tuvieron que realizar para que se haya convergido a alguna raiz, siempre y cuando no se nos haya anulado la suceción para algún valor (para evitar dividir por cero).

Por ejemplo si tomamos la función $f(z)=z^{3}-1$ y un determinado punto $z_{0}$ tendremos definida la suceción como MATH

Así podemos definir la siguiente función que tiene como parámetros de entrada al punto que se toma como condición inicial para el método de Newton, y un valor de escape del bucle itera_max que será el número máximo de términos de la suceción a calcular. Y retornará el valor correspondiente al número de iteraciones necesarias para que haya convergencia a una raíz de en caso de no haber superado el número máximo de iteraciones previsto, y que el valor de ningún término de la suceción se anule.

int Newton(complejo z,int itera_max)
{ complejo w;
    int k=0;
    w=z;
    double m=Cmod( CsumaR(Cnpot(w,3),1) );
    while ((m>0.001) && (k<=itera_max) && (w.x!=0) && (w.y!=0) )
    {      w=Cdiv( CsumaR(CprodR(Cnpot(w,3),2),-1), CprodR(Cnpot(w,3),2) );
            m=Cmod(w);
            k++;
     }
    return(k);
}

donde Cmod es la función que retorna el módulo o valor absoluto de un número complejo

double Cmod(complejo a)
{ double m=a.x*a.x+a.y*a.y;
return(m);
}

CsumaR es la función que devuelve la suma de un número complejo y uno real

complejo CsumaR(complejo a, double r)
{    complejo c;
      c.x=a.x+r;
      c.y=a.y;
      return(c);
}

CprodR es la función que retorna el producto de un número complejo y uno real

complejo CprodR(complejo a, double r)
{     complejo c;
      c.x=r*a.x;
      c.y=r*a.y;
      return(c);
}

y Cnpot devuelve la enésima potencia natura de un complejo dado

complejo Cnpot(complejo a, int n)
{    complejo c,u;
      bool nNegativo;
      c=a;
      if (n<0) {nNegativo=true;n=-n;};
      if (n==0) {c.x=1; c.y=0;}
      else if (n==1) {c=a;}
      else
      {       for(int k=1;k<n;k++)
               { c=Cprod(a,c);
                }
      }
      if (nNegativo==true){
            u.x=1;
            u.y=0;
            Cdiv(u,c);
      }
      return(c);
}

Algunas optimizaciones

Puede que al graficar este fractal con la función recién descripta se note una cierta lentitud, esto es porque la cantidad de cálculos que se realizan es alto.Podemos implementar algunas optimizaciones para ganar un poco de velocidad al realizar esas mismas tareas, valiéndose de las propiedades de los números complejos:

MATH

escribiendo

MATH

Entonces ahora podemos reescribir la función de la siguiente manera

int NewtonOpt(complejo z,int itera_max)
{    complejo w;
      int k=0;
      w=z;
      double tempx,d,m=Cmod( CsumaR(Cnpot(w,3),1) );
      while ((m>0.001) && (k<=itera_max) && (w.x!=0) && (w.y!=0) )
      {        d=3*((w.x*w.x-w.y*w.y)*(w.x*w.x-w.y*w.y)+4*w.x*w.x*w.y*w.y);
               tempx=w.x;
               w.x=2*w.x/3 + (w.x*w.x-w.y*w.y)/d;
               w.y=2*w.y/3-2*tempx*w.y/d;
               m=Cmod(w);
               k++;
      }
      return(k);
}

La imagen que obtenemos al programar esto es similar a la siguiente

Dejá un comentario, sugerencia u opinión