6.3 Sistemas EDO. modelo depredador-presa

Referencia: Chapra 28.2 p831 pdf855, Rodriguez 9.2.1 p263
https://es.wikipedia.org/wiki/Ecuaciones_Lotka%E2%80%93Volterra
https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

Modelos depredador-presa y caos. Ecuaciones Lotka-Volterra. En el sistema de ecuaciones:

\frac{dx}{dt} = ax - bxy \frac{dy}{dt} = -cy + dxy

variables
x
= número de presas
y = número de depredadores
t = tiempo de observación
coeficientes
a = razón de crecimiento de la presa, (0.5)
c = razón de muerte del depredador (0.35)
b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)

Los términos que multiplican xy hacen que las ecuaciones sean no lineales.

Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:

f =  a x - b x y
g = -c y + d x y

con valores iniciales: t = 0, x = 2, y = 1, h = 0.5

Observe que la variable tiempo no se encuentra en las expresiones f y g, h se aplica a tiempo.


Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en

Runge-Kutta d2y/dx2

Al ejecutar el algoritmo se obtienen los siguientes resultados:

 [ ti, xi, yi] 
[[  0.         2.         1.      ]
 [  0.5        1.754875   1.16975 ]
 [  1.         1.457533   1.302069]
 [  1.5        1.167405   1.373599]
 [  2.         0.922773   1.381103]
 [  2.5        0.734853   1.33689 ]
 [  3.         0.598406   1.258434]
 ...
 [ 49.         1.11309    1.389894]
 [ 49.5        0.876914   1.38503 ]
 [ 50.         0.698717   1.331107]
 [ 50.5        0.570884   1.246132]]
>>> 

Los resultados se pueden observar de diferentes formas:

a) Cada variable xi, yi versus ti, es decir cantidad de animales de cada especie durante el tiempo de observación

b) Independiente de la unidad de tiempo, xi vs yi, muestra la relación entre la cantidad de presas y predadores. Relación que es cíclica y da la forma a la gráfica.


Algoritmo con Python

Las instrucciones del algoritmo en Python usadas en el problema son:

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,3),dtype=float)
    tabla[0] = [t0,x0,y0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi]
    tabla = np.array(tabla)
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35

# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y

# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1

# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi]')
print(tabla)

Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:

# Grafica tiempos vs población
import matplotlib.pyplot as plt

plt.plot(ti,xi, label='xi presa')
plt.plot(ti,yi, label='yi predador')

plt.title('Modelo predador-presa')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('Modelo presa-predador [xi,yi]')
plt.xlabel('x presa')
plt.ylabel('y predador')
plt.grid()
plt.show()

Tarea: Añadir la animación de la gráfica usando la variable tiempo para mostrar un punto que marque el par ordenado[x,y] al variar t.