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

De MateWiki
Saltar a: navegación, buscar
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;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
    k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
    k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
    xx=xx+h*(k1+k2)/2;
    x1(n+1)=xx(1);
    x2(n+1)=xx(2);
    x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'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(x1,x3,'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(x2,x3,'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(x1,x2,'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')


EvolucionEcosistema100.jpg

En la gráfica obtenida se muestra la evolución en el tiempo de las poblaciones de presas y predadores. 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. De igual forma, se observa que una de las poblaciones de presas, que comienza teniendo una mayor población, tiende a disminuir, mientras que la que empezó con menor población aumenta. Los máximos y los mínimos que se dan en la población de predadores se mantienen constantes en el tiempo.

EvolucionPresas1Predadores100.jpg EvolucionPresas2Predadores100.jpg EvolucionPresas1Presas2100.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. Su forma, sin embargo, no nos permite decir si estamos ante un ecosistema estable.

En la tercera gráfica se aprecia la evolución de la población de un tipo de presas en función de la población del otro tipo de presas.

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;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
    k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
    k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
    xx=xx+h*(k1+k2)/2;
    x1(n+1)=xx(1);
    x2(n+1)=xx(2);
    x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'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(x1,x3,'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(x2,x3,'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(x1,x2,'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')


EvolucionEcosistema300.jpg EvolucionPresa1Predador300.jpg EvolucionPresa2Predador300.jpg EvolucionPresa1Presa2300.jpg

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

La evolución de las presas en el tiempo continúa la tendencia que se aprecia en los primeros 100 años. Tenderá a extinguirse un tipo de presas debido a que el depredador se alimentará más de ellas, provocando esto que el otro tipo de presas aumente progresivamente su población.

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.

Para esta anchura de paso obtenemos la siguiente tabla de máximos y mínimos asociada:

TablaEuler12.JPG

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

Se presenta en la primera grafica de nuevo la evolución en el tiempo de los dos tipos de presas y predador. Se aprecia lo indicado en la gráfica obtenida con Euler modificado: aumentan los dos tipos de presas, teniendo en cuenta que uno de ellos dispone de muchos más ejemplares que el otro, los predadores comienzan entonces a aumentar dado que disponen de más alimento para sobrevivir y esto provoca que disminuyan las presas. Se observa que éstos alcanzan un máximo donde las presas disminuyen hasta su mínimo y por tanto apenas existe alimento. Es en ese momento cuando empiezan a competir en ellos y disminuye su población.

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, ya que permite que pueda aumentarse la longitud de paso sin acumular tanto error.

En este caso obtenemos la siguiente tabla de máximos y mínimos:

TablaEuler01.JPG

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;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
    k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
    k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
    k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
    k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
    xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
    x1(n+1)=xx(1);
    x2(n+1)=xx(2);
    x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'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


EvolucionEcosistemaRK2.jpg

Es con este método con el que mejor se aprecia la evolución de las especies de presas (llegando a extinguirse una de ellas, la presa tipo 2) y de la especie predadora.

TablaRK.JPG

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