s3Eva_2024PAOII_T3 EDO efecto Allee en poblaciones pequeñas

Ejercicio: 3Eva_2024PAOII_T3 EDO efecto Allee en poblaciones pequeñas

literal a

Dada la ecuación diferencial ordinaria, se reemplazan los valores de las constantes:

$\frac{dx}{dt} = rx \Big(\frac{x}{K}-1 \Big) \Big(1-\frac{x}{A} \Big)$ $\frac{dx}{dt} = 0.7 x \Big(\frac{x}{10}-1 \Big) \Big(1-\frac{x}{50} \Big)$

siendo h = 0.2

$K_1 = 0.2 f(t_i,x_i) = 0.2\Big(0.7 x \Big(\frac{x}{10}-1 \Big) \Big(1-\frac{x}{50} \Big) \Big)$ $K_2 = 0.2 f(t_i+0.2, x_i + K_1)$ $= 0.2\Big(0.7(x+K_1) \Big(\frac{x+K_1}{10}-1 \Big) \Big(1-\frac{x+K_1}{50} \Big) \Big)$ $x_{i+1} = x_i + \frac{K_1+K_2}{2}$ $t_{i+1} = t_i + 0.2$

literal b

Para el desarrollo de las iteraciones se requieren valores iniciales. Para la variable independiente tiempo podría usar t = 0.

Para la variable dependiente población, según la descripción se encontraría entre el intervalo [10,50]. Si x₀<10 la población se extingue. Si la variable x₀>50 se encuentra saturada la capacidad del medio.

Por lo que se propone usar un valor mayor que el mínimo, por ejemplo x₀=11 y otro valor que se encuentre en el intervalo.

itera = 0 , t₀ = 0, x₀ = 11

$K_1 = 0.2\Big(0.7 (11) \Big(\frac{11}{10}-1 \Big) \Big(1-\frac{11}{50} \Big) \Big) = 0.1201$ $K_2 = 0.2\Big(0.7(11+0.1201) \Big(\frac{11+0.1201}{10}-1 \Big) \Big(1-\frac{11+0.1201}{50} \Big) \Big)$ $K_2= 0.1355$ $x_{i+1} = 11 + \frac{0.1201+0.1355}{2} = 11.1278$ $t_{i+1} = 0 + 0.2 = 0.2$

itera = 1 , t₀ = 0.2, x₀ = 11.1278

$K_1 = 0.2\Big(0.7 (11.1278) \Big(\frac{11.1278}{10}-1 \Big) \Big(1-\frac{11.1278}{50} \Big) \Big) = 0.1366$ $K_2 = 0.2\Big(0.7(11.1278+0.1366) \Big(\frac{11.1278+K_1}{10}-1 \Big) \Big(1-\frac{11.1278+0.1366}{50} \Big) \Big)$ $K_2= 0.1544$ $x_{i+1} = 11.1278 + \frac{0.1366+0.1544}{2} =11.2734$ $t_{i+1} = 0.2 + 0.2 = 0.4$

itera = 2 , t₀ = 0.4, x₀ = 11.2734

$K_1 = 0.2\Big(0.7 (11.2734) \Big(\frac{11.2734}{10}-1 \Big) \Big(1-\frac{11.2734}{50} \Big) \Big) = 0.1556$ $K_2 = 0.2\Big(0.7(11.2734+0.1556) \Big(\frac{11.2734+0.1556}{10}-1 \Big) \Big(1-\frac{(11.2734+0.1556)}{50} \Big) \Big)$ $K_2 = 0.1763$ $x_{i+1} = 11.2734 + \frac{0.1556+0.1763}{2} = 11.4394$ $t_{i+1} = 0.4 + 0.2 = 0.6$

literal c

Según los resultados de las tres iteraciones anteriores, la población crece. El resultado es acorde al concepto descrito en el enunciado para poblaciones mayores al valor mínimo de extinción. En el ejercicio se usó 11 como valor inicial.

Esto se puede comprobar usando el algoritmo y teniendo los resultados presentados en el literal d

literal d

Los resultados tabulados con el algoritmo son:

EDO con Runge-Kutta 2 Orden
 [ti, xi, K1, K2]
[[ 0.         11.          0.          0.        ]
 [ 0.2        11.12785958  0.12012     0.13559915]
 [ 0.4        11.2734037   0.13660392  0.15448432]
 [ 0.6        11.43943235  0.15566412  0.17639319]
 [ 0.8        11.629282    0.17778585  0.20191345]
 [ 1.         11.84695218  0.20356688  0.23177349]
...
[ 5.6        48.42689825  1.26594886  0.67489419]
 [ 5.8        49.04060051  0.81966583  0.4077387 ]
 [ 6.         49.41989811  0.5143157   0.2442795 ]]

La gráfica del ejercicio para 30 muestras es:

Instrucciones en Python

# 3Eva_2024PAOII_T3 EDO efecto Allee en poblaciones pequeñas
# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0, entrega tabla[xi,yi,K1,K2]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    # Runge Kutta de 2do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (1/2)*(K1+K2)
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2]
    return(tabla)
# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
# d1y = y' = f, d2y = y'' = f'
r = 0.7
K = 10
A = 50
d1x = lambda t,x: r*x*(x/A-1)*(1-x/K)
t0 = 0
x0 = 11
h  = 0.2
muestras = 30

# PROCEDIMIENTO
tabla = rungekutta2(d1x,t0,x0,h,muestras)
xi = tabla[:,0]
yiRK2 = tabla[:,1]

# SALIDA
# np.set_printoptions(precision=4)
print( 'EDO con Runge-Kutta 2 Orden')
print(' [ti, xi, K1, K2]')
print(tabla)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi,yiRK2)
plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[t0,x0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='[ti,xi] Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.grid()
plt.show() #comentar para la siguiente gráfica