Mandelbrot en Mathematica: tres líneas de código

Sí, hay a quien le gusta...

Entre las imágenes más sugerentes en matemáticas están sin duda los fractales. Se tratarán con frecuencia en este blog, y hoy vamos a empezar a ver cómo se puede dibujar con Mathematica el más famoso de todos ellos: el conjunto de Mandelbrot.







Existen varios sitios donde puedes leer excelentes explicaciones de este conjunto, como esta entrada de gaussianos o Mu-Ency, la enciclopedia del conjunto de Mandelbrot, donde podemos encontrar cosas tan curiosas como este conjunto de Mandelbrot hecho con texto:

vL , '*m-` -m/**\a, ... _,/#, ]),., ., '#F-.F~*^' '`'*~*eae/: . -__/* '`_* )_. ic,_ ./- T\a 7F*-~~*a, /` dL \_,\F^ '\` .*` .,___\____._/*^* R2a .m~ ` ' :*r(, . ~e ` -\r._____/`L @e _-**^)~*' ]c .^_ -~ ` .)L_, _-`' .]F_.,, ,.-e_@, -*#^ *~*v- '-/-~* '`]- d@,,..(, ^'*m-' i*^^` '^`



Sin  embargo no me parece bien empezar a hablar de cómo se dibuja una cosa sin explicar al menos de qué va dicha cosa, así que vamos a intentar contarlo por encima.

Consideremos la siguiente sucesión iterativa: ${z}_{n+1}={{z}_{n}}^2+c$ en el plano complejo.

Para cada valor de c fijo de antemano tenemos una sucesión iterativa distinta. Si se la aplicamos a varios valores complejos distintos, es decir, arrancamos la secuencia en valores distintos, obtendremos que para algunos de esos valores la secuencia de valores que se obtiene converge a un valor concreto, mientras que para otros diverge a infinito. La frontera en el plano complejo entre los valores para los que converge la secuencia y los valores para los que diverge se llama conjunto de Julia para el valor de c tomado de antemano.

c=i

c=-0.5+0.7i

Aquí podemos ver dos conjuntos de Julia para valores de c diferentes, y como se puede ver uno es conexo (de una sola pieza digamos) mientras que el otro no .

Pues bien, el conjunto de Mandelbrot se puede definir como el conjunto de los valores de c para los que el conjunto de Julia asociado es conexo. Afortunadamente no hay que probar a construir cada conjunto de Julia y ver qué pasa, ya que está demostrado que cada conjunto de Julia será conexo si para el valor 0 la sucesión formada por las iteraciones es convergente, es decir, será conexo para un cierto valor de c si la sucesión

$$\begin{matrix}
{z}_{0}=0 \\
{z}_{n+1}={{z}_{n}}^2+c \\
\end{matrix}$$

es convergente. Así que los valores complejos de c para los que no diverge a infinito en la sucesión anterior formaran el conjunto de Mandelbrot.

Una vez explicado en qué consisten los conjuntos de Julia y el de Mandelbrot, vamos a ver cómo podemos dibujarlos en Mathematica hasta de tres maneras distintas. En este post, para que no quede muy largo vamos a ver la primera de ellas, que no es más que una aproximación al conjunto muy bestia, pero con un código muy sencillo.

Empezamos con el código. Estas dos primeras líneas dejan fijo el número de iteraciones que realizaremos (así, a lo basto, sin condición de parada) y la sucesión iterativa, que se aplicará recurrentemente llamándose a sí misma: la iteración 15 llamará a la 14, ésta a la 13, y así hasta la inicial $z_0$ que vale 0.

Niter = 15;
f[c_, it_] := If[it > 0, f[c, it - 1]^2 + c, 0];

Y sólo necesitamos una tercera línea para dibujar los puntos del rectángulo complejo formado por los números con partes real e imaginaria entre -2 y 2 con distinto color según el módulo que hayan alcanzado después de las iteraciones.

En 2D podemos utilizar el comando DensityPlot, que hace exactamente lo que queremos, usar un color según la densidad alcanzada por cada punto, según lo grande que sea su módulo en este caso.
Le aplicamos al módulo una exponencial negativa para reducir los valores y que entren en el rango computable por Mathematica. Lo de ColorFunction->Hue es sólo para que se use una gama de colores diferente a la usada por defecto; cuestión de gustos.

DensityPlot[Exp[-Abs[f[(x + y*I), Niter]]], {x, -2, 2}, {y, -2, 2},  ColorFunction -> Hue]

Y así, en tres líneas de código conseguimos esta representación del conjunto de Mandelbrot:

Nos podemos aventurar con el pseudo-3D cambiando el comando DensityPlot por Plot3D:

Plot3D[Exp[-Abs[f[(x + y*I), Niter]]], {x, -2, 2}, {y, -2, 2}, ColorFunction -> Hue]

Ya avisé que era aventurarse...

Es fácil modificar el código de la segunda línea para conseguir los conjuntos de Julia asociados a distitos valores de c, al menos aproximaciones bastante burdas (a decir verdad lo que vemos son los conjuntos de valores para los que las iteraciones convergen, lo que se denomina conjuntos de Fatou, pero bueno, sus fronteras serían los de Julia).

Para c=i

Niter = 7;
f[z_, it_] := If[it > 0, f[z, it - 1]^2 + I, z];
DensityPlot[Exp[-Abs[f[(x + y*I), Niter]]], {x, -2, 2}, {y, -2, 2}, ColorFunction -> Hue]

Para c=-0.5+0.7i

Niter = 10;
f[z_, it_] := If[it > 0, f[z, it - 1]^2 - 0.5 + 0.7*I, z];
DensityPlot[Exp[-Abs[f[(x + y*I), Niter]]], {x, -2, 2}, {y, -2, 2}, ColorFunction -> Hue]

Nótese cómo se ha tenido que reducir el número de iteraciones para no salirnos del rango de valores manejable por el comando DensityPlot. Esto se solucionará cuando afinemos el código e introduzcamos una condición de parada para la función iterativa... para ello nos ayudara el hecho de que en una función iterativa del tipo ${z}_{n+1}={{z}_{n}}^2+c$ que iniciemos con un cierto $z_0$, si nos tropezamos con un valor cuyo módulo sea mayor que $\left| c \right|$ y que 2, implica que tenemos asegurada la divergencia para dicho $z_0$, por lo que no hay que seguir iterando.

Pero eso será en otra entrada; de momento en tres líneas de código tenemos un instrumento para visualizar (toscamente) conjuntos de Julia y el conjunto de Mandelbrot... no está mal.