Diferencia entre revisiones de «Modelo predador-presa (Grupo 16B)»

Revisión del 14:40 5 mar 2014

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

1 Interpretación

Vito Volterra (1860-1940), físico y matemático italiano. Volterra destacó por desarrollar la solución a las ecuaciones integrales de límites variables que lleva su nombre. Alfred James Lotka (1880-1949), químico, demógrafo y matemático norteamericano. Trabajó con la misma ecuación logística que Volterra pero con el fin de describir una reacción química en la cual las concentraciones oscilan. De esta forma se estableció el modelo conocido hoy como modelo de Lotka-Volterra y que sigue representando la base de los estudios teóricos acerca de la dinámica de poblaciones.

Según nuestro caso, se supone que cantidad de comida no es un aspecto limitante pues es suficiente para la alimentación de las presas y que éstas tienen una tasa de natalidad según Maltus y una mortalidad dependiente de los depredadores:

[math] X_p'(t) = AX_p(t) [/math]:

[math] X_p´(t)=-BX_d(t)X_p(t) [/math]

Así tenemos la tasa de crecimiento de la presa con relación a la presa: [math] X_p'(t) = AX_p(t)-BX_d(t)X_p(t) [/math]

De la misma forma, los depredadores dependerán del desarrollo y la existencia de las presas, por eso se desarrollarán de la forma: [math] F'(t)=-CF(t)+DR(t)F(t) [/math]

En términos de la dinámica de población nos encontramos con que tenemos que analizar la coexistencia entre las siguientes 3 especies:

• [math]x_1[/math]: Presa
• [math]x_2[/math]: Predador 1
• [math]x_3[/math]: Predador 2

[math] \left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_2-A_3x_1x_3, t\gt0\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_1x_2-Dx_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+Dx_2x_3\\ x_1(0)=p_{0},x_2(0)=d_{0},x_3(0)=e_{0}\end{matrix}\right. [/math]

[math]A_1 [/math]= natalidad de la presa ([math]x_1[/math]).

[math]A_2[/math] y [math]A_3[/math] = mortalidad de la presa en función de la interacción con los predadores.

[math]B_1[/math] y [math]C_1[/math] = tasa de mortalidad de los predadores.

[math]B_2[/math] y [math]C_2[/math] = tasa de natalidad de los predadores dependiente de la interacción con la presa.

[math]D [/math] = descenso de población de predadores debido a su interacción.

[math]x_1[/math]crecerá en función de su tasa de natalidad, y decrecerá dependiendo de las muertes ocasionadas por ambos tipos de predadores.

[math]x_2[/math] y [math]x_3[/math] por su parte crecerán debido a su natalidad, la cual dependerá de su alimentación (las presas [math]x_1[/math]), y decrecerán tanto por su mortalidad natural como por la mortalidad provocada por el otro predador.

2 Resolución del modelo mediante el método de Euler modificado

Resolvemos numéricamente a continuación mediante Euler Modificado para periodos de 100 y 300 años. Se introducirán las constantes dadas, y las condiciones iniciales, que para nuestro caso son de 3.5 millones de presas, 1 millón de depredadores tipo 1 y 1.2 millones de depredadores tipo 2. Con todo esto se obtienen las gráficas con las que se interpretará el problema.

2.1 Código

t0=0;
tf=100;
A1=0.4;
A2=0.3;
A3=0.35;
B1=0.3;
B2=0.05;
D=0.1;
C1=0.28;
C2=0.045;
p0=3.5;
d0=1;
e0=1.2;
Y0=[p0 d0 e0]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización 
Y(:,1)=Y0;
YY=Y0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*YY(1)-A2*YY(1)*YY(2)-A3*YY(1)*YY(3);-B1*YY(2)+B2*YY(1)*YY(2)-D*YY(2)*YY(3);-C1*YY(3)+C2*YY(1)*YY(3)-D*YY(2)*YY(3)];
% Calculamos k2=f(t_n+h,y_n+h*k1)
k2=[A1*(YY(1)+h*k1(1))-A2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-A3*(YY(1)+h*k1(1))*(YY(3)+h*k1(3));-B1*(YY(2)+h*k1(2))+B2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3));-C1*(YY(3)+h*k1(3))+C2*(YY(1)+h*k1(1))*(YY(3)+h*k1(3))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3))];
YY=YY+(h/2)*(k1+k2);
Y(:,n+1)=YY;
end

figure (1)
hold on
plot(t,Y(1,:),'r')
plot(t,Y(2,:),'b')
plot(t,Y(3,:),'g')
legend('Presas', 'Primer tipo depredadores', 'Segundo tipo depredadores')
hold off

figure (2)
plot(Y(1,:),Y(2,:))
legend('Trayectoria presas, primer tipo depredadores','Location','best')

figure (3)
plot(Y(1,:),Y(3,:))
legend('Trayectoria presas, segundo tipo depredadores','Location','best')

figure(4)
plot(Y(2,:),Y(3,:))
legend('Trayectoria primer tipo depredadores, segundo tipo depredadores','Location','best')

2.2 Gráficas

Al inicio de la modelización, en el momento en el coexistían las 2 especies de depredadores con la de presas se podía apreciar un período de inestabilidad que lleva a la especie de presas a ver su población reducir de forma muy acelerada. Una vez extinguida la segunda especie de depredadores se necesita poco tiempo para que las dos especies restantes, presa y depredador tipo 1 se estabilicen en el tiempo y puedan convivir sin problema. Se nota claramente la diferencia para los dos tiempos de retorno calculados, 100 y 300 años. En el primer caso, 100 años no se puede observar una estabilidad afianzada en el tiempo, al contrario que sucede para 300 años, donde vemos que las especies de presas y depredador tipo 1 conviven, en ausencia de la especie de depredadores tipo 2, durante tiempo de forma estabilizada.

3 Resolución del modelo mediante el método de Euler

El metodo de Euler se realiza ahora para h=1 y h=0.1. En este caso, se mantendrán las condiciones iniciales, a excepción de la cantidad de depredadores tipo 2, que pasará a ser de 0.2 millones. A su vez las constantes de coexistencia que afectan a las especies también se verán alteradas. Se comprobará la estabilidad del sistema para un período de 300 años, para lo que se extraerán las gráficas convenientes.

3.1 Código

t0=0;
tf=300;
A1=0.5;
A2=0.25;
A3=0.3;
B1=0.4;
B2=0.07;
D=0.05;
C1=0.38;
C2=0.045;
p0=3.5;
d0=1;
e0=0.2;
X0=[p0 d0 e0]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización 
X(:,1)=X0;
XX=X0;
% Iteraciones
for n=1:N
XX=XX+h*[A1*XX(1)-A2*XX(1)*XX(2)-A3*XX(1)*XX(3);-B1*XX(2)+B2*XX(1)*XX(2)-D*XX(2)*XX(3);-C1*XX(3)+C2*XX(1)*XX(3)-D*XX(2)*XX(3)];
X(:,1+n)=XX;
end

figure (5)
hold on
plot(t,X(1,:),'r')
plot(t,X(2,:),'b')
plot(t,X(3,:),'k')
legend('Presas', 'Primer tipo depredadores', 'Segundo tipo depredadores')
hold off

pm=max(X(1,:))
d1m=max(X(2,:))
d2m=max(X(3,:))

pi=min(X(1,:))
d1i=min(X(2,:))
d2i=min(X(3,:))


figure (6)
plot(X(1,:),X(2,:))
legend('Trayectoria presas, primer tipo depredadores','Location','best')

figure (7)
plot(X(1,:),X(3,:))
legend('Trayectoria presas, segundo tipo depredadores','Location','best')

3.2 Gráficas

3.3 Interpretación

A la hora de comparar la modelización del sistema bajo el método de Euler y método de Euler Modificado se puede apreciar las siguientes diferencias:

• El modelo de Euler modificado, con el que se está trabajando ahora aporta una idea mas aproximada del ecosistema ya que a la hora de la modelización es un método de 2º orden, orden mayor que el Euler, el cual nos da una idea menos precisa de la estabilidad del sistema.

4 Resolución del modelo mediante el método de Runge-Kutta de cuarto orden

El Runge Kutta de 4º orden se trabajará pensando en un período horizonte de hasta 500 años, discretizaremos con h=0.1 y nuevamente las constantes inherentes a cada especie se verán modificadas, teniendo que trabajar para dos situaciones diferentes:

• 3.5 millones de presas; 0.001 millones de depredadores tipo 1; 0.0002 millones de depredadores tipo 2
• 3.5 millones de presas; 0.00001 millones de depredadores tipo 1; 0.2 millones de depredadores tipo 2

4.1 Código

% Datos del problema con los primeros valores iniciales
t0=0;
tf=500;
A1=0.4; A2=0.3; A3=0.4;
B1=0.37; B2=0.05;
D=0.1;
C1=0.28; C2=0.07;
%Condiciones iniciales
p01=3.5;
d01=0.001;
e01=0.0002;
Y0=[p01 d01 e01]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización 
Y(:,1)=Y0;
YY=Y0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*YY(1)-A2*YY(1)*YY(2)-A3*YY(1)*YY(3);-B1*YY(2)+B2*YY(1)*YY(2)-D*YY(2)*YY(3);-C1*YY(3)+C2*YY(1)*YY(3)-D*YY(2)*YY(3)];
% Calculamos k2=f(t_n+h/2,y_n+(1/2)*h*k1)
k2=[A1*(YY(1)+(h/2)*k1(1))-A2*(YY(1)+(h/2)*k1(1))*(YY(2)+(h/2)*k1(2))-A3*(YY(1)+(h/2)*k1(1))*(YY(3)+(h/2)*k1(3));-B1*(YY(2)+(h/2)*k1(2))+B2*(YY(1)+(h/2)*k1(1))*(YY(2)+(h/2)*k1(2))-D*(YY(2)+(h/2)*k1(2))*(YY(3)+(h/2)*k1(3));-C1*(YY(3)+(h/2)*k1(3))+C2*(YY(1)+(h/2)*k1(1))*(YY(3)+(h/2)*k1(3))-D*(YY(2)+(h/2)*k1(2))*(YY(3)+(h/2)*k1(3))];
% Calculamos k3=f(t_n+h/2,y_n+(1/2)*h*k2)
k3=[A1*(YY(1)+(h/2)*k2(1))-A2*(YY(1)+(h/2)*k2(1))*(YY(2)+(h/2)*k2(2))-A3*(YY(1)+(h/2)*k2(1))*(YY(3)+(h/2)*k2(3));-B1*(YY(2)+(h/2)*k2(2))+B2*(YY(1)+(h/2)*k2(1))*(YY(2)+(h/2)*k2(2))-D*(YY(2)+(h/2)*k2(2))*(YY(3)+(h/2)*k2(3));-C1*(YY(3)+(h/2)*k2(3))+C2*(YY(1)+(h/2)*k2(1))*(YY(3)+(h/2)*k2(3))-D*(YY(2)+(h/2)*k2(2))*(YY(3)+(h/2)*k2(3))];
% Calculamos k4=f(t_n+h,y_n+h*k3)
k4=[A1*(YY(1)+h*k3(1))-A2*(YY(1)+h*k3(1))*(YY(2)+h*k3(2))-A3*(YY(1)+h*k3(1))*(YY(3)+h*k3(3));-B1*(YY(2)+h*k3(2))+B2*(YY(1)+h*k3(1))*(YY(2)+h*k3(2))-D*(YY(2)+h*k3(2))*(YY(3)+h*k3(3));-C1*(YY(3)+h*k3(3))+C2*(YY(1)+h*k3(1))*(YY(3)+h*k3(3))-D*(YY(2)+h*k3(2))*(YY(3)+h*k3(3))];
YY=YY+(h/6)*(k1+2*k2+2*k3+k4);
Y(:,n+1)=YY;
end
figure (8)
hold on
plot(t,Y(1,:),'r','LineWidth',1.4)
plot(t,Y(2,:),'b','LineWidth',1.4)
plot(t,Y(3,:),'g','LineWidth',1.4)
hold off

% Datos del problema con los segundos valores iniciales
p02=3.5;
d02=0.00001;
e02=0.2;
Z0=[p02 d02 e02]';
% Inicialización 
Z(:,1)=Z0;
ZZ=Z0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*ZZ(1)-A2*ZZ(1)*ZZ(2)-A3*ZZ(1)*ZZ(3);-B1*ZZ(2)+B2*ZZ(1)*ZZ(2)-D*ZZ(2)*ZZ(3);-C1*ZZ(3)+C2*ZZ(1)*ZZ(3)-D*ZZ(2)*ZZ(3)];
% Calculamos k2=f(t_n+h/2,y_n+(1/2)*h*k1)
k2=[A1*(ZZ(1)+(h/2)*k1(1))-A2*(ZZ(1)+(h/2)*k1(1))*(ZZ(2)+(h/2)*k1(2))-A3*(ZZ(1)+(h/2)*k1(1))*(ZZ(3)+(h/2)*k1(3));-B1*(ZZ(2)+(h/2)*k1(2))+B2*(ZZ(1)+(h/2)*k1(1))*(ZZ(2)+(h/2)*k1(2))-D*(ZZ(2)+(h/2)*k1(2))*(ZZ(3)+(h/2)*k1(3));-C1*(ZZ(3)+(h/2)*k1(3))+C2*(ZZ(1)+(h/2)*k1(1))*(ZZ(3)+(h/2)*k1(3))-D*(ZZ(2)+(h/2)*k1(2))*(ZZ(3)+(h/2)*k1(3))];
% Calculamos k3=f(t_n+h/2,y_n+(1/2)*h*k2)
k3=[A1*(ZZ(1)+(h/2)*k2(1))-A2*(ZZ(1)+(h/2)*k2(1))*(ZZ(2)+(h/2)*k2(2))-A3*(ZZ(1)+(h/2)*k2(1))*(ZZ(3)+(h/2)*k2(3));-B1*(ZZ(2)+(h/2)*k2(2))+B2*(ZZ(1)+(h/2)*k2(1))*(ZZ(2)+(h/2)*k2(2))-D*(ZZ(2)+(h/2)*k2(2))*(ZZ(3)+(h/2)*k2(3));-C1*(ZZ(3)+(h/2)*k2(3))+C2*(ZZ(1)+(h/2)*k2(1))*(ZZ(3)+(h/2)*k2(3))-D*(ZZ(2)+(h/2)*k2(2))*(ZZ(3)+(h/2)*k2(3))];
% Calculamos k4=f(t_n+h,y_n+h*k3)
k4=[A1*(ZZ(1)+h*k3(1))-A2*(ZZ(1)+h*k3(1))*(ZZ(2)+h*k3(2))-A3*(ZZ(1)+h*k3(1))*(ZZ(3)+h*k3(3));-B1*(ZZ(2)+h*k3(2))+B2*(ZZ(1)+h*k3(1))*(ZZ(2)+h*k3(2))-D*(ZZ(2)+h*k3(2))*(ZZ(3)+h*k3(3));-C1*(ZZ(3)+h*k3(3))+C2*(ZZ(1)+h*k3(1))*(ZZ(3)+h*k3(3))-D*(ZZ(2)+h*k3(2))*(ZZ(3)+h*k3(3))];
ZZ=ZZ+(h/6)*(k1+2*k2+2*k3+k4);
Z(:,n+1)=ZZ;
end
figure (9)
hold on
plot(t,Z(1,:),'r','LineWidth',1.4)
plot(t,Z(2,:),'b','LineWidth',1.4)
plot(t,Z(3,:),'g','LineWidth',1.4)
hold off

4.2 Gráficas

Resultado del modelo con valores iniciales [math] p_0 =3.5[/math] millones de presas [math] d_0=0,001[/math] millones de predadores de un tipo [math]e_0=0,0002[/math]millones de predadores de otro tipo

Resultado del modelo con valores iniciales [math] p_0 =3.5[/math] millones de presas [math] d_0=0,00001[/math] millones de predadores de un tipo [math]e_0=0,2[/math]millones de predadores de otro tipo

4.3 Interpretación

Llevando a comparar este método de Runge Kutta 4º orden con el método de Euler se puede llegar a las conclusiones siguientes:

• Como Euler Modificado es un método de segundo orden y Runge-Kutta, con el que trabajamos, de 4º orden, Euler Modificado es mas inestable, y aporta ideas menos precisas en cuanto a la estabilidad del ecosistema.