Ver el código de «Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)»

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== 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

<br />

* 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

<br />

* 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

<br />

* 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

<br />

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

<br />

:<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>0</math>

== Resolución numérica ==

=== 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:

{{matlab|codigo=
%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')
}}

[[Archivo: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.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo: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á:

{{matlab|codigo=
%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')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo: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.

=== 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:

{{matlab|codigo=
%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')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo: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:

[[Archivo:TablaEuler12.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%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')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo: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:

[[Archivo:TablaEuler01.JPG]]

=== 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:

{{matlab|codigo=
%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
}}

[[Archivo: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.

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]