jueves , 6 agosto 2020
¿Sabes que significa? - Diagrama de bifurcacion
¿Sabes que significa? - Diagrama de bifurcacion

Aplicación logística y diagrama de bifurcación

En este tutorial vamos a analizar una manera de graficar el diagrama de bifurcación de la aplicación logística. Tradicionalmente la aplicación logística surge de considerar la ecuación $$x_{n+1}=rx_n(1-x_n)$$

Esta ecuación se planteo originalmente como una manera de describir el comportamiento de una población. Por ejemplo, si una población en un momento dado esta descrita como $x_n$ y $r$ es una taza de crecimiento, entonces la población en un momento posterior en el tiempo $x_{n+1}$ crecerá en proporción a la población actual $rx_n$. Ya que en la mayoría de los escenarios existe un límite para que tanto puede crecer una población (limitaciones de recursos, espacio limitado, competición entre los miembros de la misma población, etc.) el factor $(1-x_n)$ representa que la población no puede crecer mas allá de un máximo poblacional, que en esta ecuación esta representado por $1$. Se entiende que $x_n$ representa una población normalizada respecto al máximo poblacional, por lo que naturalmente $x_n \in [0,1]$. El modelo pretendía que para poblaciones pequeñas el crecimiento fuera prácticamente $rx_n$ (es decir, el crecimiento de la población es prácticamente libre y todos los integrantes de la población dejan descendencia) pero que la población se estanque al alcanzar el límite de población.

Conviene destacar que este modelo nunca pretendió ser una descripción exacta del comportamiento de una población, simplemente era una manera sencilla de representarlo, sin embargo Feigenbaum encontró en 1975 que las conclusiones obtenidas para este modelo eran válidas no solo para este, sino para múltiples modelos que solo tenían en común estar descritos por una ecuación iterativa del tipo: $$x_{n+1}=f(x_n)$$ donde $f(x)$ es una función con un solo máximo cuadrático parametrizada por el parámetro de bifurcación $r$ (que para los propósitos de este post, es simplemente una constante). Feigenbaum también descubrió que todas las funciones $f(x)$ que cumplían las condiciones anteriores daban origen a las ahora llamadas constantes de Feigenbaum $\delta= 4.669201609102990671853203821578… $ y $\alpha= 2.502907875095892822283902873218… $ que en teoría de bifurcaciones, se consideran constantes de importancia análoga a la constante $\pi$ de geometría y a la constante $e$ del cálculo.

En este post vamos a aprovechar el resultado de Feigenbaum para graficar primero el diagrama de bifurcación de la aplicación logística tradicional ($x_{n+1}=rx_n(1-x_n)$), y después, bajo la premisa de que el mismo resultado es aplicable para todas las relaciones $x_{n+1}=f(x_n)$, vamos a usar el mismo código para generar los diagramas de bifurcación de otras funciones conocidas. Esto nos ayudará a cimentar la idea de que los diagramas de bifurcación surgen espontáneamente de multitud de funciones y no son en absoluto exclusivos de la función $x_{n+1}=rx_n(1-x_n)$.

Convergencia de la iteración

Para entender el significado del diagrama de bifurcación, comencemos por graficar primero las iteraciones de $x_{n+1}=rx_n(1-x_n)$ para $x_n \in [0,1]$ y $r \in [0,4]$. Aprovecho para aclarar en este punto que para $r>4$, $x_n$ diverge. Por otro lado $r<0$ no tiene un significado físico, por lo que no lo consideraremos (el lector puede intentar el ejercicio de graficar para $r<0$, generar los diagramas de bifurcación y ponderar su posible significado). Para hacer esta gráfica, primero vamos a partir de un valor de $x_n$ pequeño (una población inicial poco saturada, lo que permite su libre expansión y observar su evolución en el tiempo).

import matplotlib.pyplot as plt
rs=[0.5*i for i in range(1,8)]  #Valores de r
rs.extend([3.75, 3.9])  #Agregamos un par de valores extra.
n=100   #Graficamos 100 puntos para cada r.
x0=0.01   #Valor inicial de xn

fig=plt.figure(figsize=(15,7))
ordenadas=range(n)
for j, r in enumerate(rs):
    xns=[x0]
    for i in range(n-1):
        xns.append(r*xns[i]*(1-xns[i]))
    plt.subplot(int('33'+str(j+1)))
    plt.plot(ordenadas, xns, label=f'r = {r}')
    plt.xlabel('Iteraciones')
    plt.ylabel('$x_n$')
    plt.legend()

plt.subplots_adjust(hspace=0.5, wspace=0.25)
plt.show()

Probando con los valores de $r=0.5, 1, 1.5, 2, 2.5, 3, 3.5, 3.75, 3.9$ obtenemos los siguientes comportamientos para $x_0=0.01$.

Comportamiento de la iteración de $x_n$ partiendo de $x_0=0.01$ para múltiples valores de $r$.

A groso modo, para valores de $r \leq 1$, la población $x_n$ tiende a desaparecer. Para valores $1<r<3$ la población experimenta un rápido crecimiento seguido de una estabilización hacia un valor fijo. Para $r>3$ no hay estabilidad pero $x_n$ se mantiene oscilando de manera más o menos periódica alrededor de la misma serie de valores.

Parecería razonable suponer que la elección de $x_0$ cambiaría dramáticamente como evoluciona la población. Para comparar, ahora probamos con $x_0=0.95$ (un valor muy cercano al límite poblacional) con lo que obtenemos el siguiente comportamiento para los mismos valores de $r$.

Comportamiento de la iteración de $x_n$ partiendo de $x_0=0.95$ para múltiples valores de $r$.

Curiosamente, de nuevo se observa que para $r \leq 1$ la sucesión decae a 0. Para $1<r<3$ la sucesión converge a un valor fijo que es igual al valor al que convergen las gráficas para $x_0=0.01$. Para $r>3$ de nuevo se observa que $x_n$ se mantiene oscilando de manera semiperiódica y que los valores entre los que oscila la sucesión son aparentemente similares a los valores entre los que oscila la sucesión para $x_0=0.01$.

Este comportamiento nos permite vislumbrar que de hecho, la sucesión converge a un valor o valores fijos y que los valores a los que la sucesión converge son independiente del valor inicial de $x_0$. A partir de este punto asumiremos que esto es cierto, el lector puede usar el código provisto para generar las sucesiones con otros valores de $x_0$ y verificar que sus respectivas sucesiones convergen al mismo valor.

Diagrama de bifurcación de $rx_n(1-x_n)$

Sabiendo que las sucesiones convergen a un valor fijo independiente de $x_0$, resulta válido preguntarnos hacia que valor o valores converge la sucesión si conocemos el valor de $r$. Para responder esta incógnita, vamos a graficar en función de $r$ el valor o valores hacia los que converge la sucesión. Para esto, vamos a asignar a $x_0$ un valor cualquiera ($x_0=0.01$) y vamos a calcular 1,000 términos de la sucesión. Luego vamos a tomar los últimos 50 términos y los colocaremos en una misma gráfica.

Estaremos tomando 50 términos y no solo el último elemento de la sucesión porque, de lo observado en nuestras anteriores gráficas, para $r>3$ la sucesión parece oscilar de manera periódica o semi-periódica entre una serie de valores y queremos registrar todos estos valores en nuestro gráfico. Es claro que cuando la sucesión converge a un solo valor, los 50 términos caerán en el mismo punto del gráfico, pero cuando la sucesión comience a oscilar entre varios valores, lo mejor será haber guardando los suficientes puntos para que el gráfico nos muestre un panorama completo de la dinámica de la sucesión.

x0=0.01  #Población inicial
n=1000  #Numero total de iteraciones
ultimos=50   #Términos que consideraremos.

r=0
fig=plt.figure()
while r<4: 
    xns=[x0]   #Lista donde se guardan los elementos de la sucesión.
    for i in range(n):
        xns.append(r*xns[i]*(1-xns[i]))
    r+=0.005
    xfinales=xns[n-ultimos:]    #Ultimos 50 valores de la sucesión.
    rs=[r for xs in xfinales]    #Lista con r repetido 50 veces.
    plt.scatter(rs, xfinales, c='b', marker='o', s=0.1)
    
plt.xlabel('$r$')
plt.ylabel('$x_n$\'s convergentes')
plt.show()

Con lo que obtenemos el gráfico:

Resultado de graficar los valores hacia los que $x_n$ converge para $r \in [0,4]$.

El gráfico resultante es tradicionalmente conocido como diagrama de bifurcación de la aplicación logística y nos permite apreciar que, dependiendo del valor de $r$, la sucesión puede converger a 1, 2, 4 o más valores. El gráfico también nos permite vislumbrar que la impresión que tuvimos antes acerca de que la sucesión se mantenía oscilando de manera periódica entre los mismos valores era correcta. Para valores de $r$ cercanos a 4, la cantidad de valores entre los que oscila la sucesión es tan grande que en la práctica $x_n$ nos da la impresión de saltar de punto a punto de manera caótica, aunque el presente gráfico nos permite apreciar que existe un orden subyacente.

Definición de una clase

Como se mencionó antes, Feigenbaum demostró que un comportamiento similar puede esperarse para $$x_{n+1}=f(x_n)$$ Si este es el caso, podemos seguir un procedimiento análogo para generar los diagramas de bifurcación de otras funciones $f(x)$.

Ya que queremos un código que funcione para funciones $f(x)$ en general, en vez de escribir nuevo código para cada función que grafiquemos, optaremos por definir una clase que haga todo el trabajo difícil y genere las gráficas deseadas con solo indicarle la función que queremos graficar. Esto lo podemos hacer con la siguiente clase.

class secuenciaLogistica:
    """
    Modela el sistema dinámico descrito por: x_{n+1} = f(r, x_n)
    """
    def __init__(self, fx):
        """
        Inicializa con x y r = 0 para una función fx arbitraria
        """
        self.x, self.r, self.fx = 0, 0, fx
    
    def actualizar(self):
        """
        Aplica el mapeo de acuerdo a la función fx.
        Se asume que fx es una función de 2 argumentos: r y x.
        """
        self.x = self.fx(self.r, self.x)
    
    def generarSecuencia(self, n):
        "Genera una secuencia de longitud n para la función fx"
        xns=[]
        for i in range(n):
            xns.append(self.x)
            self.actualizar()
        return xns
    
    def diagBifurcacion(self, inicio, fin, paso, x0, n, ultimos, titulo='', aspecto=(10,7)):
        """
        Genera el diagrama de bifurcación de fx.
        r en el intervalo [inicio, fin] dando pasos de tamaño `paso`
        Usa una sucesión de `n` elementos y toma los 
        `ultimos` elementos de la sucesión para incluirlos en la gráfica.
        Asigna un `titulo` a la gráfica.
        Emplea una relacion de aspecto `aspecto`.
        """
        import matplotlib.pyplot as plt
        self.r=inicio
        fig=plt.figure(figsize=aspecto)
        while self.r<fin:
            self.x=x0
            xns=self.generarSecuencia(n)
            self.r+=paso
            xfinales=xns[n-ultimos:]
            rs=[self.r for xs in xfinales]     
            plt.scatter(rs, xfinales, c='b', marker='o', s=0.1)
        plt.title(titulo)
        plt.xlabel('$r$')
        plt.ylabel('$x_n$\'s convergentes')
        plt.show()

Con la que podemos generar los diagramas de bifurcación de cualquier $f(r,x)$ con que se inicialice la clase.

Diagramas de bifurcación de $x_{n+1}=f(x_n)$

Probemos nuestra nueva clase para reproducir el diagrama de bifurcación anterior, esto lo hacemos de la siguiente manera.

aplLogistica = lambda r, x: r*x*(1-x)   #Define la función.
secOriginal=secuenciaLogistica(aplLogistica)   #Inicializamos la clase
secOriginal.diagBifurcacion(0,4,0.005,0.01,1000,50, titulo='$f(x_n)=rx_n(1-x_n)$')   #Llamamos al método que grafica el diagrama de bifurcación.
Diagrama de bifurcación de $f(x_n)=rx_n(1-x_n)$ generado usando la clase secuenciaLogistica.

Interesántemente, 3 líneas de código nos permitieron replicar todos los resultados anteriores. Aplicando la misma clase para otras funciones $f(x)$ representativas obtenemos los siguientes resultados.

Regla de iteración del fractal de Mandelbrot $f(x_n)=x_n^2+c$
reglaMandelbrot = lambda c,z: z**2+c 
secMandelbrot=secuenciaLogistica(reglaMandelbrot)
secMandelbrot.diagBifurcacion(-2,1/4,0.005,0,1000,50, titulo='$f(x_n)=x_n^2+c$')
Diagrama de bifurcación de $f(x_n)=x_n^2+c$ generado usando la clase secuenciaLogistica.

Puedes hallar información adicional sobre este fractal, incluyendo el código fuente para graficarlo en la entrada del blog: Generación del fractal de Mandelbrot en Octave/Matlab

Regla de iteración $f(x_n)=rx_n(1-tanh(x_n))$
from numpy import tanh
reglaTanh = lambda r,x: r*x*(1-tanh(x)) 
secTanh=secuenciaLogistica(reglaTanh)
secTanh.diagBifurcacion(0,17,0.01,0.01,1000,50, titulo='$f(x_n)=rx_n(1-tanh(x_n))$')
Diagrama de bifurcación de $f(x_n)=rx_n(1-tanh(x_n))$ generado usando la clase secuenciaLogistica.
Regla de iteración $f(x_n)=sin(rx_n)$
from numpy import sin
reglaSin = lambda r,x: sin(r*x) 
secSin=secuenciaLogistica(reglaSin)
secSin.diagBifurcacion(-10,0,0.005,0.01,1000,50, titulo='$f(x_n)=sin(rx_n)$')
Diagrama de bifurcación de $f(x_n)=sin(rx_n)$ generado usando la clase secuenciaLogistica.

Adviértase que una vez que la instancia de la clase está definida, podemos emplear los métodos de la clase con total libertad. Por ejemplo, podemos ampliar el rango de la gráfica anterior simplemente invocándola con nuevos parámetros.

secSin.diagBifurcacion(-5,5,0.0025,0.01,1000,50,titulo='$f(x_n)=sin(rx_n)$', aspecto=(14,7))
Diagrama de bifurcación de $f(x_n)=sin(rx_n)$ en el rango ampliado $[-5,5]$.

Conclusiones

Aunque normalmente los diagramas de bifurcación se asocian con la aplicación logística, los resultados anteriores nos permiten apreciar que este tipo de diagramas surgen de funciones $f(x)$ bastante arbitrarias. Naturalmente, el rango en que cada función se puede iterar sin diverger es particular de cada función, por lo que cada $f(x)$ ha de analizarse caso por caso para decidir el rango sobre el que generar un diagrama de bifurcación es posible.

Otra conclusión importante de estos diagramas es que nos permiten vislumbrar que incluso para reglas iterativas relativamente simples, el comportamiento de las sucesiones puede experimentar una dinámica caótica. Este argumento a menudo se emplea como forma de insinuar que, independientemente de cuan profundo sea nuestro conocimiento de un sistema físico, nuestra habilidad para calcular el estado final de dicho sistema físico siempre estará limitada ya que el comportamiento de los sistemas es intrínsecamente caótico. Esto puede verse como una expresión del popularmente conocido como Efecto Mariposa y nos ilustra que pequeños cambios en las condiciones iniciales (o en nuestro caso, cambios en el valor de $r$) pueden resultar en condiciones finales completamente diferentes.

Código fuente

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

GitHub

También puedes descargarlo directamente de aquí:

Descarga

Versión extendida con gráficas interactivas (no muy eficiente):

GitHub

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

Estima Pi con random() en Python

¿Podemos calcular $\pi$ con números aleatorios?

Tutorial en Python (Jupyter) que ilustra de manera gráfica el procedimiento para obtener una aproximación de Pi por el Método de Montecarlo. Se incluye el código fuente.

Deja un comentario

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