Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

De MateWiki
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. Dinámica de poblaciones. Grupo 4
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Planteamiento e interpretación

Se conoce como modelo de Volterra-Lotka el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

  • Las tres poblaciones son homogéneas
  • Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas [math]x_1(t)[/math] y [math]x_2(t)[/math], con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:
[math]A_1x_1[/math]
[math]B_1x_2[/math]
Siendo [math]A_1[/math] y [math]B_1[/math] constantes posistivas


  • La alimentación de la población de predadores [math]x_3(t)[/math] depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:
[math]C_2 x_1 x_3 + C_3 x_2 x_3[/math]
Siendo [math]C_2[/math] y [math]C_3[/math] constantes positivas


  • La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:
[math]A_2 x_1 x_3[/math]
[math]B_2 x_2 x_3[/math]
Siendo [math]A_2[/math] y [math]B_2[/math] constantes positivas


  • Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:
[math]C_1 x_3[/math]
Siendo [math]C_1[/math] una constante positiva


Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:


[math] \left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right. [/math]
Siendo [math]t\gt0[/math]

2 Resolución numérica

2.1 Método de Euler modificado

Resolvemos mediante el método de Euler modificado con [math]h=0.1[/math] en un intervalo de tiempo de 0 a 100 años y tomando [math]A_1=0.35[/math], [math]A_2=0.6[/math], [math]B_1=0.3[/math], [math]B_2=0.5[/math], [math]C_1=0.37[/math], [math]C_2=0.04[/math] y [math]C_3=0.035[/math]. Dadas las condiciones iniciales [math]p_{0}=2[/math] millones de presas [math]x_1[/math], [math]q_{0}=1.4[/math] millones de presas [math]x_2[/math] y [math]d_{0}=1[/math] millón de predadores, el código en Matlab es el siguiente: 

%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;

x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
    k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
    k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
    xx=xx+(h/2)*(k1+k2);
    x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')


EvolucionEcosistema.jpg

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

EvolucionPresas1Predadores.jpg EvolucionPresas2Predadores.jpg EvolucionPresas1Presas2.jpg

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

FALTA EXPLICAR LA TERCERA!!!

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será: 

%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
    k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
    k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
    xx=xx+(h/2)*(k1+k2);
    x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')


EvolucionEcosistema2.jpg EvolucionPresa1Predador2.jpg EvolucionPresa2Predador2.jpg EvolucionPresa1Presa22.jpg

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

2.2 Método de Euler

Resolvemos mediante el método de Euler con [math]h=1[/math] y con [math]h=0.1[/math], en un intervalo de tiempo de 0 a 300 años y tomando [math]A_1=0.35[/math], [math]A_2=0.6[/math], [math]B_1=0.3[/math], [math]B_2=0.5[/math], [math]C_1=0.37[/math], [math]C_2=0.04[/math] y [math]C_3=0.035[/math]. Dadas las condiciones iniciales [math]p_{0}=0.8[/math] millones de presas [math]x_1[/math], [math]q_{0}=2.4[/math] millones de presas [math]x_2[/math] y [math]d_{0}=0.2[/math] millones de predadores, el código en Matlab es el siguiente: 

%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
    xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
    x(:,n+1)=xx;   
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off 
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')


EvolucionEcosistemaEuler.jpg EvolucionPresa1PredadorEuler.jpg EvolucionPresa2PredadorEuler.jpg

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Tomando [math]h=0.1[/math] el codigo será: 

%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
    xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
    x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off 
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')


EvolucionEcosistema01.jpg EvolucionPresa1Predador01.jpg EvolucionPresa2Predador01.jpg

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior.

FALTA METER LAS TABLAS!!!

2.3 Método de Runge-Kutta

Resolvemos mediante el método de Runge-Kutta con [math]h=0.1[/math], en un intervalo de tiempo de 0 a 500 años y tomando [math]A_1=0.4[/math], [math]A_2=0.7[/math], [math]B_1=0.2[/math], [math]B_2=0.4[/math], [math]C_1=0.37[/math], [math]C_2=0.04[/math] y [math]C_3=0.03[/math]. Dadas las condiciones iniciales [math]p_{0}=3.5[/math] millones de presas [math]x_1[/math], [math]q_{0}=0.2[/math] millones de presas [math]x_2[/math] y [math]d_{0}=0.4[/math] millones de predadores, el código en Matlab es el siguiente: 

%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
    k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
    k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
    k3=[A1*(xx(1)+(h/2)*k2(1))-A2*(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(1)),B1*(xx(2)+(h/2)*k2(2))-B2*(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(2)),-C1*(xx(3)+(h/2)*k2(3))+C2*(xx(1)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))+C3*(xx(2)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))]';
    k4=[A1*(xx(1)+h*k3(1))-A2*(xx(1)+h*k3(1))*(xx(3)+h*k3(1)),B1*(xx(2)+h*k3(2))-B2*(xx(2)+h*k3(2))*(xx(3)+h*k3(2)),-C1*(xx(3)+h*k3(3))+C2*(xx(1)+h*k3(3))*(xx(3)+h*k3(3))+C3*(xx(2)+h*k3(3))*(xx(3)+h*k3(3))]';
    xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
    x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off


EvolucionEcosistemaRK.jpg

La situación que nos muestra esta gráfica es diferente a las anteriores. Se aprecia que la población [math]x_2[/math] tiende a extinguirse. Esto tiene dos consecuencias, el aumento de los predadores y la estabilizacion de la presa [math]x_1[/math] que hasta el momento de la desaparicion de [math]x_2[/math] había estado creciendo

FALTA METER LOS MÁXIMOS Y MÍNIMOS!!!

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelo_Predador-Presa._Dinámica_de_poblaciones_(Grupo_4)&oldid=9327»