Mandelbrot en Mathematica: afinando el código.

Es difícil encontrar una imagen de

Benoît Mandelbrot en la que no

parezca el Maestro Joda.

En esta entrada comenzamos a hablar del conjunto de Mandelbrot, describimos en qué consiste y vimos una forma bastante burda de representarlo en Mathematica, si bien dicho código tiene la ventaja de ser muy sencillo, así como ser válido para "vislumbrar" conjuntos de Julia.







Una de las principales pegas que tiene esa forma de representar el conjunto de Mandelbrot es que no tiene condición de parada ninguna, por lo que para cada valor del plano complejo se realiza la totalidad de las iteraciones fijadas de antemano. Para solucionar ese problema utilizaremos que está demostrado que al aplicarle el proceso iterativo
$$\begin{matrix}

{z}_{0}=0 \\
{z}_{n+1}={{z}_{n}}^2+c \\
\end{matrix}$$
a un cierto valor de c, si obtenemos un valor con módulo superior a 2, entonces la divergencia está asegurada, por lo que sabemos que dicho valor no pertenece al conjunto de Mandelbrot sin necesidad de seguir iterando.


Vamos a utilizar un código que a diferencia del anterior es muy descriptivo y aunque es largo muestra muy bien el proceso de construcción del conjunto; es en un 80% obra de Novak Petrovic.


Comenzamos. Los comentarios los pongo entre signos de (* *) de forma que se puede copiar y pegar del tirón todo en un notebook de Mathematica sin que los comentarios impidan la ejecución del código.

En primer, lugar puesto que vamos a comprobar qué valores de c pertenecen al conjunto de Mandelbrot de entre los contenidos en el rectángulo complejo cuya diagonal está definida por -2-2i y 2+2i, definimos así los esquinas inferior izquierda y superior derecha del rectángulo.

cpaso determinará el número de puntos para el que comprobaremos si pertenece al conjunto; marcará la diferencia entre cada parte real o imaginaria.
rmax definirá el valor límite del módulo para seguir iterando.
Nmax será el número máximo de iteraciones a realizar en cualquier caso.


esqinfizq = -2.0 - 2.0 *I;
esqsupder = 2.0 + 2.0 *I;
cpaso = 0.05;
rmax = 2.0;
Nmax = 50;

(* definimos ahora la sucesión iterativa que aplicaremos a cada valor del rectángulo complejo*)

znmas1[zn_, c_] := zn^2 + c;

(* En resultados almacenaremos los valores a dibujar, para después usar el comando DensityPlot. Necesitaremos también una variable auxiliar que llamaremos resultadosaux *)

resultados = {};

(* Vamos ahora tomando valores para la parte imaginaria de c, que llamaremos y, de entre los contenidos en nuestro rectángulo definido arriba por dos de sus esquinas*)

For[cpartey = Im[esqinfizq], cpartey <= Im[esqsupder],  cpartey += cpaso,

(* Ahora en resultadosaux, guardaremos todos los valores correspondientes a las iteraciones para cada valor de c con una cierta parte imaginaria y, lo que supondrá una fila para valor de y compuesta por todas las posibles partes reales con las que combinarse, que llamaremos x*)

resultadosaux = {};

(* Para cada valor de y vamos tomando todos los posibles de x en nuestra rectángulo, avanzado al paso definido al principio por la variable cpaso *)

For[cpartex = Re[esqinfizq], cpartex <= Re[esqsupder], cpartex += cpaso,

(* Bien, en este punto tenemos un cierto valor de c, con su parte real x y su imaginaria y; lo llamaremos c y será el valor de la primera iteración. Ponemos el contador de las iteraciones para este valor de c a 0 *)

contador = 0;
c = cpartex + I*cpartey;
zvalor = c;

(*La variable control la usaremos para interrumpir las iteraciones si se cumple una cierta condición*)

control = False;
While[!control,

(* Comprobamos si el módulo del valor obtenido en una iteración es superior al definido al principio, o bien si hemos superado el número máximo de iteraciones*)

If[Abs[zvalor] > rmax || contador > Nmax,

(* En ese caso, la variable contador tomará el valor cierto, habremos acabado con ese valor de c y añadimos el número de iteraciones empleado a la fila de valores de x correspondiente de resultadosaux*)

control = True;
AppendTo[resultadosaux, contador],

(*En caso contrario, pasamos a la siguiente iteración e incrementamos en 1 el contador de iteraciones*)

zvalor = znmas1[zvalor, c];
      contador++;
   
      ];
    ];
  ];

(* Una vez hallamos recorrido todos los valores de x del rango correspondientes a un cierto valor de y, añadimos la fila de valores obtenidos tras las iteraiones a la variable resultados, y pasaremos al siguiente valor de y*)

AppendTo [resultados, resultadosaux];

  ];

(* Por último dibujamos los valores contenidos en resultados, que serán las iteraciones llevadas a cabo para cada valor del plano complejo definido de antemano, tomándo más o menos valores según el cambio de valor indicado por la variable cpaso. Usaremos la función de color Hue que queda bastante bien, eso es lo de menos*)

ListDensityPlot[resultados, ColorFunction -> Hue]

Obteniendo la siguiente representación del conjunto:

Un código laborioso ahora, pero a diferencia del de la entrada anterior describe muy bien lo que estamos haciendo para obtener la representación del conjunto de Mandelbrot. El resultado es mejor por otra parte: si un punto pertenece al conjunto, pertenece y ya está, no se muestra de un color distinto según el valor del módulo obtenido tras las iteraciones, y si un punto no pertenece al conjunto se muestra de un cierto color según las iteraciones que hayan sido necesarias para obtener un módulo superior a 2.

En el término medio está la virtud dicen, y eso es lo que haremos en la tercera y última entrada dedicada al tratamiento del conjunto de Mandelbrot con Mathematica.