Caos y conjuntos fractales

F. Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Introducción

En la primera parte de esta práctica se estudia la manera de "compilar funciones" en Mathematica e investigamos las iteraciones de la función logística y su relación con el "caos matemático". Los ejemplos que se ponen deberían ser suficientes para que aprendas a "compilar funciones" sencillas. En las dos últimas partes deberás definir funciones compiladas para representar gráficamente conjuntos fractales obtenidos por iteración de funciones complejas.

Compilar funciones. La función logística. Caos matemático

Cuando definimos una función de la forma usual en Mathematica, por ejemplo f[x_] := x* Sin[1/x], dicha definición queda almacenada de forma que cuando le damos a x un valor particular Mathematica sustituye dicho valor en x*Sin[1/x]. El valor que le damos a x puede ser de muy variadada naturaleza; puede ser exacto o aproximado, entero, real o un complejo, una expresión simbólica o una lista... Mathematica debe tener en cuenta todas estas posibilidades con el consiguiente coste en tiempo de cálculo. Recuerda que en lenguajes como C o Fortran hay que declarar los tipos de variables que vamos a usar en cada proceso.  En Mathematica también podemos hacer esto usando "Compile[{{x,tipox},{y,tipoy},{z,tipoz},...},expr]" que crea una "función compilada" que contiene una sucesión de instrucciones muy simples escritas en un lenguaje próximo al código máquina y evalúa "expr" para valores dados de los argumentos x,y,z,... de los cuales se especifica en cada caso su tipo. Esta forma de proceder es útil para definir funciones que deben realizar una gran cantidad de cálculos númericos sencillos repetitivos o iterativos. Se trata de un procedimiento que Mathematica usa internamente en comandos que tienen que evaluar repetidamente una función en muchos valores como, por ejemplo, Plot[ ] o NIntegrate[ ]. Estos comandos usan la opción "Compile". Veamos un bonito ejemplo del uso de "Compile" con la función logística que también es un buen ejemplo de cómo se trabaja con listas en Mathematica.

La función logística

logistica = 
      Compile[{{a, _Real}, {x, _Real}}, a*x*(1 - x)]

Observa que las variables de la función se ponen a la derecha encerradas entre llaves; la expresión {{a,_Real},{x,_Real}} le dice a Mathematica que "logistica" es una función de dos variables que son de tipo "Real". Esta función se usa igual que cualquier otra función de dos variables.

logistica[ .2, 2.3]

Plot[logistica[3. , x], {x, 0, 1}];

Plot3D[logistica[a, x], {a, 0, 4}, {x, 0, 1}, 
        Ticks -> None];

Hemos definido la "función logística" logistica[a, x] = a x (1 - x) donde se considera que a es un parámetro con valores 0 < a ≤ 4 y la variable x se mueve en el intervalo [0, 1]. Vamos a hacer iteraciones de esa función para valores fijos del parámetro a en diferentes puntos del intervalo [0,1]. Fijafo un valor de a, la función x -> logistica[a, x] puede representarse en la forma logistica[a, #] &. Así evitamos tener que definir de forma explícita la función x -> logistica[a, x].

prmtrfijo[a_Real][m_List]:=Map[logistica[a, #] &, 
        m]

prmtrfijo[3.4][Range[0. , 1. ,  .01]]  (* 
      dividimos [0, 
            1] en 100 partes iguales y calculamos los valores de la función logística para 
a = 
        3.4 en cada punto de la partición obtenida *)

Para calcular las iteraciones podemos usar el comando NestList[f,z,n] cuya salida es una lista cuyos elementos son las iteradas sucesivas de f sobre z, empezando en el valor inicial z y terminando en (f^([n]))(z) (el resultado de aplicar f n veces al valor z). Nuestro valor inicial de partida será una lista de valores de x a la que aplicaremos la función prmtrfijo[ ] n veces. Nos vamos a quedar solamente con los últimos m elementos de la lista y le añadimos la lista inicial al principio.

iteraciones[a_Real, L_List, n_Integer, m_Integer] := 
      Prepend[Take[NestList[prmtrfijo[a], L, n], -m], L] /; 
        n > m

Aquí puedes ver cómo funciona.

\(\(\(tmp = 
      iteraciones[3.4, Range[0. , 1. ,  .2], 105, 5]

    tmp\  // 
      TableForm[#, 
          TableHeadings -> {{\*""(Subscript[x, 0])"", \*""(Subscript[x, 101])"", \*""(Subscript[x, 102])"", \*""(Subscript[x, 103])"", \*""(Subscript[x, 104])"", \*""(Subscript[x, 105])""}, {}}] &\)\(\ \)\( (*\ 
      además\ de\ la\ lista\ inicial, Range[0. , 1. ,  .2], 
      obtenemos\ otras\ 5\ listas\ con\ las\ iteraciones\ de\ órdenes\ 101, 102, 103, 
      104\ y\ 105\ sobre\ cada\ punto\ de\ la\ lista\ inicial\ para\ el\ valor\ a = 3.4\ *) \)\)\)

Vamos a transformar esta lista en otra apropiada para representarla gráficamente con ListPlot[ ]. Con el comando Transpose[ ] obtenemos una lista de listas con las iteradas de cada punto de la lista inicial.

tmp = Transpose[tmp]
\(tmp // 
      TableForm[#, 
          TableHeadings -> {{}, {\*""(Subscript[x, 0])"", \*""(Subscript[x, 101])"", \*""(Subscript[x, 102])"", \*""(Subscript[x, 103])"", \*""(Subscript[x, 104])"", \*""(Subscript[x, 105])""}}] &\)

Ahora convertimos cada lista en una lista de pares cuya primera componente es el primer elemento de dicha lista. Para ello aplicamos una regla de sustitución apropiada. Observa que tenemos que aplicar una regla de sustitución diferida. Después disponemos todos los pares obtenidos en una lista (Flatten[lista,1]).

tmp = 
      Flatten[tmp /. 
          x_?VectorQ :> Table[{x[[1]], x[[k]]}, {k, 2, Length[x]}], 
        1]

Finalmente, eliminamos los pares que se repiten y ordenamos de menor a mayor con Union[ ] y obtenemos una lista de pares cuyos primeros elementos son los de la lista inicial y los segundos elementos son sus iteradas de órdenes 101,102,103,104 y 105 (aunque no necesariamente en ese orden).

tmp = Union[tmp]//(NumberForm[#, 9] & )

Podemos combinar los comando anteriores en uno solo

listafinal[a_, n_, m_, p_] := 
      Union[Flatten[
          Transpose[iteraciones[a, Range[0. , 1. , N[1/p]], n, m]] /. 
            x_?VectorQ :> Table[{x[[1]], x[[k]]}, {k, 2, Length[x]}], 
          1]]

Comprobamos este comando.

(listafinal[3.4, 105, 5, 5] // NumberForm[#, 9] & ) == 
      tmp

Representamos la lista obtenida.

ListPlot[listafinal[3.4, 105, 5, 5], 
        PlotStyle -> {PointSize[0.015], Hue[0]},  
        PlotRange -> All];

Las sucesiones obtenidas, salvo la primera y la última que son la sucesión nula, son (a partir de un término en adelante) periódicas con período 2 y repiten los valores 0.451963 y 0.842154 indefinidamente; y este comportamiento es independiente del punto de partida. Para convencernos de que esto es siempre así vamos a repetir lo anterior haciendo más divisiones.

ListPlot[listafinal[3.4, 200, 100, 100], 
      PlotStyle -> {PointSize[0.01], Hue[0]},  
      PlotRange -> All];

La sucesión de las iteradas para el valor de a = 3.4 resulta ser (a partir de un término en adelante) una sucesión periódica con período igual a 2 independientemente del valor inicial de partida 0 < x < 1. Probemos a ver lo que pasa con otro valor próximo del parámetro a. Pongamos a = 3.5.

ListPlot[listafinal[3.5, 200, 100, 100], 
        PlotStyle -> {PointSize[0.01], Hue[0]}, 
        PlotStyle -> AbsolutePointSize[0.0001],  
        PlotRange -> All];

La sucesión de las iteradas para el valor de a = 3.5 resulta ser (a partir de un término en adelante) una sucesión periódica con período igual a 4 independientemente del valor inicial de partida 0 < x < 1. Veamos lo que pasa si cambiamos un poquito el valor de a y tomamos a = 3.7.

ListPlot[listafinal[3.7, 200, 100, 100], 
        PlotStyle -> {PointSize[0.01], Hue[0]}, 
        PlotStyle -> AbsolutePointSize[0.0001],  
        PlotRange -> All];

¡Sorpresa! La sucesión de iteradas ahora parece estar distribuida por todo el intervalo [0.25, 0.95]de una forma caótica independientemente del valor inicial de partida 0 < x < 1. Esto se llama "caos" matemático.

El diagrama de bifurcación

La siguiente celda contiene algunos comandos que he definido para representar el "diagrama de bifurcación" que muestra las sucesiones de iteradas para distintos valores de a. Se trata de una curiosidad matemática y no es mi propósito que estudies los comandos que siguen. Tómatelo como un juego. Ejecuta la celda y sigue adelante.

bifurcacion[n_Integer, m_Integer, p_Integer, q_Integer] := 
      Table[{a, Drop[iteraciones[a, Range[0. , 1. , N[1/p]], n, m], 1]}, {a, 
            0. , 4. , N[1/q]}] // Chop
diferentes[L_List] := 
      Module[{lista = {}}, 
        For[k = 1, k < Length[L], k++, 
          lista = Append[lista, 
              If[MemberQ[lista, L[[k]]] == False, 
                L[[k]], {}]]]; Flatten[lista]]
lista[n_Integer, m_Integer, p_Integer, q_Integer] := 
      Map[diferentes, 
        Partition[
          Flatten[bifurcacion[n, m, p, q], 3], (p + 1)*m + 
            1]]
puntos[n_Integer, m_Integer, p_Integer, q_Integer] := 
      Union[Flatten[
          lista[n, m, p, q] /. 
            x_?VectorQ :> 
              If[Length[x] > 1, 
                Table[{x[[1]], x[[k]]}, {k, 2, 
                    Length[x]}], {{x[[1]], x[[1]]}}], 1]]

Aquí puedes ver el diagrama de bifurcación.

ListPlot[puntos[75, 25, 15, 50], 
        PlotStyle -> {PointSize[0.001], Hue[0]},  
        PlotRange -> All, 
        AxesLabel -> {"a", "iteradas"}];

Lo que estás viendo es una "gráfica". Estamos representando en el eje de abscisas valores del parámetro a y para cada valor de a representamos en ordenadas los valores de la sucesión de iteradas correspondientes a dicho valor pero calculadas para distintos valores iniciales de \(TraditionalForm\`x ∈ \(\(\([0, 1]\) . \ 
            La\)\(\ \)\(gráfica\)\(\ \)\(nos\)\(\ \)\(dice\)\(\ \)\(que\)\(\ \)\)\)para valores del parámetro a ∈ [0, 1]las iteradas convergen a 0 (independiente del valor inicial de partida en [0,1]); para a∈[1,2.9] la sucesión de iteradas converge a un límite (el mismo independientemente del valor inicial de partida en [0,1]); para valores de a ∈ [3.1, 3.4] la sucesión de iteradas resulta ser (a partir de un término en adelante) una sucesión periódica con período igual a 2 (independientemente del valor inicial de partida 0 < x < 1). Si amplías la gráfica anterior podrás encontrar intervalos de valores de a para los que la sucesión de iteradas resulta ser (a partir de un término en adelante) una sucesión periódica con período igual a 4 (independientemente del valor inicial de partida 0 < x < 1). Para valores de a > 3.7el comportamiento de la sucesión, independientemente del valor inicial de partida 0 < x < 1, es caótio.

Compilar funciones

Compilar la función logística ha permitido hacer los cálculos en muy poco tiempo.

iteraciones[3.4, Range[0. , 1. ,  .01], 500, 100]; // 
      Timing   

    iteraciones[3.4, Range[0. , 1. ,  .001], 500, 100]; // 
      Timing

Si no compilamos el tiempo necesario es casi el doble.

\(\(\(Clear[logistica]

    logistica[a_, x_] = a*x*\((1 - x)\); \  (*\ 
      definimos\ la\ función\ sin\ compilar\ *) 
\(\(iteraciones[3.4, Range[0. , 1. ,  .01], 500, 
          100];\)\(//\)\(Timing\)\(\ \)\( (*\ 
        estamos\ usando\ la\ función\ sin\ compilar\ 100  x500 = 
          50000\ veces\ *) \)\)

    \(iteraciones[3.4, Range[0. , 1. ,  .001], 500, 100];\) // 
      Timing\)\(\ \)\( (*\ 
      estamos\ usando\ la\ función\ sin\ compilar\ 1000  x500 = 
        500000\ veces\ *) \)\)\)

Clear["Global`*"]

"Compile" supone por defecto que sus argumentos son números reales con precisión máquina. Podemos compilar cuatro tipos de variables _Integer, _Real, _Complex y True | False. Aquí tienes un ejemplo con una variable booleana.

sumdifcond=Compile[{{x, 
            True | False}, {y, _Integer}, {z, _Integer}}, 
        If[x, z + y, z - y]]

sumdifcond[True, 5, 8]

Los siguientes ejemplos permiten apreciar la diferencia entre compilar o no una función. El primer comando calcula la suma de Riemann de una función f en un intervalo [a,b] que se divide en n partes iguales evaluando la función en el punto medio de cada intervalo. El segundo comando hace lo mismo con la diferencia de que se compila la función antes de evaluarla y también se compila el propio comando. Aplicamos estos comandos a la función g(x) = (1/x) cos(1/x) en el intervalo [.01,1] al que dividimos en 100000 partes iguales. Observa la notable diferencia en los tiempos de ejecución.

f[x_] = Cos[1/x]/x;
fc = Compile[{{x, _Real}}, Cos[1/x]/x];
Sriemann[a_, b_, n_] := 
        Module[{inc = N[(b - a)/n]}, 
          Sum[f[a + 0.5*inc + k*inc], {k, 0, n - 1}]*
            inc];
SCPriemann = 
        Compile[{{a, _Real}, {b, _Real}, {n, _Integer}}, 
          Module[{inc = N[(b - a)/n]}, 
            Sum[fc[a + 0.5*inc + k*inc], {k, 0, n - 1}]*
              inc]];
Sriemann[0.01, 1. , 100000] // Timing
SCPriemann[0.01, 1. , 100000] // Timing

El uso de "Compile" está especialmente indicado cuando hay que hacer operaciones repetitivas sencillas que no se prestan al manejo de listas. La siguiente función calcula el primer natural n tal que la suma 1 + 1/2 + 1/3 +  ...  1/n es mayor que un número positivo dado x. Como esas sumas crecen muy lentamente hay que realizar un gran número de cálculos. Observa la gran diferencia entre compilar el código o no hacerlo.

primerasuma[x_] := 
      Module[{n = 1, s = 1. }, 
        While[s ≤ x, n = n + 1; s = s + N[1/n]]; n]   

    CPprimerasuma = 
          Compile[{{x, _Real}}, 
            Module[{n = 1, s = 1. }, 
              While[s ≤ x, n = n + 1; s = s + N[1/n]]; 
              n]] ;

primerasuma[12] // Timing
CPprimerasuma[12] // Timing

Los conjuntos de Julia y de Mandelbrot

Los fractales más conocidos son los conjuntos de Julia y de Mandelbrot. Es muy fácil generar dichos conjuntos usando números complejos. Para cada número complejo c definamos la función \(\(f\_c\) : ℂ⟶ℂ\) por (Subscript[f, c])(z) = (z^2) + c para todo z  ∈  ℂ. El conjunto de Julia  asociado al número complejo c se define como el conjunto J(c) de los números complejos z tales que la sucesión definida por (Subscript[z, 0]) = z,  
          (Subscript[z, n]) = (Subscript[z, n - 1]^2) + c = (Subscript[f, c])((Subscript[z, n - 1])) está acotada (algunos autores llaman conjunto de Julia a la frontera de J(c)). El conjunto de Mandelbrot es el conjunto de los números complejos c tales que la sucesión definida por (Subscript[z, 0]) = 0,  
          (Subscript[z, n]) = (Subscript[z, n - 1]^2) + c = (Subscript[f, c])((Subscript[z, n - 1])) está acotada. En otros términos, el conjunto de Mandelbrot es el conjunto de los números complejos c tales que 0 ∈ J(c). 
Para representar la sucesión (Subscript[z, 0]) = z,  
          (Subscript[z, n]) = (Subscript[z, n - 1]^2) + c = (Subscript[f, c])((Subscript[z, n - 1])) suele usarse la notación (Subscript[z, n]) = (Subscript[f, c]^((n)))(z), con el convenio de que (Subscript[z, 0]) = (Subscript[f, c]^((0)))(z) = z. Esta sucesión se llama la órbita de z por (Subscript[f, c]) y sus elementos se obtienen iterando sucesivamente la función (Subscript[f, c]) partiendo del valor inicial z. La notación (Subscript[f, c]^((n)))representa la composición de la función (Subscript[f, c]) consigo misma n veces.
Es sabido que si |c|≤2  y z∈ J(c), entonces la órbita de  z por (Subscript[f, c]) está contenida en el disco cerrado de centro 0 y radio 2 y, en particular, |z|≤2. También es sabido que el conjunto de Mandelbrot está contenido en el disco cerrado de centro 0 y radio 2 y si un punto c está en el conjunto de Mandelbrot, entonces  la órbita de  0 por (Subscript[f, c]) está contenida en el disco cerrado de centro 0 y radio 2.

Ejercicio 1

Ejercicio 2

Aproximación de raíces complejas por el método de Newton

Como seguramente ya sabes, el método de Newton es una técnica frecuentemente usada para calcular soluciones de una ecuación f(z) = 0. Para ello, se elige un valor inicial (Subscript[z, 1]) y, a partir de él, se forma la sucesión definida por (Subscript[z, n + 1])=(Subscript[z, n])-((f((Subscript[z, n])))/(f ' ((Subscript[z, n])))). Si dicha sucesión converge, su límite es una solución de la ecuación f(z) = 0. Naturalmente, cuando hay más de una solución, la elección del valor inicial determina a cuál de ellas converge el proceso. 
Mathematica tiene un comando "FindRoot[f(z)=0,{z,a}]" que utiliza precisamente el método de Newton para calcular soluciones de  la ecuación f(z) = 0 a partir del valor inicial z = a. Evalúa la siguiente celda.

FindRoot[(z^2) + 1 == 0, {z, 0.5}]

¿Sabes por qué ha fallado el método en este caso particular? Fíjate en cómo se van obteniendo las aproximaciones sucesivas (Subscript[z, n + 1])=(Subscript[z, n])-((f((Subscript[z, n])))/(f ' ((Subscript[z, n])))). Es evidente que, para f(z)=(z^2)+1 partiendo de un valor real obtenemos siempre valores reales por lo que el método no puede proporcionarnos las raíces de (z^2) + 1  = 0 que son complejas. Repite con una ligera variante.

FindRoot[(z^2) + 1 == 0, {z, 0.5 + 0.1  I}]

Ahora Mathematica ha calculado la solución z = i con gran precisión. Observa ahora el efecto de cambiar un poco el valor inicial.

FindRoot[(z^2) + 1 == 0, {z, 0.5 - 0.1  I}]

Ejercicio 3

Como puedes ver Mathematica aplica el método de Newton para calcular raíces tanto reales como complejas. Ahora bien, la intrepretación del método en el caso del cálculo de raíces reales no tiene sentido en el caso complejo. En el caso real, el punto (Subscript[z, n + 1]) se obtiene como el punto de corte de la recta tangente a la curva y = f(z) en el punto ((Subscript[z, n]), f((Subscript[z, n]))) con el eje de abscisas y se conocen condiciones precisas que garantizan la convergencia del método. Lo que ocurre en el caso complejo puede ser mucho más complicado como vas a ver.

Ejercicio 4