jueves , 24 septiembre 2020
¿Que lo origina? - La belleza de la complejidad.
¿Que lo origina? - La belleza de la complejidad.

Generación del fractal de Mandelbrot en Octave/Matlab

En este tutorial vamos a analizar el significado del conjunto de Mandelbrot y vamos a hacer un pequeño script en Octave/Matlab para su graficación.

Números complejos

Motivados por el deseo de poder escribir soluciones explícitas para polinomios de grado $n$ del tipo $$p_n(z)=c_0+c_1z+c_2z^2+c_3z^3+…+c_nz^n=\sum_{k=0}^n c_kz^k$$ los matemáticos desarrollaron el concepto de unidad imaginaria que axiomáticamente se define como $i=\sqrt{-1}$. De esta manera la raíz cuadrada de cualquier número negativo por ejemplo $-p$ puede expresarse en términos de $i$ como $$\sqrt{-p}=\sqrt{-1}\sqrt{p}=ip$$ Los matemáticos emplearon esta construcción para demostrar que todos los polinomios $p_n(z)$ del tipo expresado arriba tienen exactamente $n$ raíces del tipo $z=a+ib$ con $a, b \in \Re$, lo que constituye el teorema fundamental del álgebra (que tan solo por el nombre, nos da una idea de la importancia del teorema y de por qué los matemáticos optarían por inventar un número imposible con tal de hallar una demostración). En general, todos los números del tipo $z=a+ib$ son llamados números complejos y se define $a$ como la parte real del número $z$, en tanto que $b$ es la parte imaginaria de $z$.

Derivación del fractal de Mandelbrot de primeros principios.

El plano complejo

Normalmente es muy fácil representar los números reales en una recta numérica, sin embargo, ya que los números complejos están formados por dos componentes, los matemáticos optaron por representarlos como puntos en un plano cartesiano, como se ilustra en la siguiente imagen de wikipedia.

Representación del $z=a+ib$ en el plano complejo. Origen: Wikipedia.

Desde este punto de vista, cada número complejo $z$ es similar a un vector, incluso la suma de números complejos es equivalente a la suma vectorial ya que si $z_1=a+ib$ y $z_2=c+id$ entonces la suma de ambos se expresa como $$z=z_1+z_2=(a+c)+i(b+d)$$ que es equivalente a la suma de dos vectores cartesianos $$\overrightarrow{\rm (a,b)}+ \overrightarrow{\rm (c,d)} = \overrightarrow{\rm (a+c,b+d)}$$

Sin embargo, a diferencia de los vectores cartesianos, el producto de dos números complejos tiene una definición natural que surge simplemente de efectuar las operaciones algebraicas. Desarrollando el producto de la manera convencional, el producto de $z_1$ y $z_2$ toma la forma: $$z=z_1z_2=(ac-bd)+i(ad+bc)$$

Continuando con el tratamiento de $z$ cual si vectores cartesianos se tratase, si $z=a+ib$ se define la magnitud de $z$ como $$|z|=\sqrt{a^2+b^2}$$

El conjunto de Mandelbrot

El conjunto de Mandelbrot surge de considerar lo siguiente, si se define la regla iterativa: $$z_{n+1}=z_n^2+c$$ con $z_0=0$ y siendo $c$ todos los puntos del plano complejo, algunos puntos tendrán la tendencia a crecer ilimitadamente conforme las iteraciones progresan. Por otro lado, así como hay puntos que crecen ilimitadamente, habrá otros puntos que, independientemente de cuanto progresen las iteraciones, siempre estarán acotados.

Para aclarar este punto consideremos un ejemplo de ambos casos. Si $c=1$, entonces los primeros términos de la iteración son los siguientes: $$z_0=0$$ $$z_1=0^2+1=1$$ $$z_2=1^2+1=2$$ $$z_3=2^2+1=5$$ $$z_4=5^2+1=26$$ $$z_5=26^2+1=677$$ y así sucesivamente. Con esto ya podemos vislumbrar que el punto $c=1$ del plano complejo es un punto con tendencia a diverger.

Si ahora consideramos el caso de $c=-1$, entonces obtenemos lo siguiente: $$z_0=0$$ $$z_1=0^2-1=-1$$ $$z_2=(-1)^2-1=0$$ $$z_3=0^2-1=-1$$ $$z_4=(-1)^2-1=0$$ de donde podemos vislumbrar que el punto $c=-1$ del plano complejo jamás alcanzará una magnitud mayor que $1$.

En general el conjunto de Mandelbrot se define como el conjunto de todos los puntos $c$ del plano complejo para los que la iteración definida arriba no diverge. De aquí que podemos asegurar que $c=-1$ pertenece al conjunto de Mandelbrot en tanto que $c=1$ no pertenece. Lo relevante aquí es que esta operación debemos efectuarla para todos los puntos del plano complejo para determinar si los puntos pertenecen o no al conjunto. Resulta claro que calcular a mano punto por punto es una tarea imposible, por lo que necesariamente, para hallar una representación del conjunto con buena definición necesitamos de una herramienta que realice la comprobación punto por punto de manera automática.

Cabe mencionar que matemáticamente se sabe que todos los puntos para los que $|c|>2$ son puntos que divergen, por lo que, por la misma razón, cualquier $c$ para la que en algún punto de la iteración $|z_n|>2$ también serán puntos que diverjan y necesariamente no son parte del conjunto de Mandelbrot.

Algoritmo

Sobre estas bases, la implementación más sencilla para representar el conjunto de Mandelbrot consiste en identificar los puntos que NO pertenecen al conjunto de Mandelbrot. Para hacer esto seguimos los siguientes pasos:

  1. Asumimos inicialmente que todos los puntos podrían ser parte del conjunto de Mandelbrot.
  2. Ejecutamos el algoritmo de iteración: $z_{n+1}=z_n^2+c$ con $z_0=c$ y $c$ siendo los puntos del plano complejo.
  3. Determinamos aquellas $c$ para las que $|z_{n+1}|>2$ y los descartamos como candidatos. Los puntos $c$ que sobreviven son nuestra n-esima aproximación del conjunto de Mandelbrot.
  4. Continuamos iterando desde el paso 2 con aquellos puntos $c$ que aún son candidatos.

Idealmente el conjunto de Mandelbrot está formado por los puntos que sobreviven un número infinito de iteraciones. En la práctica tenemos que cortar el cálculo para algún valor $N$. En nuestro caso vamos a definir esta $N$ como aquella iteración para la que ya no observemos ninguna diferencia entre la cantidad de pixeles que pertenecían al Mandelbrot en la iteración anterior y los pixeles que pertenecen al Mandelbrot en la presente iteración.

Implementación

Como ya se mencionó, esta implementación se realizará en Octave, que es una plataforma equivalente a Matlab en un 99%, con la ventaja de que se trata de un software gratuito y de código libre que se puede descargar para cualquier plataforma de la página de Octave.

Para comenzar, definimos el número de pixeles de la imagen usando la variable lado, y definimos la imagen como una matriz cuadrada de tamaño lado*lado. Usamos este tamaño para definir un par de matrices inicializadas a cero donde almacenaremos los valores z de cada punto (zn) y los colores de la imagen del fractal (mandelbrot).

lado=720; 
tamano=[lado, lado];
zn=zeros(tamano);
mandelbrot=zeros(tamano);

Ya que sabemos que todos los valores $|z|>2$ divergen, es suficiente con graficar los puntos del plano en el intervalo $[-2,2]$ y $[-2i,2i]$ para la representación del fractal. Sabiendo esto, definimos con el siguiente código vectorizado una matriz c con las coordenadas del plano complejo en el citado intervalo y con el origen de coordenadas en el centro de la imagen.

paso=4/lado;
[X,Y]=meshgrid(1:lado,1:lado);
c=((X-lado/2)+(lado/2-Y)*i)*paso;

Una ventaja de Octave es que nativamente maneja números complejos y que identifica a i como la unidad imaginaria, por lo que la matriz c definida arriba será una matriz poblada con números complejos. Del mismo modo, al saber Octave que se trata de números complejos, todas las operaciones matemáticas que efectúe Octave con estas matrices automáticamente obedecerán las reglas de los números complejos sin necesidad de que nosotros debamos hacer ningún tratamiento especial, lo que facilita enormemente las operaciones.

Después listamos todos los pixeles como potenciales candidatos del conjunto de Mandelbrot como lo definimos en el paso 1 de nuestro algoritmo. Para esto usamos la función find() que encuentra los indices de todos los valores distintos de 0. Para esto simplemente sumamos 1 a la matriz c para asegurar que la función find() encontrara solo valores distintos de 0.

pixelesAntes=0;   
dentro=find(c+1);   %Paso 1
pixelesDespues=length(dentro); 
iteracion=0;

Después ejecutamos el bucle que efectuará las iteraciones. El bucle compara el número de pixeles que formaban parte del conjunto de Mandelbrot antes y despues de la iteración y prosigue solo si el número de pixeles experimenta cambios.

Las iteraciones se efectuan únicamente sobre los pixeles que hayan sobrevivido la iteración anterior (Paso 2). Esto se logra almacenando en dentro los índices de los puntos que aún son candidatos. Después de iterar, se descartan aquellos puntos que fallan en satisfacer que su magnitud sea menor que 2 (Paso 3). Finalmente, a aquellos puntos que hayan sobrevivido les asignamos un color en términos del número de iteraciones transcurridas con mandelbrot(dentro)=iteracion.

while( pixelesAntes!=pixelesDespues )
  iteracion=iteracion+1;
  pixelesAntes=pixelesDespues;   
  zn(dentro)=zn(dentro).^2+c(dentro);   %Paso 2.
  dentro=find(zn<2);                   %Paso 3.
  pixelesDespues=length(dentro);   
  mandelbrot(dentro)=iteracion;    
endwhile

Cuando el número de pixeles dentro del conjunto de Mandelbrot deja de cambiar (lo que básicamente implica que los detalles finos del fractal se vuelven menores que el tamaño del pixel) el ciclo se detiene. En este punto la matriz mandelbrot contiene una serie de curvas de nivel representando la n-esima aproximación del conjunto de Mandelbrot, cada una con un tono distinto de gris en el intervalo [0, iteracion] (el último valor tomado por la variable iteracion). Con esto es preferible renormalizar al intervalo [0,1] para facilitar su visualización

mandelbrot=mandelbrot/iteracion;

Ya en el intervalo apropiado, el fractal resultante se visualiza fácilmente en una escala de grises con el comando imshow().

figure(1)
imshow(mandelbrot);
Fractal de Mandelbrot resultante de 138 iteraciones con una imagen de 720×720 pixeles en el intervalo [-2,2]x[-2i,2i].

La imagen resultante se puede guardar fácilmente desde el menú archivo de la ventana.

Colorización

Un detalle atractivo para la representación del fractal es asignar colores diferentes en función del valor de gris de la imagen. Aunque esto puede hacerse desarrollando un algoritmo propio para la colorización de la imagen, la manera más sencilla es simplemente tomar la imagen anterior en niveles de gris y asignar colores a la imagen empleando la herramienta Mapa de Gradiente de Photoshop que nos provee un menú interactivo para la selección de los colores. Optativamente también podemos usar la herramienta Mapa de Degradado del editor gratuito Gimp con el mismo propósito.

Menú Mapa de Gradiente en Photoshop, útil para colorizar nuestro fractal.

Curiosidades

Existe una relación interesante entre el conjunto de Mandelbrot y los diagramas de bifurcación. La relación existe porque la ecuación que define la iteración del Mandelbrot ($z^2+c$) puede transformarse en la ecuación que define la aplicación logística $rx(1-x)$ empleando el cambio de coordenadas adecuado. Puedes aprender más sobre la aplicación logística y generar el mapa de bifurcación de la ecuación del Mandelbrot en esta entrada del blog. También checa esta entrada de Wikipedia para aprender sobre la transformación que relaciona ambos conceptos.

Relación entre la aplicación logística y el Mandelbrot. Origen: Wikipedia, Buddhabrot.

Código fuente

Si conoces GitHub, puedes hallar el script de Octave con los comandos de este ejercicio en el link:

GitHub

También puedes descargarlo directamente de aquí:

Descarga

Acerca de itsgaraet

Físico especializado en aplicaciones ópticas. En la prepa me llamó la atención aprender lenguaje C y desde entonces comencé a usar programación para resolver problemas que lo ameritaran. Después de aprender Fortran, Matlab y LabVIEW para automatizar mediciones, realizar simulaciones o calcular propiedades físicas, me di cuenta de que a menudo solo buscaba pretextos para resolver problemas empleando programación ya que la experiencia de definir un problema y plantear su solución es una experiencia que disfruto bastante. Esto me llevó a intentar profundizar en el campo del software, principalmente en lo relativo a Ciencia de Datos y Desarrollo Web. También soy amante de la fotografía, planear una nueva imagen y terminar de llevarla a la realidad en Photoshop siempre es algo interesante.

Checa también

Paquetes incluidos en Anaconda

Instalación de Python 3.7 para Ciencia de Datos con Anaconda

Tutorial que expone el procedimiento para instalar Python con Anaconda en Windows y Linux (Ubuntu).

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *