sábado, 18 de marzo de 2017

Más conjuntos de Mandelbrot


En la anterior entrada vimos por encima qué era el conjunto de Mandelbrot y cómo se podía generar el gráfico. Recordamos que todo se basaba en la función generadora z2+c. Se puede cambiar la función y repetir el proceso con otras. Estudiemos otras variantes:

  • (z+c)(z-c)
  • z3+c
  • z4+c

Para implementar estas nuevas funciones, dado que hay que hacer algunas operaciones un tanto tediosas, pensé en utilizar alguna librería Javascript que facilitase el trabajo con números complejos. Probé Mathjs, pero es demasiado lenta, así que se ha quedado todo como Javascript plano.


De nuevo, obtenemos unos bonitos fractales.



(z+c)(z-c)


z3+ c



z4+c






En esta página tenemos el código modificado:

jueves, 16 de marzo de 2017

El conjunto de Mandelbrot para 'dummies'




¿Quién no se ha sentido atraído por la belleza de esta figura? La representación gráfica del conjunto de Mandelbrot es el típico ejemplo de figura fractal y es, probablemente, la más popular.


Pero, ¿qué es exactamente, y cómo se llega a esta figura?

Definiendo el conjunto de Mandelbrot


Pertenecen al conjunto de Mandelbrot todos aquellos números complejos que cumplen una propiedad: que el valor de una sucesión determinada, que ahora veremos, para un determinado número, es convergente, es decir, no tiende a infinito.

La sucesión de la que hablamos es recursiva, y  es esta:

zn+1 = z2 + C

donde C es el número que estamos tratando de comprobar si pertenece al conjunto o no y z se va calculando recursivamente, esto es, tomando como entrada de cada iteración el valor obtenido anteriormente.

Recordemos que la aritmética de números complejos involucra operaciones con binomios del tipo a+bi, donde i es la unidad compleja. Recordemos que las operaciones en las que interviene i son un poco especiales, en este caso debemos recordar que i2 es -1.

Vamos a hacer unas cuentas con la ayuda del intérprete de Python, que tiene soporte de serie para números complejos (la unidad imaginaria se denota con j en Python).

Primero, definimos la función que aplicaremos de forma recursiva:
>>> def fm(z, c):
...     return z**2 + c

Ahora, vamos a probar qué pasa con algunos números al aplicarle la función. Empezaremos con el número 1 + i, en notación Python habría que escribirlo como 1 + 1j. Aplicaremos la función anterior cinco veces seguidas, en la primera iteración tanto z como c serán 1 + i y en las siguientes, el resultado de la función se usa como entrada para la siguiente. Es más difícil decirlo que hacerlo:
>>> z = c = 1+1j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(1+3j)
(-7+7j)
(1-97j)
(-9407-193j)
(88454401+3631103j)

Vemos que tras apenas 5 iteraciones, los valores se van disparando, así que, claramente, la sucesión de valores es divergente, se va a infinito.

Probemos con 0.5 + 0.3j:
>>> z = c = 0.5 + 0.3j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(0.66+0.6j)
(0.5756+1.092j)
(-0.36114864+1.5571104j)
(-1.79416445761-0.82469660658j)
(3.03890160806+3.25928267968j)

También parece que diverge, pero no tan rápido como el caso anterior.

Otro caso más, con el número 0.2 - 0.1j:


>>> z = c = 0.2 - 0.1j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(0.23-0.14j)
(0.2333-0.1644j)
(0.22740153-0.17670904j)
(0.220485371029-0.180367812122j)
(0.216081251188-0.179536927955j)

Tras cinco iteraciones, no parece que se vaya a infinito. Probemos otras cinco iteraciones:


>>> for x in range(10):
...     z = fm(z, c)
...     print z
... 
(0.214457598616-0.177589128054j)
(0.214454163201-0.176170675885j)
(0.214954481072-0.175561069755j)
(0.21538373972-0.175475277291j)
(0.215598582395-0.175589042903j)
(0.215651236743-0.175713497468j)
(0.215630222717-0.175785666083j)
(0.215595792549-0.175809404656j)
(0.215572598999-0.175807535868j)
(0.215563255771-0.175798574862j)

 Parece claro que este número, 0.2 - 0.1j, al aplicarle la función del principio, la sucesión que se obtiene es convergente.

Se dice que un número pertenece al conjunto de Mandelbrot cuando la sucesión es convergente, es decir, no se va a infinito. Se puede probar que en cuanto que uno de los resultado de la sucesión se sale del rango de ± 2 ± 2i, la sucesión será divergente, esto es, el número no pertenecería al conjunto. La forma más rápida de calcular esto es comprobando, para cada uno de los resultados intermedios, si el módulo del número complejo es > 2, o lo que es lo mismo, si la suma de los cuadrados de sus partes real e imaginaria es > 4. De esta última forma llegamos al mismo resultado ahorrándonos hacer una raíz cuadrada, que siempre es más costosa computacionalmente.

Sabiendo esto, vamos a empezar a comprobar todos los posibles números imaginarios que hay entre -2 - 2i y 2 + 2i, y viendo si la sucesión que se deriva de cada uno de ellos diverge. Obviamente, hay infinitos números en esta región, así que tenemos que ir calculando el mayor número de ellos con pequeños incrementos.

De la misma forma, tendremos que hacer unas cuantas iteraciones para saber si la sucesión diverge o no. Por ejemplo, para el número 0.5 + 0.3j que vimos antes tuvimos que hacer 5 iteraciones o calcular 5 términos de la sucesión antes de darnos cuenta de que realmente divergía.

Empecemos de forma sencilla: moviéndonos en pasitos de 0.1 y haciendo 3 iteraciones. Tendríamos esta representación gráfica:



Vemos que tenemos muy poco nivel de detalle, incrementemos paulatinamente el nivel de detalle (haciendo los incrementos más pequeños y aumentando el número de iteraciones para "cazar" más números que están fuera del conjunto), ahí van unas cuantas muestras:





En este punto, ya tenemos un detalle bastante fino, pero es probable que, como hemos hecho pocas iteraciones, haya números que pensemos que pertenecen al conjunto, pero en realidad, no pertenecen. Aumentemos las iteraciones.





Hasta ahora nos hemos limitado a colorear en rojo los números que, al aplicar la función, la sucesión diverge; en negro, los que converge.

El siguiente paso en detalle, y aquí es donde aparece la "magia" de este conjunto, es asignar diferentes colores a la velocidad con la que los números que no pertenecen al conjunto, al evaluar la sucesión, divergen. Esto es, si asignamos colores diferentes al número de iteraciones necesarias para determinar si el número no pertenece al conjunto, empezamos a obtener unos bonitos gráficos.

La gracia en este caso es encontrar una gradación de colores y asignación al número de iteraciones que quede bonita. Por ejemplo, una sencilla, basada en amarillos:



Un último ejemplo, con mayor resolución, utilizando una gama más o menos aleatoria de colores:


Implementación


Es bastante sencillo implementar una representación gráfica del fractal. Esta que se muestra está hecha en HTML Canvas y Javascript, usando "plain Javascript", sin plugins ni librerías.

Lo primero que debemos hacer es definir un canvas en nuestro HTML. También hemos agregado unos input para poder jugar con el incremento y el número de iteraciones a probar.


<canvas id="canvas"></canvas>
<form>
<p>Increment: <input name="increment" value="0.01"></p>
<p>Iterations: <input name="iterations" value="1000"></p>
<p><input type="button" value="Run!" id="run"></p>
<p id="message"></p>
</form>

A continuación, escribimos el código Javascript que hará los cálculos y representará los puntos en pantalla. El primer paso, es acceder e inicializar el elemento canvas:


var width = 600;
var height = 600;
var c = document.getElementById("canvas");
c.setAttribute('width', width);
c.setAttribute('height', height);
var ctx = c.getContext("2d");

Ahora definimos algunas constantes y funciones:
  • adjust_scale(i, j) => centra los puntos correctamente en el gráfico. Nuestro sistema de coordenadas es el mátemático (0, 0) en el centro, en un canvas el (0, 0) es la esquina superior izquierda.
  • get_color(nIter) => devuelve un color en función del número de iteraciones. En esta función se puede ser creativo y experimentar para probar diferentes gamas de colores y gradaciones. En nuestro ejemplo nos limitamos a una gama de rojos.
  • mandelbrot(z, c) => calcular un término de la sucesión aplicando la fórmula zn+1 = z2 + C, ya explicada.
  • in_mandelbrot_set(c, maxIter) => determina si un número C está en el conjunto de Mandelbrot utilizando un máximo de iteraciones. Devuelve un objeto con dos claves: si la sucesión es convergente y el número de iteraciones necesarias para llegar a dicha conclusión.
  • run() => lanza el gráfico.

El código completo es este:

<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<canvas id="canvas"></canvas>
<form>
<p>Increment: <input name="increment" value="0.01"></p>
<p>Iterations: <input name="iterations" value="1000"></p>
<p><input type="button" value="Run!" id="run"></p>
<p id="message"></p>
</form>

<script>
var width = 600;
var height = 600;
var c = document.getElementById("canvas");
c.setAttribute('width', width);
c.setAttribute('height', height);
var ctx = c.getContext("2d");

var mandelbrot_limit = 1.8;
var mandelbrot_max_module = 4;



function adjust_scale(i, j) {
    return [(width/(2*mandelbrot_limit))*i + (width/2), (height/(2*mandelbrot_limit))*j + (height/2)];
}

function get_color(nIter) {
    var factor = 50 + nIter;
    factor = parseInt(factor);
    if (factor > 255) factor = 255;
    return "rgb("+factor+",0,0)";
}

function mandelbrot(z, c) {
    var x = z[0]; 
    var y = z[1];
    var real = x*x - y*y + c[0];
    var imag = 2*x*y + c[1];
    return [real, imag];
}

function in_mandelbrot_set(c, maxIter) {
    var x = 0;
    var res = c
    while (x < maxIter) {
        res = mandelbrot(res, c);
        if ((res[0]*res[0] + res[1]*res[1]) > mandelbrot_max_module) {
            return {'status': false, 'nIter': x}; // No está en el conjunto, hemos tardado x interaciones en saberlo
        }
        x++;
    }
    return {'status': true, 'nIter': x};
}

function run() {
    ctx.fillStyle = "#FFFFFF";
    ctx.fillRect(0, 1, width, width);

    var mandelbrot_inc = parseFloat(document.forms[0].elements.increment.value);
    var mandelbrot_test_iterations = parseInt(document.forms[0].elements.iterations.value);
    var point_width = 1;
    if (mandelbrot_inc >= 0.0009) point_width = 2;
    if (mandelbrot_inc >= 0.009) point_width = 8;       
    if (mandelbrot_inc >= 0.09) point_width = 16;       

    document.getElementById("message").innerHTML = "Increment: "+mandelbrot_inc+", Iterations: "+mandelbrot_test_iterations;
    for(var i=-mandelbrot_limit; i<=mandelbrot_limit; i += mandelbrot_inc) {
        for (var j=-mandelbrot_limit; j<=mandelbrot_limit; j += mandelbrot_inc) {
            check = in_mandelbrot_set([i, j], mandelbrot_test_iterations);
            if (check.status) {
                ctx.fillStyle = "#000000";
            } else {
                ctx.fillStyle = get_color(check.nIter);
            }
            var point = adjust_scale(i, j);
            ctx.fillRect(point[0], point[1], point_width, point_width);
        }
    }
}

document.getElementById("run").addEventListener('click', run);

</script>