1. Ejercicio - Modelo predador-presa
Referencia: Chapra 28.2 p831 pdf855, Rodríguez 9.2.1 p263, Stewart 9.6 p627
Los modelos predador-presa y caos descritos en las ecuaciones Lotka-Volterra tienen como variables:
x = número de presas
y = número de predadores
t = tiempo de observación
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 predador-presa sobre la muerte de la presa (0.7)
c = razón de muerte del predador (0.35)
d = efecto de la interacción predador-presa sobre el crecimiento del predador, (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 se empieza a llenar la primera fila de la tabla.
| i | ti | xi | yi |
|---|---|---|---|
| 0 | 0 | 2 | 1 |
| ... |
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
| i | ti | xi | yi |
|---|---|---|---|
| 0 | 0 | 2 | 1 |
| 1 | 0.5 | 1.7548 | 1.1697 |
itera = 1
t = 0.5, x = 1.7548, y = 1.1697, h = 0.5
K1x = 0.5 \Big( 0.5 (1.7548) - 0.7 (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.1101..
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[1] + h = 0.5 +0.5 = 1
| i | ti | xi | yi |
|---|---|---|---|
| 0 | 0 | 2 | 1 |
| 1 | 0.5 | 1.7548 | 1.1697 |
| 2 | 1 | 1.4575 | 1.3020 |
itera=2
t = 1, x = 1.4575, y = 1.3020, h = 0.5
K1x = 0.5 \Big( 0.5 (1.4575) - 0.7 (1.4575)(1.3020) \Big) = -0.2998 K1y = 0.5 \Big(-0.35(1.3020) + 0.35(1.4575)(1.3020) \Big) = 0.1042 K2x = 0.5 \Big( 0.5 (1.4575+(-0.2998)) - 0.7 (1.4575+(-0.2998))(1.3020+0.1042) \Big) =-0.2804 K2y = 0.5 \Big(-0.35(1.3020+0.1042) + 0.35(1.4575+(-0.2998))(1.3020+0.1042) \Big) = 0.0388 x[3] = 1.4575 + \frac{-0.2998+(-0.2804)}{2} = 1.1674 y[3] = 1.3020 + \frac{0.1042+0.0388}{2} = 1.3735 t[3] = t[2] + h = 1 +0.5 = 1.5
| i | ti | xi | yi |
|---|---|---|---|
| 0 | 0 | 2 | 1 |
| 1 | 0.5 | 1.7548 | 1.1697 |
| 2 | 1 | 1.4575 | 1.3020 |
| 3 | 1.5 | 1.1674 | 1.3735 |
y así se continúa con las iteraciones hasta completar el intervalo que se requiere observar.
Se procede con la elaboración del algoritmo de Runge-Kutta para sistemas EDO.
3. Algoritmo en Python para Sistemas EDO con Runge-Kutta
Planteamiento que se ingresan al algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en
EDO d2y/dx2 Runge-Kutta con Python
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 se definen en el planteamiento analítico del ejercicio mediante las funciones f y g con condiciones iniciales. Lo planteado se incluye en el bloque de ingreso.:
# 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 RK2')
##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:])
Al optar por mostrar los resultados parciales dentro de la función rungekutta2_fg() mediante el parámetro vertabla, el bloque de salida ya no se usa para mostrar los resultados y se lo comenta con '#'.
4. Gráfica en Python
En el caso de la gráfica el primer bloque se usa para definir los títulos y etiquetas a usar en las gráficas.
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'
etiq_x = 'x = tiempo'
etiq_y = 'y = presa'
etiq_z = 'z = predador'
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(etiq_y)
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(etiq_x)
plt.ylabel(etiq_z)
plt.legend()
plt.grid()
plt.tight_layout()
#plt.show() #comentar para la siguiente gráfica
# gráfica yi vs zi
fig_yz, graf3 = plt.subplots()
graf3.plot(yi,zi)
graf3.set_title(titulo)
graf3.set_xlabel(etiq_y)
graf3.set_ylabel(etiq_z)
graf3.grid()
plt.show() #comentar para la siguiente gráfica
Un ejercicio por revisar sobre sistemas EDO:
El Problema de los Tres Cuerpos, una Visualización del CAOS del Cosmos. Mates Mike, 8 febrero 2024.