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

Revisión actual del 11:15 18 may 2014

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

2.2.1 Gráfica para 100 años

Evolución de las tres especies a lo largo del tiempo

Relación entre presas y predadores del tipo 1

Relación entre presas y predadores del tipo 2

Relación entre los dos tipos de predadores

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.

2.2.2 Gráfica para 300 años

Evolución de las tres especies a lo largo del tiempo

Relación entre presas y predadores del tipo 1

Relación entre presas y predadores del tipo 2

Relación entre los dos tipos de predadores

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.2.1 Gráficas correspondientes a [math]h=0.1 [/math]

Evolución de las tres especies a lo largo del tiempo

Relación entre presas y depredadores del tipo 1

Relación entre presas y depredadores del tipo 2

En el primer caso, gráfica 1, no se alcanzará la estabilidad , ya que una de las especies se reducirá al mínimo ; mientras que las otras dos crecerán de manera exponencial.

Por otro lado, la estabilidad entre la especie de presas y depredadores tipo 1 mantiene una estabilidad cíclica a lo largo del tiempo, dónde a medida que la población de uno disminuye la del otro ve su número de ejemplares aumentar y viceversa.

En cambio entre los depredadores tipo 2, los que se extinguen, y las presas se observa que la dependencia de los depredadores ante la falta de presas que comer les llevaría a la extinción como así finalmente sucede

3.2.2 Gráficas correspondientes a [math]h=1 [/math]

Evolución de las tres especies a lo largo del tiempo

Relación entre presas y depredadores del tipo 1

Relación entre presas y depredadores del tipo 2

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;
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')
plot(t,Y(2,:),'b')
plot(t,Y(3,:),'g')
legend('Presas', 'Primer tipo depredadores', 'Segundo tipo depredadores')
hold off


pm0=max(Y(1,:))
d1m0=max(Y(2,:))
d2m0=max(Y(3,:))

pi0=min(Y(1,:))
d1i0=min(Y(2,:))
d2i0=min(Y(3,:))

% 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')
plot(t,Z(2,:),'b')
plot(t,Z(3,:),'g')
legend('Presas', 'Primer tipo depredadores', 'Segundo tipo depredadores')
hold off

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

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

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

Bajo la modelización mediante Runge Kutta se aprecia con mayor precisión la estabilidad del ecosistema con las dos especies finales. Para este método se ha planteado dos situaciones diferentes que evalúan de diferente forma el ecosistema.

Después, con las presas en el mismo número de individuos y una especie de depredadores reducida a una milésima parte y la otra aumentada mil veces la cantidad anterior llega a apreciarse como la estabilidad del ecosistema es alcanzada de forma mas rápida sin que la especie de menor número interfiera en la dominante.

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.