Modelo predador-presa (Grupo 16B)
 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 | |
Contenido
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
% Datos del problema
t0=0;
tf=300;
%Coeficientes
A1=0.4; A2=0.3; A3=0.35;
B1=0.3; B2=0.05;
D=0.1;
C1=0.28; C2=0.045;
%Condiciones iniciales
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')
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')
3 Resolución del modelo mediante el método de Euler
t0=0;
tf=300;
%Variables
A1=0.5; A2=0.25; A3=0.3;
B1=0.4; B2=0.07;
D=0.05;
C1=0.38; C2=0.045;
%Condiciones iniciales
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')
hold off
pm=max(X(1,:));
d1m=max(X(2,:));
d2m=max(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')
4 Resolución del modelo mediante el método de Runge-Kutta de cuarto orden
% 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
