Modelo predador-presa (Grupo 16B)

De MateWiki
Revisión del 14:13 5 mar 2014 de Enrique941 (Discusión | contribuciones) (→‎Resolución del modelo mediante el método de Euler)

Saltar a: navegación, buscar
Warning.png 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.
Trabajo realizado por estudiantes
Título Modelo predador-presa. Grupo 16
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Javier Díez Olaya 121
Javier Lozano Aragoneses 248
Enrique Martínez Mur 271

Begoña Bigeriego Alvarez 637

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


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

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

border border border

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:

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:

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelo_predador-presa_(Grupo_16B)&oldid=10233»