1. Ejercicio - Modelo depredador-presa
Referencia: Chapra 28.2 p831 pdf855, Rodríguez 9.2.1 p263
wikipedia.org - Ecuaciones Lotka Volterra ,wikipedia.org - Lotka Volterra equations
Los modelos depredador-presa y caos descritos en las ecuaciones Lotka-Volterra tienen como variables:
t = tiempo de observación
x = número de presas
y = número de depredadores
Los coeficientes en la ecuación se describen como:
a = razón de crecimiento de la presa, (0.5)
b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
c = razón de muerte del depredador (0.35)
d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)
Para el ejercicio, considere como puntos iniciales en la observación de las especies:
t=0, x=2, y=1, h=0.5
Los términos que multiplican xy hacen que las ecuaciones sean no lineales.
Observe que la variable tiempo no se encuentra en las expresiones f y g, h se aplica a la variable tiempo.
Referencia: Can You Spot This Leopard Before Its Prey Does? | Nat Geo Wild. Nat Geo Animals. 1 octubre-2017
2. Desarrollo analítico - paso a paso
Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:
f(t,x,y) = 0.5 x - 0.7 xy g(t,x,y) = -0.35y + 0.35xyLas expresiones se adaptan al método de Runge-Kutta para primeras derivadas por cada variable de población. Se deben usar de forma simultánea para cada tiempo t.
K1x = h f(t,x,y) = 0.5 \Big( 0.5 x - 0.7 xy \Big) K1y = h g(t,x,y) = 0.5 \Big(-0.35y + 0.35xy \Big)..
K2x = h f(t+h,x+K1x,y+K1y) = 0.5 \Big( 0.5 (x+K1x) - 0.7 (x+K1x)(y+K1y) \Big) K2y = h g(t+h,x+K1x,y+K1y) = 0.5 \Big(-0.35(y+K1y) + 0.35(x+K1x)(y+K1y) \Big)..
x[i+1] = x[i] + \frac{K1x+K2x}{2} y[i+1] = y[i] + \frac{K1y+K2y}{2} t[i+1] = t[i] + hcon lo que se puede aplicar al ejercicio en cada iteración,dadas las condiciones iniciales.
itera = 0
t = 0, x = 2, y = 1, h = 0.5
K1x = 0.5 f(0,2,1) = 0.5 \Big( 0.5 (2) - 0.7 (2)(1) \Big) = -0.2 K1y = 0.5 g(0,2,1) = 0.5 \Big(-0.35(1) + 0.35(2)(1) \Big) =0.175..
K2x = 0.5 f(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big( 0.5 (2+(-0.2)) - 0.7 (2+(-0.2))(1+0.175) \Big) = -0.29025 K2y = 0.5 g(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big(-0.35(1+0.175) + 0.35(2+(-0.2))(1+0.175) \Big) = 0.1645..
x[1] = x[0] + \frac{K1x+K2x}{2} = 2 + \frac{-0.2+(-0.29025)}{2} = 1.7548 y[1] = y[0] + \frac{K1y+K2y}{2} = 1 + \frac{0.175+0.1645}{2}= 1.1697 t[1] = t[0] + h = 0 +0.5 = 0.5
itera = 1
t = 0.5, x = 1.7548, y = 1.1697, h = 0.5
K1x = 0.5 \Big( 0.5 (0,1.7548) - 0.7 (0,1.7548)(1.1697) \Big) = -0.2797 K1y = 0.5 \Big(-0.35(1.1697) + 0.35(1.7548)(1.1697) \Big) = 0.1545..
K2x = 0.5 \Big( 0.5 (1.7548+(-0.2797)) - 0.7 (1.7548+(-0.2797))(1.1697+0.1545) \Big) =-0.3149 K2y = 0.5 \Big(-0.35(1.1697+0.1545) + 0.35(1.7548+(-0.2797))(1.1697+0.1545) \Big) = 0.1645..
x[2] = 1.7548 + \frac{-0.2797+(-0.3149)}{2} = 1.4575 y[2] = 1.1697 + \frac{0.1545+0.1645}{2} = 1.3020 t[2] = t[0] + h = 0.5 +0.5 = 1
itera=2
t = 1, x = 1.4575, y = 1.3020, h = 0.5
continuar como tarea ...
3. Algoritmo en Python para Sistemas EDO con Runge-Kutta
Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en
Al ejecutar el algoritmo se obtienen los siguientes resultados:
Runge-Kutta Segunda derivada
i [ xi, yi, zi ]
[ K1y, K1z, K2y, K2z ]
0 [0. 2. 1.]
[0. 0. 0. 0.]
1 [0.5 1.754875 1.16975 ]
[-0.2 0.175 -0.29025 0.1645 ]
2 [1. 1.457533 1.302069]
[-0.279749 0.154528 -0.314935 0.11011 ]
3 [1.5 1.167405 1.373599]
[-0.29985 0.104254 -0.280406 0.038807]
4 [2. 0.922773 1.381103]
[-0.26939 0.040241 -0.219874 -0.025233]
5 [2.5 0.734853 1.33689 ]
[-0.215362 -0.018665 -0.160478 -0.069761]
6 [3. 0.598406 1.258434]
[-0.160133 -0.062033 -0.11276 -0.09488 ]
... Los resultados de la tabla se muestran parcialmente, pues se usaron mas de 100 iteraciones.
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.
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
# 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
# Algoritmo como función
def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
vertabla=False, precision=6):
''' solucion a EDO d2y/dx2 con Runge-Kutta 2do Orden,
f(x,y,z) = z #= y'
g(x,y,z) = expresion d2y/dx2 con z=y'
tambien es solucion a sistemas edo f() y g()
x0,y0,z0 son valores iniciales, h es tamano de paso,
muestras es la cantidad de puntos a calcular.
'''
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,3+4),dtype=float)
# incluye el punto [x0,y0,z0,K1y,K1z,K2y,K2z]
tabla[0] = [x0,y0,z0,0,0,0,0]
xi = x0 # valores iniciales
yi = y0
zi = z0
for i in range(1,tamano,1):
K1y = h * f(xi,yi,zi)
K1z = h * g(xi,yi,zi)
K2y = h * f(xi+h, yi + K1y, zi + K1z)
K2z = h * g(xi+h, yi + K1y, zi + K1z)
yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
xi = xi + h
tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
if vertabla==True:
np.set_printoptions(precision)
print('EDO f,g con Runge-Kutta 2 Orden')
print('i ','[ xi, yi, zi',']')
print(' [ K1y, K1z, K2y, K2z ]')
for i in range(0,tamano,1):
txt = ' '
if i>=10:
txt = ' '
print(str(i),tabla[i,0:3])
print(txt,tabla[i,3:])
return(tabla)
# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,
muestras,vertabla=True)
# SALIDA
print('Sistemas EDO: Modelo presa-predador')
##print('i ','[ xi, yi, zi',']')
##print(' [ K1y, K1z, K2y, K2z ]')
##for i in range(0,tamano,1):
## txt = ' '
## if i>=10:
## txt = ' '
## print(str(i),tabla[i,0:3])
## print(txt,tabla[i,3:])
4. Gráfica en Python
Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:
# GRAFICA ---------------------
import matplotlib.pyplot as plt
titulo = 'Sistemas EDO Runge-Kutta 2ord - Modelo predador-presa'
i = muestras # iteración en gráfica
titulo = titulo+', i='+str(i)
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
K1y = tabla[:,3]
K1z = tabla[:,4]
K2y = tabla[:,5]
K2z = tabla[:,6]
plt.subplot(211)
plt.plot(xi[0],yi[0],'o',
color='red', label ='[t0,y0]')
if i<20: # evita muchos puntos
plt.plot(xi[1:i+2],yi[1:i+2],'o',
color='green', label ='[t[i],y[i]]')
plt.plot(xi[0:i+2],yi[0:i+2],
color='blue',label='y(t)')
if i<muestras: # gráfica para una iteración
plt.plot(xi[i+1],yi[i+1],'o',color='orange',
label ='[x[i+1],y[i+1]]')
plt.plot(xi[i:i+3],yi[i:i+3],'.',color='gray')
plt.plot(xi[i:i+2],[yi[i],yi[i]], color='orange',
label='h',linestyle='dashed')
plt.plot([xi[i+1]-0.02*h,xi[i+1]-0.02*h],
[yi[i],yi[i]+K1y[i+1]],
color='green',label='K1y',linestyle='dashed')
plt.plot([xi[i+1]+0.02*h,xi[i+1]+0.02*h],
[yi[i],yi[i]+K2y[i+1]],
color='magenta',label='K2y',linestyle='dashed')
plt.plot([xi[i+1]-0.02*h,xi[i+1]+0.02*h],
[yi[i]+K1y[i+1],yi[i]+K2y[i+1]],
color='magenta')
if np.min(yi[0:i+1])<0: # linea 0
plt.axhline(0, color='red')
plt.ylabel('y = presa')
plt.title(titulo)
plt.legend()
plt.grid()
plt.tight_layout()
plt.subplot(212)
plt.plot(xi[0],zi[0],'o',
color='red', label ='[t0,z0]')
if i<20: # evita muchos puntos
plt.plot(xi[1:i+2],zi[1:i+2],'o',
color='green', label ='[t[i],v[i]]')
plt.plot(xi[0:i+2],zi[0:i+2],
color='green',label='z(t)')
if i<muestras: # gráfica para una iteración
plt.plot(xi[i+1],zi[i+1],'o',color='orange',
label ='[t[i+1],z[i+1]]')
plt.plot(xi[i:i+3],zi[i:i+3],'.',color='gray')
plt.plot(xi[i:i+2],[zi[i],zi[i]], color='orange',
label='h',linestyle='dashed')
plt.plot([xi[i+1]-0.02*h,xi[i+1]-0.02*h],
[zi[i],zi[i]+K1z[i+1]],
color='green',label='K1z',linestyle='dashed')
plt.plot([xi[i+1]+0.02*h,xi[i+1]+0.02*h],
[zi[i],zi[i]+K2z[i+1]],
color='magenta',label='K2z',linestyle='dashed')
plt.plot([xi[i+1]-0.02*h,xi[i+1]+0.02*h],
[zi[i]+K1z[i+1],zi[i]+K2z[i+1]],
color='magenta')
plt.xlabel('x = tiempo')
plt.ylabel('z = predador')
plt.legend()
plt.grid()
plt.tight_layout()
#plt.show() #comentar para la siguiente gráfica
# gráfica yi vs zi
titulo = 'Sistemas EDO Runge-Kutta 2ord - presa vs predador'
titulo = titulo+', i='+str(i)
fig_yz, graf3 = plt.subplots()
graf3.plot(yi,zi)
graf3.set_title(titulo)
graf3.set_xlabel('presa')
graf3.set_ylabel('predador')
graf3.grid()
plt.show() #comentar para la siguiente gráfica
El Problema de los Tres Cuerpos, una Visualización del CAOS del Cosmos. Mates Mike, 8 febrero 2024.