Fractales: Mandelbrot

El tercer y último ejemplo de fractales que daremos es el famoso fractal de Mandelbrot. Este fractal se basa en una ecuación y un método iterativo. La imagen que mostramos no es una gráfica de una ecuación sino un conjunto de valores.

Partimos de la ecuación: Zn+1 = Zn2 + C, donde Z y C son números complejos y n ≥ 0. Inicialmente, Z0 = C. Z es la variable, C una constante, y n el índice para la secuencia y por tanto para la iteración. Brevemente, los números complejos se componen de dos partes distintas: la parte real x y la parte imaginaria y. Podemos escribir un número complejo como z = x + iy, donde x e y son números reales e i es el número complejo definido como i = √-1. También se puede reescribir usando la siguiente notación: (x, y).

Para la ecuación de Mandelbrot, podemos reescribirla en la forma de sus componentes:

(xn+1,yn+1) = (xn2-yn2+xc, 2xnyn+yc)

Visto de otro modo:

xn+1 = xn2 - yn2 + xc,
yn+1 = 2xnyn + yc

donde Zn = (xn, yn) y C = (xc, yc).

El método para obtener el fractal, basada en esta ecuación, trata de averiguar si un número complejo cualquiera para C pertenece al conjunto de Mandelbrot. Si C pertenece al fractal, entonces se colorea de negro (por ejemplo), y si no, pues de blanco. La forma de saber si C pertenece o no al fractal de Mandelbrot es calculando la siguiente iteración y comprobando su distancia a un valor limitante, r. Esto es, |Zn| < r, donde |Zn| es la distancia o longitud del número complejo.

Esto se puede reescribir como:

|Zn| = |√(xn2+yn2)|< r

Generalmente, se usa un radio (valor limitante) r = 2 para el fractal de Mandelbrot.

Ahora que tenemos las condiciones impuestas, deberemos iterar el proceso. Sin embargo, necesitamos otro valor limitante para saber cuándo terminar. Este valor deberá ser grande para dejar suficiente tiempo al proceso con el fin de comprobar que el valor de Zn no se salga fuera del radio r. Podemos usar un nivel de tolerancia de 3.000, 10.000, o incluso 20.000. Pero la verdad es que este valor depende del área del fractal en que nos fijemos. Si magnificamos una parte del fractal es posible que perdamos algunos detalles porque nuestro nivel de tolerancia no era muy grande. La causa de esto es que algunos valores de Zn sean malinterpretados como valores no pertenecientes al conjunto de Mandelbrot, cuando en verdad sí pertenecen.

Cambio de Coordenadas

Según vimos en el capítulo 2, debemos realizar un cambio de coordenadas, ya que los sistemas de coordenadas que usamos no son idénticos. Como nuestro proceso usa la fórmula de Mandelbrot, nos encontramos en el sistema de coordenadas complejas. Sin embargo, podemos pasarnos al sistema cartesiano con facilidad. El problema está en que al terminar nuestros cálculos, debemos representar algo en la pantalla. Por este motivo debemos cambiar del sistema cartesiano al sistema gráfico de la pantalla.

En el caso de generar el fractal de Mandelbrot, iremos píxel a píxel comprobando si éstos pertenecen al conjunto de Mandelbrot o no. Para hacer esto, debemos conocer previamente las dimensiones de nuestra vista. Como dijimos en el capítulo 2, no podemos representar todos los valores, porque éstos son infinitamente muchos, cuando sólo disponemos de la pantalla que es finito: limitado. Por lo tanto, debemos escoger las dimensiones para los valores de X y para los valores de Y. Luego, debemos convertir tales valores a sus equivalentes en píxeles.

Digamos que queremos representar la imagen entre x = [-6, +6] e y = [-4, +4] a una resolución de 500 x 500. Necesitamos saber qué representa cada píxel (xp,yp) en términos de estas dimensiones, ya que cada pareja de valores corresponderá a un valor de C en nuestra fórmula de Mandelbrot. Siguiendo nuestro ejemplo, obtenemos que,

        |xui - xuf|     |-6 - 6|
dxup = ------------ = --------- =
         anchurap         500

       12 unidades
    = ------------- = 0,024 unidades/píxel
       500 píxeles

        |yui - yuf|     |-4 - 4|
dyup = ------------ = --------- =
         alturap          500

        8 unidades
    = ------------- = 0,016 unidades/píxel
       500 píxeles

En nuestro ejemplo, cada píxel representa 0,024 unidades en el eje X y 0,016 unidades en el eje Y del plano cartesiano. Los valores de dxup y dyup son nuestros incrementos para comprobar cada píxel para C en nuestra fórmula de Mandelbrot.

Algoritmo

De nuevo, separaremos el algoritmo para dibujar el fractal de Mandelbrot en dos partes principales:

  • La primera parte, Iniciar_Mandelbrot(), es el cálculo para cambiar de píxeles a unidades cartesianas y la invocación a la función iterativa, Mandelbrot().
  • La segunda parte es el algoritmo del proceso iterativo, Mandelbrot().

Para Iniciar_Mandelbrot(), el algoritmo es:

Iniciar_Mandelbrot()
  1. Inicializar valores para: xui, xuf, yui, yui, xpi, xpf, ypi, ypf,
    num_max_iteraciones, radio, C
  2. anchurau Fracta |xuf - xui|
  3. alturau ← |yuf - yui|
  4. anchurap ← |xpf - xpi + 1|
  5. alturap ← |ypf - ypi + 1|
  6. dxup ← anchurau / anchurap
  7. dyup ← alturau / alturap
  8. Bucle: xp ← xpi hasta xpf con incremento de 1
  9. Bucle: yp ← ypi hasta ypf con incremento de 1
  10. C.x ← xui + xp * dxup
  11. C.y ← yui - yp * dyup
  12. Color ← Mandelbrot( C, num_max_iteraciones, radio )
  13. PonPíxel( xp, yp, Color )
  14. Terminar
  • C puede almacenar una coordenada en unidades cartesianas - un número complejo.
  • Necesitamos dos bucles anidadas porque necesitamos recorrer y comprobar todos los píxeles en pantalla.
  • La función PonPíxel() hace referencia a una función de la biblioteca o API gráfica para establecer un color para un píxel determinado, según las coordenadas dadas.

Para Mandelbrot(), el algoritmo es:

Color Mandelbrot( C, max_iter, r )
  1. bEsMandelbrot ← Verdadero
  2. Z ← C
  3. Z_Sig ← Ecuación( Z, C )
  4. Bucle: k ← 1 hasta max_iter con incremento de 1 y
    mientras que bEsMandelbrot = Verdadero
  5. Si Distancia( Z_Sig ) < r entonces
  6. Z ← Z_Sig
  7. Z_Sig ← Ecuación( Z, C )
  8. Si no, entonces
  9. bEsMandelbrot ← Falso
  10. Si bEsMandelbrot = Verdadero, entonces
  11. Terminar( Negro )
  12. Si no, entonces
  13. Terminar( Blanco )

C, Z, y Z_Sig pueden almacenar una coordenada en unidades cartesianas, cada una.

En este algoritmo, tenemos dos condiciones para terminar el bucle:

  1. k > max_iter, y
  2. bEsMandelbrot = Falso

Esta última condición se basa en que la distancia de Z ≥ r. Es decir, Z se dispara hacia el infinito.

Al final, la función Mandelbrot() terminará devolviendo o bien el color negro, si el valor de C pertenece al conjunto de Mandelbrot, o bien blanco, si el valor de C no pertenece a dicho conjunto.

Para Ecuación(), el algoritmo es:

Complejo Ecuación( Z, C )
  1. Resultado.x ← Z.x*Z.x - Z.y*Z.y + C.x
  2. Resultado.y ← 2*Z.x*Z.y + C.y
  3. Terminar( Resultado )

Aplicamos la fórmula Zn+1 = Zn + C, pero descomponiéndola en partes reales e imaginarias para usarse en el plano cartesiano.

Para Distancia(), el algoritmo es:

Real Distancia( Z )
  1. Resultado ← Raíz_Cuadrada_de( Z.x*Z.x + Z.y*Z.y )
  2. Terminar( Resultado )

Observaciones

En este método, comprobamos cada píxel para saber si tal coordenada en el plano complejo-cartesiano corresponde al conjunto de Mandelbrot. Esto se realiza en base a su fórmula y según las condiciones limitantes. Este algoritmo utiliza métodos vistos en el trazado de una ecuación en el plano cartesiano. Sin embargo, en el trazado, sólo dibujábamos líneas que pertenecían a la línea. En el caso del fractal de Mandelbrot, debemos comprobar cada píxel.

En el algoritmo anterior, hemos tenido en cuenta la orientación del eje-Y de la pantalla (en píxeles). Esto lo hemos realizado con la siguiente fórmula:

xu = xui + xp * dxup,
yu = yui - yp * dyup

Fijémonos en el signo positivo para el valor de xu, y el signo negativo para yu. El signo positivo para x indica una orientación idéntica, mientras que el negativo de y indica una orientación contraria al sistema de coordenadas. Esto suele ser lo habitual en sistemas gráficos, pero si esto no es el caso, con simplemente cambiar los signos podemos elegir la orientación que nos sirva.

Otra observación es la complejidad del tiempo de ejecución del algoritmo. Con el algoritmo descrito anteriormente, observaremos una ralentización en su ejecución. Esto se debe a que el algoritmo sigue comprobando si cada valor permanece dentro de las restricciones dadas, hasta completar todas las iteraciones. Sin embargo, podríamos aplicar ciertos criterios según la "naturaleza" de la función base. Realizando comprobaciones de periocidad y dirección, podemos optimizar el proceso. El inconveniente se basa en que cada comprobación depende de la función en sí. Otro inconveniente está en que no todos los lectores saben manipular números complejos para implementar tales comprobaciones en sus programas.

Explicación Detallada

Fractales como los de Mandelbrot, Julia, y muchos más no son utilizables por la mayoría de las personas, que no tengan una buena base de matemáticas. Las imágenes basadas en tales fractales son estéticamente agradables. Otros usos se basan en la propiedad de autosemejanza. Hoy en día existen varios métodos de compresión de datos basados en fractales - esta cualidad de autosemejanza. Se puede alcanzar compresiones de hasta 5.000 á 1. Por supuesto, todo depende de la variedad de los datos al igual a poder encontrar la fórmula matemática para describir tales conjuntos de datos.

Ejercicios

Enlace al paquete de este capítulo.

  1. Escriba un programa que represente el fractal de Mandelbrot, con estas características:
    x = [-1, 1], e y = [-1, 1]
    Valor máximo de iteraciones: 3.000
    Radio: 4
    Resolución: 300 x 300
  2. Podemos usar otras funciones para crear otros tipos de fractales. Realice un programa que genere un fractal basado en cada una de las siguientes fórmulas, reemplazando el algoritmo Ecuación() en el apartado Algoritmo.

    Nota: Se aconseja usar alguna biblioteca para manipular números complejos, especialmente si no se sabe hacer usando matemáticas. En un programa de C++, se puede usar la plantilla complex<T> declarada en <complex>. Se sugiere usar la clase complex<double> para crear cada variable compleja y usar los operadores sobrecargados para llevar a cabo las operaciones necesarias: *, +, -, /, sin(), cos(), exp(), etc..

    1. Zn+1 = (conj(Zn))2 + C

      Nota: conj() se refiere al conjugado; esto es:
      z = x + iy  o  z = (x,y),
      su conjugado es:
      z' = x - iy  o  z = (x,-y),

      [IMAGEN]: Ejercicio 2a
      Ejercicio 2a
    2. Zn+1 = sen(Zn / C)
      [IMAGEN]: Ejercicio 2b
      Ejercicio 2b
    3. Zn+1 = C*exp(Zn)
      [IMAGEN]: Ejercicio 2c
      Ejercicio 2c
    4. Zn+1 = Zn3 + C
      [IMAGEN]: Ejercicio 2d
      Ejercicio 2d
    5. Zn+1 = C*Zn2
      [IMAGEN]: Ejercicio 2e
      Ejercicio 2e
    6. Zn+1 = C*senh(Zn)
      [IMAGEN]: Ejercicio 2f
      Ejercicio 2f
    7. Zn+1 = C*sen(Zn)
      [IMAGEN]: Ejercicio 2g
      Ejercicio 2g
    8. Zn+1 = C*exp(conj(Zn)) * C*exp(Zn) = C2*exp(2x)
      [IMAGEN]: Ejercicio 2h
      Ejercicio 2h
    9. Nota: x indica la parte real de Zn

    10. Zn+1 = sen(Zn) + Zn2 + C
      [IMAGEN]: Ejercicio 2i
      Ejercicio 2i
    11. Zn+1 = cosh(Zn) + Zn2 + C
      [IMAGEN]: Ejercicio 2j
      Ejercicio 2j
    12. Zn+1 = exp(Zn2) + C
      [IMAGEN]: Ejercicio 2k
      Ejercicio 2k
    13. Zn+1 = ln(Zn) + Z2 + C

      Nota: ln hace alusión al logaritmo neperiano o natural

      [IMAGEN]: Ejercicio 2l
      Ejercicio 2l

    Observación: Podemos crear cualquier función, especialmente tomando como modelo:
    Zn+1 = f(z) * Zn, donde f(z) es cualquier función trigonométrica compleja;
    por ejemplo,
    f(z) = cos(z),
    f(z) = tan(z),
    f(z) = 2*sen(z)*cos(z),
    etc.

    Advertencia: Es posible que, en algunas bibliotecas estándares, las implementaciones de algunas funciones no den buenos resultados. Por lo tanto, use las siguientes definiciones:

  3. Se aconseja crear un puntero a una función, para luego poder cambiar de función fácilmente, sin tener que reprogramar el algoritmo. Por ejemplo, el código parcial en C++, puede ser el siguiente:

    typedef Punto (*f_z)( Punto, Punto );
    
    Punto mandelbrot( Punto Z, Punto C );
    
    void fractal( f_z f, Punto Z_0, Punto C, unsigned long num_iteraciones, double r )
    \{
      ...
      for( k=2; k<=num_iteraciones && !bEsMandelbrot; k++ )
        if( (d=dist(Z_sig)) < r )
        \{
          Z = Z_sig;
          Z_sig = f( Z, C ); /* Invocamos la función general */
        }
        else bEsMandelbrot = true;
     ...
    }
    
    void Dibujar(void)
    \{
      ...
      fractal( mandelbrot, Z_0, C, num_iteraciones, r );
    }
    

    Podemos ver que sólo tenemos que pasar el puntero de la función que queramos mostrar. El algoritmo implementado en fractal() aplicará cualquier función dada.

  4. Las imágenes interesantes de los fractales se encuentran en la "costa" o borde del fractal. Aquí es donde se encuentran muchas bahías, ríos, afluentes, y deltas. Cree un programa para facilitar el engrandecimiento (o zoom) de cualquier área de un fractal. Esto se puede lograr, calculando las dimensiones originales de la imagen, un factor de engrandecimiento, y un punto central a tal área engrandecida.

    Por ejemplo,
    Si comenzamos con unas dimensiones de:
    x = [-1, +3] e y = [0, +2],

    entonces podemos hacer un zoom al área como punto central (0,5, -0,5) y con un factor de zoom de 3,0. Es decir, nuestro área nueva será 3 veces menor que la original, con un punto central de (0,5, -0,5). Esto se haría de la siguiente forma:
    1. Mudamos una de las esquinas de nuestro área rectangular al origen: (0, 0). Esto se puede calcular fácilmente:
      Elegimos la esquina inferior izquierda: (-1, 0). Ahora desplazamos las dimensiones, para que esta esquina se sitúe en el origen, (0, 0):
      x = [-1-(-1), +3-(-1)]
      y = [0-0, +2-0]

      Obtenemos,
      x = [0, +4]
      y = [0, +2]
    2. Ahora cambiamos las dimensiones según el factor de zoom: 3,0. Como se trata de un zoom hacia dentro, estamos reduciendo el tamaño original por 3,0, que es lo mismo que dividir las longitudes entre 3,0. Esto sería:
      x = [0/3,0, 4/3,0]
      y = [0/3,0, 2/3,0]

      Resultando en,
      x = [0,0, 1,3333]
      y = [0,0, 0,6666]
    3. Ya que las dimensiones han sido reducidas, ahora podemos desplazar las dimensiones de nuestro área según el punto central dado: (0,5, -0,5). Esto no es más que una suma:
      x = [0,0+0,5, 1,3333+0,5], e
      y = [0,0-0,5, 0,6666-0,5],

      obtenemos que,
      x = [+0,5, +1,8333]
      y = [-0,5, +0,1666]

      Recalculando el fractal con estas nuevas dimensiones, obtendremos una imagen de un área reducida 3 veces de la original y centrada en el punto (0,5, -0,5). Como el área es más pequeña que la original, los valores de dxup y dyup se ven reducidos. Esto implica que nuestra imagen permitirá mostrar más información al someter estas nuevas dimensiones a la misma resolución gráfica.
  5. Si el lector tiene experiencia en programar aplicaciones interactivas, entonces podemos crear un programa que permita al usuario elegir el área a investigar mediante la creación de un rectángulo pequeño en la imagen. Por ejemplo, al pinchar el botón izquierdo del ratón en la imagen, podría comenzar a crear un rectángulo arrastrándolo. Dicho rectángulo puede resultar al soltar tal botón izquierdo.

    Este rectángulo ya contiene las nuevas dimensiones de la imagen, para realizar el zoom. Lo único que tendríamos que hacer es convertir los valores de los píxeles a su representación en coordenadas (matemáticas) del plano complejo.
  6. Aún no hemos visto la forma de manipular colores; esto se trata en el siguiente capítulo 4. De todas formas, es una parte importante de la estética de la imagen. Un método popular se basa en dar un color a un píxel el cual representa un número complejo (coordenada) que no pertenece al conjunto del fractal. El criterio de dar un color u otro se basa en la distancia de tal coordenada desde su posición anterior (dentro del conjunto del fractal) a su posición actual (fuera del conjunto del fractal). Restringiendo tal distancia a un valor entre el intervalo: [0,00, 1,00], podemos usar un mapa de colores. Con este mapa podemos asignar o calcular un color, de acuerdo a nuestra gama de colores, según un valor entre el intervalo anterior. El tema de mapa de colores es tratado en el siguiente capítulo 4.

    Podemos hacer el siguiente cálculo para obtener un valor, a modo de parámetro para nuestro mapa de colores, en el intervalo [0,00, 1,00]:
    mapa( r / d ), donde d r
    d es la distancia entre las dos últimas coordenadas, y
    r es una constante que indica el radio del fractal.
    Por ejemplo,
    [IMAGEN]: Ejercicio 6
    Ejercicio 6

    Si el lector no tiene suficientes conocimientos para usar colores, entonces es recomendable pasar al siguiente capítulo 4, y luego volver a este ejercicio, cuando sepa la forma de manipular colores.