AMIGOS PARA SIEMPRE: Apuntes de cálculo

ecuaciones diferenciales

ecuaciones de Lotka-Volterra, también conocidas como ecuaciones predador-presa o presa-predador, son un par de ecuaciones diferenciales de primer orden no lineales que se usan para describir dinámicas de sistemas biológicos en el que dos especies interactúan, una como presa y otra como depredador. Las ecuaciones fueron propuestas de forma independiente por Alfred J. Lotka en 1925 y Vito Volterra en 1926. Tales ecuaciones se definen como:

$\frac{dx}{dt} = x(\alpha - \beta y)$

$\frac{dy}{dt} = - y(\gamma + \delta x)$

donde:

y es el número de algún predador (por ejemplo, un lobo);
x es el número de sus presas (por ejemplo, conejos);
dy/dt y dx/dt representa el crecimiento de las dos poblaciones en el tiempo;
t representa el tiempo; y
α, β, γ y δ son parámetros (positivos) que representan las interacciones de las dos especies.

Explicación de las ecuaciones

Usando las series de Taylor se obtiene una solución lineal a las ecuaciones:

$f(x,y) = A_0 - A_1 x - A_2 y$

$g(x,y) = B_0 + B_1 x - B_2 y$ .

Con estos coeficientes se puede estudiar los modelos de competición, enfermedad y mutualismo (biología) en un ecosistema.

Presa

$\frac{dx}{dt} = \alpha x - \beta x y$

Se asume que las presas tienen suministro de comida ilimitado por tiempo definido, y se reproducen exponencialmente a menos que exista algún predador. Estecrecimiento exponencial está representado en la ecuación por el término αx. El término de la ecuación βxy viene a representar el encuentro de las dos especies y su interacción. Si x o y son cero no existe interacción.

Se puede interpretar la ecuación como el cambio del número de presas viene dado por su propio crecimiento menos la tasa de encuentros con predadores.

Depredador

$\frac{dy}{dt} = \delta xy - \gamma y$

En esta ecuación, δxy representa el crecimiento de los depredadores (fíjese en la similitud con la ecuación para las presas, pero en este caso para el crecimiento de los depredadores es necesario usar la razón a la que se consumen las presas, x). γy representa la muerte natural de los depredadores de forma exponencial; a más depredadores es necesario que el número de víctimas o presa aumente para mantener la población.

Se puede interpretar la ecuación como el crecimiento de los depredadores por la caza de presas menos la muerte natural de éstos.

Las ecuaciones de Lotka-Volterra son un modelo biomatemático que pretende responder a estas cuestiones prediciendo la dinámica de las poblaciones de presa y depredador bajo una serie de hipótesis:

El ecosistema está aislado: no hay migración, no hay otras especies presentes, no hay plagas...
La población de presas en ausencia de depredadores crece de manera exponencial: la velocidad de reproducción es proporcional al número de individuos. Las presas sólo mueren cuando son cazadas por el depredador.
La población de depredadores en ausencia de presas decrece de manera exponencial.
La población de depredadores afecta a la de presas haciéndola decrecer de forma proporcional al número de presas y depredadores (esto es como decir de forma proporcional al número de posibles encuentros entre presa y depredador).
La población de presas afecta a la de depredadores también de manera proporcional al número de encuentros, pero con distinta constante de proporcionalidad (dependerá de cuanto sacien su hambre los depredadores al encontrar una presa).

Se trata de un sistema de dos ecuaciones diferenciales de primer orden, acopladas, autónomas y no lineales:

dxdt=αx−βxy

dydt=−γy+δyx

donde

x es el número de presas (cebras en nuestro caso) e

y es el número de depredadores (leones). Los parámetros son constantes positivas que representan:

α: tasa de crecimiento de las presas.
β: éxito en la caza del depredador.
γ: tasa de decrecimiento de los depredadores.
δ: éxito en la caza y cuánto alimenta cazar una presa al depredador.

Resolución

Resolveremos este sistema en Python usando la función odeint de scipy.integrate(puedes ver cómo funciona en el artículo El salto de Felix Baumgartner en Python) y representaremos los resultados con matplotlib. Para esto usaremos: Python 3.4, numpy 1.9.0, matplotlib 1.4.0, scipy 0.14.0.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def df_dt(x, t, a, b, c, d):
    """Función del sistema en forma canónica"""
    dx = a * x[0] - b * x[0] * x[1]
    dy = - c * x[1] + d * x[0] * x[1]
    return np.array([dx, dy])

# Parámetros
a = 0.1
b = 0.02
c = 0.3
d = 0.01

# Condiciones iniciales
x0 = 40   # Presas
y0 = 9    # Depredadores
conds_iniciales = np.array([x0, y0])

# Condiciones para integración
tf = 200
N = 800
t = np.linspace(0, tf, N)

solucion = odeint(df_dt, conds_iniciales, t, args=(a, b, c, d))

Ahora que ya tenemos la solución, sólo tenemos que pintarla para ver cómo han ido las cosas en nuestra sabana virtual.

plt.plot(t, solucion[:, 0], label='presa')
plt.plot(t, solucion[:, 1], label='depredador')

Vemos que se trata de una solución periódica en la que, como decíamos al principio, un aumento en la población de cebras va seguido del aumento del número de leones. Un gran número de depredadores merma la población de presas y a los pobres leones les toca pasar hambre una temporada. Otra forma interesante de visualizar estos datos es ver el número de presas en función del número de depredadores, en lugar de a lo largo del tiempo, es decir, podemos visualizar su mapa de fases. Podemos pintar también el campo de direcciones de nuestras ecuaciones usando la función quiver. El tamaño de las flechas se ha normalizado para que todas tengan la misma longitud y se ha usado un colormap para representar el módulo.

x_max = np.max(solucion[:,0]) * 1.05
y_max = np.max(solucion[:,1]) * 1.05

x = np.linspace(0, x_max, 25)
y = np.linspace(0, y_max, 25)

xx, yy = np.meshgrid(x, y)
uu, vv = df_dt((xx, yy), 0, a, b, c, d)
norm = np.sqrt(uu**2 + vv**2)
uu = uu / norm
vv = vv / norm

plt.quiver(xx, yy, uu, vv, norm, cmap=plt.cm.gray)
plt.plot(solucion[:, 0], solucion[:, 1])

Si nos fijamos en la línea azul, la coordenada

x en cada punto indica el número de presas y la coordenada

y el número de depredadores. La evolución a lo largo del tiempo que hemos representado antes, se obtiene al recorrer esta curva en sentido antihorario. Podemos ver también como el campo de direcciones nos señala la tendencia del sistema en cada situación. Por ejemplo, una flecha que apunta hacia arriba a la derecha indica que con ese número de cebras y leones en nuestra sabana, la tendencia será que aumenten ambos.

Llegados a este punto podemos preguntarnos qué habría ocurrido si el número inicial de cebras y leones hubiese sido otro. Como ya sabemos integrar ecuaciones diferenciales, bastaría con cambiar nuestra x0 e y0 y repetir el proceso (incluso podríamos hacer un widget interactivo). Sin embargo, se puede demostrar que a lo largo de las líneas del mapa de fases, como la que hemos pintado antes, se conserva la cantidad:

C=αlny−βy+γlnx−δx

Por tanto, pintando un contour de esta cantidad podemos obtener la solución para distintos valores iniciales del problema.

def C(x, y, a, b, c, d):
    return a * np.log(y) - b * y + c * np.log(x) - d * x

fig, ax = plt.subplots(1,2)

ax[0].plot(solucion[:, 0], solucion[:, 1], lw=2, alpha=0.8)
ax[0].scatter(c/d, a/b)
levels = (0.5, 0.6, 0.7, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.775, 0.78, 0.781)
ax[0].contour(xx, yy, constant, levels, colors='blue', alpha=0.3)

ax[1].plot(t, solucion[:, 0], label='presa')
ax[1].plot(t, solucion[:, 1], label='depredador')

Vemos que estas curvas se van haciendo cada vez más y más pequeñas, hasta que, en nuestro caso, colapsarían en un punto en torno a

(30,5). Se trata de un punto de equilibrio o punto crítico; si el sistema lo alcanzase, no evolucionaría y el número de cebras y leones sería constante en el tiempo. El otro punto crítico de nuestro sistema es el

(0,0). Analizándolos matemáticamente se obtiene que:

El punto crítico situado en (0,0) es un punto de silla. Al tratarse de un punto de equilibrio inestable la extinción de cualquiera de las dos especies en el modelo sólo puede conseguirse imponiendo la condición inicial nula.
El punto crítico situado en (γ/δ,α/β) es un centro (en este caso los autovalores de la matriz del sistema linealizado son ambos imaginarios puros, por lo que a priori no se conoce su estabilidad).

La ecuación de Poisson-Boltzmann es una ecuación diferencial que describe interacciones electrostáticas entre moléculas ensoluciones iónicas. La ecuación puede ser usada como fundamento matemático del modelo de la Doble Capa Eléctrica Interfacial deGouy-Chapman, propuesta inicialmente por L. G. Gouy en 1910 y completada por Chapman en 1913.

La ecuación es importante en los campos de la dinámica molecular y la biofísica, porque puede usarse en el modelado de disoluciones continuas, como aproximación de los efectos de los disolventes en estructuras de proteínas, ADN, ARN, y otras moléculas en disoluciones de distinta fuerza iónica. Algunas veces la ecuación Poisson–Boltzmann resulta difícil de resolver para sistemas complejos, problema que se está solucionando con el desarrollo del análisis numérico por ordenador. La ecuación puede escribirse como:

$\vec{\nabla}\cdot\left[\epsilon(\vec{r})\vec{\nabla}\Psi(\vec{r})\right] = -4\pi\rho^{f}(\vec{r}) - 4\pi\sum_{i}c_{i}^{\infty}z_{i} q \lambda(\vec{r}) \cdot \exp \left[{\frac{-z_{i}q\Psi(\vec{r})}{k_B T}}\right]$

Donde:

\vec{\nabla}\cdot

es el operador divergencia,

\epsilon(\vec{r})

representa la posicion-dependencia dieléctrica,

\vec{\nabla} \Psi(\vec{r})

representa el gradiente del potencial electrostático,

\rho^{f}(\vec{r})

representa la densidad de carga del soluto,

c_{i}^{\infty}

representa la concentración del ion i a una distancia de infinito desde el soluto,

z_{i}

es la carga del ion, q es la carga del protón,

k_B

es la constante de Boltzmann,

T es la temperatura, y :

\lambda(\vec{r})

es un factor para la acesibilidad de la posicion-dependencia de la posición r hasta los iones en la disolución. Si el potencial no es grande, la ecuación se puede linealizar para así poder ser resuelta más fácilmente, llevando a la ecuación Debye–Hückel.

ecuación de Scheil-Gulliver (también llamada ecuación de Scheil o ecuación de congelación normal¹ ) describe redistribución del soluto en lasolidificación de una aleación. Este enfoque se aproxima a la solidificación en estado de no-equilibrio, suponiendo un equilibrio termodinámico local del frente de avance de la solidificación en la interfase sólido-líquido. Esto permite el uso de diagrama de fases en equilibrio para el análisis de la solidificación.

A diferencia de la solidificación en equilibrio, el soluto no se retrodifunde de nuevo en el sólido y es rechazado por completo en el líquido. La mezcla completa del soluto en el líquido también se asume como resultado de la convección y/o de la agitación.

Derivación

La siguiente figura muestra la distribución del soluto en una solidificación en estado de no-equilibrio donde no haydifusión del soluto en el sólido ni mezcla completa del soluto con el líquido.

$\ D_S = 0$ y $\ D_L = \infty$ .

Se supone que existe equilibrio en la interfase, lo que permite el uso de un diagrama de fases en equilibrio.

Las zonas rayadas en la figura representan la cantidad de soluto en el sólido y el líquido. Teniendo en cuenta que la cantidad total de soluto en el sistema debe conservarse, las áreas se igualan de la siguiente manera: