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 \left(\frac{x}{K}-1 \right) \left(1-\frac{x}{A} \right) \frac{dx}{dt} = 0.7 x \left(\frac{x}{10}-1 \right) \left(1-\frac{x}{50} \right)

siendo h = 0.2

K_1 = 0.2 f(t_i,x_i) = 0.2\left(0.7 x \left(\frac{x}{10}-1 \right) \left(1-\frac{x}{50} \right) \right) K_2 = 0.2 f(t_i+0.2, x_i + K_1) = 0.2\left(0.7(x+K_1) \left(\frac{x+K_1}{10}-1 \right) \left(1-\frac{x+K_1}{50} \right) \right) 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\left(0.7 (11) \left(\frac{11}{10}-1 \right) \left(1-\frac{11}{50} \right) \right) = 0.1201 K_2 = 0.2\left(0.7(11+0.1201) \left(\frac{11+0.1201}{10}-1 \right) \left(1-\frac{11+0.1201}{50} \right) \right) 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\left(0.7 (11.1278) \left(\frac{11.1278}{10}-1 \right) \left(1-\frac{11.1278}{50} \right) \right) = 0.1366 K_2 = 0.2\left(0.7(11.1278+0.1366) \left(\frac{11.1278+0.1366}{10}-1 \right) \left(1-\frac{11.1278+0.1366}{50} \right) \right) 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\left(0.7 (11.2734) \left(\frac{11.2734}{10}-1 \right) \left(1-\frac{11.2734}{50} \right) \right) = 0.1556 K_2 = 0.2\left(0.7(11.2734+0.1556) \left(\frac{11.2734+0.1556}{10}-1 \right) \left(1-\frac{11.2734+0.1556}{50} \right) \right) 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: EDO dy/dx, Runge-Kutta con 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,
                vertabla=False,precision=6):
    '''solucion a EDO dy/dx, con Runge Kutta de 2do orden
    d1y es la expresi n dy/dx, tambien planteada como f(x,y),
    valores iniciales: x0,y0, tama o de paso h.
    muestras es la cantidad de puntos a calcular. 
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    tabla[0] = [x0,y0,0,0] # incluye el punto [x0,y0]
    
    xi = x0 # valores iniciales
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2]
       
    if vertabla==True:
        np.set_printoptions(precision)
        print( 'EDO dy/dx con Runge-Kutta 2 Orden')
        print('i, [xi,     yi,     K1,    K2]')
        for i in range(0,tamano,1):
            print(i,tabla[i])

    return(tabla)

# PROGRAMA PRUEBA -------------------
# 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
print( 'EDO con Runge-Kutta 2 Orden')
print(' [ti, xi, K1, K2]')
print(tabla)

# GRAFICA
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()