¿Cómo crece una población?

El modelo más simple de crecimiento poblacional de organismos es $\frac{dx}{dt}=rx$, donde $x(t)$ es la población en el tiempo $t$ y $r>0$ es la tasa de crecimiento.

Este modelo predice crecimiento exponencial $x(t)=x_0e^{rt}$ (solución de la ecuación diferencial) donde $x_0=x(0)$ es la población inicial. ¿Es esto válido?

Recordar que $\lim_{t\to\infty}x(t)=x_0\lim_{t\to\infty}e^{rt}=\infty$.

Este modelo no tiene en cuenta entonces sobrepoblación ni recursos limitados.

En realidad la tasa de crecimiento no es una constante, sino que depende de la población $\frac{dx}{dt}=\mu(x)x$. Cuando $x$ es pequeña $\mu(x)\approx r$, como antes, pero cuando $x>1$ (población normalizada) $\mu(x)<0$: la tasa de muerte es mayor a la tasa de nacimiento. Una forma matemática conveniente de modelar lo anterior es con una tasa de crecimiento $\mu(x)$ decreciendo linealmente con $x$.

Referencia:

Strogatz, Steven. NONLINEAR DYNAMICS AND CHAOS, ISBN: 9780813349107, (eBook disponible en biblioteca).

Ecuación Logística

Primero, veamos como luce $\mu(x)$ con decrecimiento lineal respecto a la población x.

Como queremos que $\mu(0)=1$ y $\mu(1)=0$, la línea recta que conecta estos puntos es...



In [48]:

    
import numpy as np                      # Librería python numérico
import matplotlib.pyplot as plt         # Librería para graficar
%matplotlib inline

# Poner números en los ejes con tamaño personalizado
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size 
mpl.rcParams['ytick.labelsize'] = label_size 

# Definimos la tasa de crecimiento lineal
def mu(r,x):
    return r*(1-x)

r = 1                                   # Tasa máxima de crecimiento
N = 200                                 # Número de puntos para graficar
x = np.linspace(0,1.2,N)                # Vector de valores de población

# Gráfico de la función mu(x)
plt.figure()
plt.plot(x,mu(r,x),linewidth=2,label="$\mu(x)=r(1-x)$")
plt.plot(np.array([0, 1.2]),np.array([0, 0]),'--k')
plt.plot(np.array([0, 0]),np.array([0, 1]),'--k')
plt.legend(loc='best')
plt.grid()
plt.xlabel("$x$", fontsize=18)
plt.ylabel("$\mu(x)$", fontsize=18)
plt.show()

Entonces, con esta elección de $\mu(x)=r(1-x)$, obtenemos la llamada ecuación lógistica, publicada por Pierre Verhulst en 1838.

$$\frac{dx}{dt} = r\; x\; (1- x)$$

Solución a la ecuación diferencial

La ecuación diferencial inicial tiene solución analítica, $$ x(t) = \frac{1}{1+ (\frac{1}{x_{0}}- 1) e^{-rt}}.$$

Graficamos varias curvas de la solución analítica para $r = \left[-1, 1\right]$.



In [49]:

    
# Definimos la solución analítica
def logi_sol(t, x0, r):
    return 1/(1 + (1/x0 - 1) * np.exp(-r * t))

t = np.linspace(0,10)                                 # Vector de tiempo
x0 = 0.05                                          # Condicion inicial
# Gráficos para diferentes r
for r in np.arange(-1, 1,.1):
    plt.plot(t, logi_sol(t, x0, r))

plt.grid()
plt.xlabel('$t$', fontsize = 18)
plt.ylabel('$x$', fontsize = 18)
plt.show()

Como podemos ver, la solución a está ecuación en el continuo nos puede ganantizar la extinción o bien un crecimiento descomunal, dependiendo del valor asignado a $r$.

Numéricamente, ¿cómo resolveríamos esta ecuación?



In [50]:

    
from scipy.integrate import odeint             # Función de solución numérica de EDO

# Definicimos la ecuación logística para integrar con odeint
def poblacion(x, t, r):
    return mu(r,x)*x

r = 1                                          # Tasa de crecimiento máxima
x0 = 0.05                                      # Población inicial
tt = t                                         # Vector de tiempo de 0 a 10
xx = odeint(poblacion, x0, tt, (r,))           # Solución numérica

# Gráfico de la solución
plt.figure()
plt.plot(tt, xx)
plt.grid()
plt.xlabel('$t$', fontsize = 18)
plt.ylabel('$x$', fontsize = 18)
plt.show()

¿Qué tan buena es la aproximación de la solución numérica?

Hay ecuaciones diferenciales ordinarias no lineales para las cuales es imposible obtener la solución exacta. En estos casos, se evalúa una solución aproximada de forma numérica.

Para el caso anterior fue posible obtener la solución exacta, lo cual nos permite comparar ambas soluciones y evaluar qué tan buena es la aproximación que nos brinda la solución numérica.

Primero veamos esto gráficamente



In [58]:

    
x_numerica = xx[:,0]                        # Extraemos el vector de la solución numérica
x_exacta = logi_sol(t,x0,r)                 # Solución exacta con r=1

# Gráfico de la comparación
plt.figure()
plt.plot(t, x_exacta, '-b', lw = 4, label = 'solución exacta')
plt.plot(t, x_numerica, '*k', lw = 4, label = 'solución numérica')
plt.grid()
plt.legend(loc='best')
plt.xlabel('$t$', fontsize = 18)
plt.ylabel('$x$', fontsize = 18)
plt.show()

Gráficamente vemos que la solución numérica está cerca (coincide) con la solución exacta. Sin embargo, con esta gráfica no podemos visualizar qué tan cerca están una solución de la otra. ¿Qué tal si evaluamos el error?



In [59]:

    
error = x_exacta-x_numerica                 # Vector de error de aproximación

# Gráfico del error
plt.figure()
plt.plot(t, error, '-r', lw = 4, label = '$e=x_{exacta}-x_{numerica}$')
plt.grid()
plt.legend(loc='best')
plt.xlabel('$t$', fontsize = 18)
plt.ylabel('$error$', fontsize = 18)
plt.show()



In [60]:

    
# Gráfico de la comparación y del error
fig, (ax1, ax2) = plt.subplots(1, 2,figsize =(13,4.5))

ax1.plot(t, x_exacta, '-r', lw = 4, label = 'solución exacta')
ax1.plot(t, x_numerica, '*k', lw = 4, label = 'solución numérica')
ax1.grid()
ax1.legend(loc='best')
ax1.set_xlabel('$t$', fontsize = 18)
ax1.set_ylabel('$x$', fontsize = 18)

ax2.plot(t, error, '-r', lw = 4, label = '$e=x_{exacta}-x_{numerica}$')
ax2.grid()
ax2.legend(loc='best')
ax2.set_xlabel('$t$', fontsize = 18)
ax2.set_ylabel('$error$', fontsize = 18)

plt.show()

Entonces, cualitativamente ya vimos que la solución numérica es suficientemente buena. De todas maneras, es siempre bueno cuantificar qué tan buena es la aproximación. Varias formas:

Norma del error: tenemos el error de aproximación en ciertos puntos (especificados por el vector de tiempo). Este error es entonces un vector y le podemos tomar su norma 2

$$||e||_2=\sqrt{e[0]^2+\dots+e[n-1]^2}$$



In [52]:

    
np.linalg.norm(error)









    Out[52]:





2.7486108664240486e-07

Error cuadrático medio: otra forma de cuantificar es con el error cuadrático medio

$$e_{ms}=\frac{e[0]^2+\dots+e[n-1]^2}{n}$$



In [61]:

    
np.mean(np.sum(error**2))









    Out[61]:





7.5548616950243596e-14

Integral del error cuadrático: evalúa la acumulación de error cuadrático. Se puede evaluar cabo con la siguiente aproximación rectangular de la integral

$$e_{is}=\int_{0}^{t_f}e(t)^2\text{d}t\approx \left(e[0]^2+\dots+e[n-1]^2\right)h$$

donde $h$ es el tamaño de paso del vector de tiempo.



In [62]:

    
h = t[1]-t[0]
np.sum(error**2)*h









    Out[62]:





1.5418085091886449e-14

Comentarios del modelo logístico

El modelo no se debe tomar literalmente. Más bien se debe interpretar metefóricamente como que la población tiene una tendencia a crecer hasta su tope, o bien, desaparecer.

La ecuación logística fue probada en experimentos de laboratorio para colonias de bacterias en condiciones de clima constante, abastecimiento de comida y ausencia de predadores. Los experimentos mostraron que la ecuación predecía muy bien el comportamiento real.

Por otra parte, la predicción no resultó tan buena para moscas que se alimentan de frutas, escarabajos y otros organismos con ciclos de vida complejos. En esos casos se observaron fluctuaciones (oscilaciones) inmensas de la población.

Mapa logístico

La ecuación logística (curva de crecimiento logístico) es un modelo del crecimiento continuo en el tiempo. Una modificación de la ecuación continua a una ecuación de recurrencia discreta conocida como mapa logistico es muy usada.

Referencia:

Si reemplazamos la ecuación logísitica por la ecuación a diferencias:
$$x_{n+1} = r\; x_{n}(1- x_{n}),$$
donde $r$ es la razón de crecimiento máximo de la población y $x_{n}$ es la n-ésima iteración. Entonces, lo que tenemos que programar es la siguiente relación recursiva
$$x_{n+1}^{(r)} = f_r(x_n^{(r)}) = rx_n^{(r)}(1-x_n^{(r)})$$

El siguiente gif muestra las primeras 63 iteraciones de la anterior ecuación para diferentes valores de $r$ variando entre 2 y 4.

Tomado de https://upload.wikimedia.org/wikipedia/commons/1/1f/Logistic_map_animation.gif.

Note que:

Para $2<r<3$ el las soluciones se estabilizan en un valor de equilibrio.
Para $3<r<1+\sqrt{6}\approx 3.44949$ el las soluciones oscilan entre dos valores.
Para $3.44949<r<3.54409$ las soluciones oscilan entre cuatro valores.
Para $r>3.54409$ las soluciones exhiben un comportamiento caótico.

Caos: comportamiento determinista aperiódico muy sensible a las condiciones iniciales. Es decir, pequeñas variaciones en dichas condiciones iniciales pueden implicar grandes diferencias en el comportamiento futuro

¿Cómo podemos capturar este comportamiento en una sola gráfica?



In [18]:

    
# Definición de la función mapa logístico
def mapa_logistico(r, x):
    return r * x * (1 - x)

# Para mil valores de r entre 2.0 y 4.0
n = 1000
r = np.linspace(2.0, 4.0, n)

# Hacemos 1000 iteraciones y nos quedamos con las ultimas 100 (capturamos el comportamiento final)
iteraciones = 1000
ultimos = 100

# La misma condición inicial para todos los casos. 
x = 1e-5 * np.ones(n)

# Gráfico
plt.figure(figsize=(7, 5))
for i in np.arange(iteraciones):
    x = mapa_logistico(r, x)
    if i >= (iteraciones - ultimos):
        plt.plot(r, x, ',k', alpha=.2)
plt.xlim(np.min(r), np.max(r))
plt.ylim(-.1, 1.1)
plt.title("Diagrama de bifurcación", fontsize=20)
plt.xlabel('$r$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.show()



In [23]:

    
fig, (ax1, ax2) = plt.subplots(1, 2, sharex='col', sharey='row',figsize =(13,4.5))
r = np.linspace(.5, 4.0, n)
for i in np.arange(iteraciones):
    x = mapa_logistico(r, x)
    if i >= (iteraciones - ultimos):
        ax1.plot(r, x, '.k', alpha=1,  ms = .1)
r = np.linspace(2.5, 4.0, n)
for i in np.arange(iteraciones):
    x = mapa_logistico(r, x)
    if i >= (iteraciones - ultimos):
        ax2.plot(r, x, '.k', alpha=1, ms = .1)
ax1.set_xlim(.4, 4)
ax1.set_ylim(-.1, 1.1)
ax2.set_xlim(2.5, 4)
ax2.set_ylim(-.1, 1.1)
ax1.set_ylabel('$x$', fontsize = 20)
ax1.set_xlabel('$r$', fontsize = 20)
ax2.set_xlabel('$r$', fontsize = 20)
plt.show()



In [25]:

    
fig, (ax1, ax2) = plt.subplots(1, 2, sharex='col', sharey='row',figsize =(13,4.5))
r = np.linspace(.5, 4.0, n)
for i in np.arange(iteraciones):
    x = mapa_logistico(r, x)
    if i >= (iteraciones - ultimos):
        ax1.scatter(r, x,  s = .1, cmap= 'inferno', c = x, lw = 0)
r = np.linspace(2.5, 4.0, n)
for i in np.arange(iteraciones):
    x = mapa_logistico(r, x)
    if i >= (iteraciones - ultimos):
        ax2.scatter(r, x, s = .1, cmap = 'inferno', c = x, lw = 0)
ax1.set_xlim(.4, 4)
ax1.set_ylim(-.1, 1.1)
ax2.set_xlim(2.5, 4)
ax2.set_ylim(-.1, 1.1)
ax1.set_ylabel('$x$', fontsize = 20)
ax1.set_xlabel('$r$', fontsize = 20)
ax2.set_xlabel('$r$', fontsize = 20)
plt.show()

Esta última opción ya no es tan buena, porque son demasiados puntos y el uso de memoria aumenta.

Big data, insane!

Para esto veremos la siguiente clase como graficar lo anterior en mapas circulares...

preguntar por reposición de clase de viernes 22

Gráficando el mapeo logístico de forma circular

Se gráfica el seno y el cosenode un ángulo que va aumentando gradualmente, multiplicado por el valor de x, el cual también va cambiando con cada iteración. Vamos a necesitar millones de puntos.

\begin{align} X &= x \,R \cos(\theta)\\ Y &= x \,R \sin(\theta) \end{align}

Vamos a considerar un ángulo inicial $\theta_0 = 0 $, $R =1$ y un incremento $\delta \theta = 0.2$.

Los valores de $r$ utilizados para contruir las figuras de arriba son los siguientes (de izquierda a derecha e iniciando con el panel superior):

$$r = 0.9,\, 2.5,\, 3.2,\, 3.46,\, 3.57,\, 3.59,\, 3.64,\, 3.83,\, 3.9,\, 3.99$$

A continuación se describe la construcción de las figuras anteriores.

https://github.com/bokeh/datashader

conda install -c bokeh datashader



In [480]:

    
import pandas as pd
import datashader as ds
from datashader import transfer_functions as tf
from datashader.colors import Greys9, inferno, viridis
from datashader.utils import export_image
from functools import partial
background = "black"



In [481]:

    
img_map = partial(export_image, export_path="circular_maps", background=background)



In [482]:

    
def circular_map(r1, theta, dtheta, radius, iterar):
    r = np.array([r1])
    x = np.array([1e-5])
    x_list, y_list = [], []
    for i in range(iterar):
        x = mapa_logistico(r, x)
        X = x * radius * np.cos(theta)
        Y = x * radius * np.sin(theta)
        x_list.append(X[0])
        y_list.append(Y[0])
        theta = theta + dtheta
    return x_list, y_list

Los valores de $r$ los vamos a elegir en los puntos de transición en el diagrama de bifurcación.



In [483]:

    
x_list, y_list = circular_map(.9, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img = tf.dynspread(img, threshold=1, max_px= 5)
img_map(img,"map_dead")









    Out[483]:



In [484]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img = tf.dynspread(img, threshold=1, max_px= 5)
img_map(img,"map_dead_1")









    Out[484]:



In [485]:

    
x_list, y_list = circular_map(2.5, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
#img = tf.dynspread(img, threshold=1, max_px= 1)
img_map(img,"map_alive")









    Out[485]:



In [486]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
#img = tf.dynspread(img, threshold=1, max_px= 1)
img_map(img,"map_alive_1")









    Out[486]:



In [487]:

    
x_list, y_list = circular_map(3.2, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive2")









    Out[487]:



In [488]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive2_1")









    Out[488]:



In [489]:

    
x_list, y_list = circular_map(3.46, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive4")









    Out[489]:



In [490]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive4_1")









    Out[490]:



In [491]:

    
x_list, y_list = circular_map(3.57, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos1")









    Out[491]:



In [492]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos1_1")









    Out[492]:



In [493]:

    
x_list, y_list = circular_map(3.59, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos2")









    Out[493]:



In [494]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos2_1")









    Out[494]:



In [495]:

    
x_list, y_list = circular_map(3.64, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos3")









    Out[495]:



In [496]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos3_1")









    Out[496]:



In [497]:

    
x_list, y_list = circular_map(3.83, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1,1), y_range=(-1,1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_stable")









    Out[497]:



In [498]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_stable_1")









    Out[498]:



In [499]:

    
x_list, y_list = circular_map(3.9, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1.1,1.1), y_range=(-1.1,1.1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos4")









    Out[499]:



In [500]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos4_1")









    Out[500]:



In [501]:

    
x_list, y_list = circular_map(3.99, 0, 0.2, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1.1,1.1), y_range=(-1.1,1.1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos5")









    Out[501]:



In [ ]:



In [502]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos5_1")









    Out[502]:



In [510]:

    
x_list, y_list = circular_map(3.9, 0, 0.1, 1, 3000000)
d = {'r': x_list, 'x': y_list}
df1 = pd.DataFrame(d)
cvs1 = ds.Canvas(plot_width=300, plot_height=300, x_range=(-1.1,1.1), y_range=(-1.1,1.1))
agg1 = cvs1.points(df1, 'r', 'x')
img = tf.shade(agg1, cmap = inferno, how='eq_hist')
img_map(img,"map_alive_chaos6")









    Out[510]:



In [511]:

    
img = tf.shade(agg1, cmap = viridis, how='eq_hist')
img_map(img,"map_alive_chaos6_1")









    Out[511]: