Modelo predador-presa (Grupo 7)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo predador-presa (Grupo 7) |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2015-16 |
| Autores | Gómez Apiñániz, Adrián
Herranz Rodas, Ignacio Ragolta Villarroya, Ana Reinoso Muñoz, Cristina |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En esta página se discute un modelo matemático que trata de resolver la dinámica de poblaciones de especies competidoras, teniendo en cuenta el estudio realizado a principios del siglo pasado por Alfred J. Lotka y Vito Volterra, el cual consiste en modelizar la situación depredador-presa con ayuda de un sistema de ecuaciones diferenciales de primer orden no lineales, correspondientes a la variación de ambas poblaciones por unidad de tiempo.
Contenido
- 1 . Interpretación del modelo
- 2 . Resolución del problema mediante el método de Euler modificado
- 3 . Resolución del problema mediante el método de Euler
- 4 . Resolución del problema mediante el método de Runge-Kutta
- 4.1 . Para [math]p_{0}=3,5 millones[/math] de presas, [math]d_{0}=0,001 millones[/math] de predadores de la primera especie, [math]e_{0}=0,0002 millones[/math] de predadores de la segunda especie
- 4.2 . Para [math]p_{0}=3,5 millones[/math] de presas, [math]d_{0}=0,00001 millones[/math] de predadores de la primera especie, [math]e_{0}=0,2 millones[/math] de predadores de la segunda especie
- 5 . Comparación de resultados
1 . Interpretación del modelo
El problema de valor inicial a estudiar, es el adjunto:
[math] \left\{\begin{matrix} \frac{dx_{1}}{dt}=A_{1}x_{1}-A_{2}x_{1}x_{2}-A_{3}x_{1}x_{3}\\ \frac{dx_{2}}{dt}=-B_{1}x_{2}+B_{2}x_{1}x_{2}-Dx_{2}x_{3}\\ \frac{dx_{3}}{dt}=-C_{1}x_{3}+C_{2}x_{1}x_{3}-Dx_{2}x_{3}\\x_{1}(0)=p_{0},x_{2}(0)=d_{0},x_{3}(0)=e_{0}\\t\gt 0 \end{matrix}\right. [/math]
A priori, se puede observar que en esta ocasión el modelo de Lotka-Volterra ha sido modificado, al haber añadido una tercera población que actuará como depredador. Por tanto, el problema a estudiar presenta una única presa frente a dos depredadores.
Además, teniendo en cuenta las ecuaciones adjuntas y la hipótesis de que la cantidad de alimentos para la presa es siempre suficiente, se ha realizado la siguiente interpretación:
[math]x_{1}[/math]: número de presas
[math]x_{2}[/math]: número de depredadores de la primera especie
[math]x_{3}[/math]: número de depredadores de la segunda especie
[math]A_{1}[/math]: constante de la que depende la tasa de natalidad de las presas
[math]A_{2}[/math]: constante de interacción entre las presas y el primer depredador
[math]A_{3}[/math]: constante de interacción entre las presas y el segundo depredador
[math]B_{1}[/math]: constante de la que depende la tasa de mortalidad de los depredadores de la primera especie
[math]B_{2}[/math]: constante de interacción entre el primer depredador y las presas
[math]D[/math]: constante de interacción entre ambos depredadores
[math]C_{1}[/math]: constante de la que depende la tasa de mortalidad de los depredadores de la segunda especie
[math]C_{2}[/math]: constante de interacción entre el segundo depredador y las presas
[math]p_{0}[/math]: número inicial de presas
[math]d_{0}[/math]: número inicial de depredadores de la primera especie
[math]e_{0}[/math]: número inicial de depredadores de la segunda especie
Siendo todas las constantes positivas.
La primera ecuación, representa la variación del número de presas a lo largo del tiempo, la cual depende del número de nacimientos menos el número de muertes. El número de nacimientos se rige por una tasa de natalidad que según la ley malthusiana depende a su vez de un coeficiente [math]A_{1}[/math] y del número de presas existentes. Por otro lado, el número de muertes depende de las interacciones entre cada uno de los depredadores y las presas, modelado con ayuda de los coeficientes [math]A_{2}[/math] y [math]A_{3}[/math].
La segunda ecuación, representa la variación del número de depredadores de la primera especie a lo largo del tiempo. Esta variación a su vez depende de tres factores:
-La tasa de mortalidad: que queda reflejada según la constante [math]B_{1}[/math] y el número de depredadores de la primera especie existentes.
-La interacción con la presa: que queda reflejada según la constante [math]B_{2}[/math] y el número tanto de depredadores de la primera especie como de presas. Se puede ver que a mayor número de presas, mayor crecimiento de la población de los depredadores.
-La interacción con la otra especie depredadora: que varía según la constante [math]D[/math] y el número de depredadores de cada una de las especies. Esto refleja que cuantos más depredadores existan de ambas especies, menor será el crecimiento de cada una de ellas, pero además se puede ver que cuanto mayor sea la diferencia entre el número de individuos de ambas especies predadoras, menor será el número de muertes.
Análogamente, la tercera ecuación obtiene la tasa de crecimiento para la segunda especie depredadora, modificando las constantes [math]B_{1}[/math] por [math]C_{1}[/math] y [math]B_{2}[/math] por [math]C_{2}[/math], pero manteniendo el factor de interacción entre especies.
2 . Resolución del problema mediante el método de Euler modificado
Teniendo en cuenta el modelo planteado, se ha resuelto por el método de Euler modificado para [math]T=[0,100][/math] y [math]T=[0,300][/math] años. Para ambos casos se ha utilizado un único programa que, al iniciarlo, solicita el instante final para el que se quiere obtener la resolución y obtiene cuatro gráficas, siendo la primera la variación del número de individuos de las tres poblaciones respecto al tiempo y siendo las tres restantes la variación de las poblaciones entre sí, dos a dos. El programa es el siguiente:
%Datos del problema
ne=3;
t0=0;
tn=input('Introduce el instante final: '); %Serán 100 o 300 años
y0=[3.5;1;1.2];
ftexto='[0.4*y(1)-0.3*y(1)*y(2)-0.35*y(1)*y(3);-0.3*y(2)+0.05*y(1)*y(2)-0.1*y(2)*y(3);-0.28*y(3)+0.045*y(1)*y(3)-0.1*y(2)*y(3)]';
f=inline(ftexto,'t','y');
%Así se convierte el texto en una función
%Discretización
h=0.1;
N=round((tn-t0)/h);
%Creación del vector t
t=t0:h:tn;
%Creación el vector y de ceros
y=zeros(ne,N+1);
%Se adjudia el primer valor de y
y(:,1)=y0;
%Aplicación del esquema numérico
for i=1:N;
K=f(t(i),y(:,i));
y(:,i+1)=y(:,i)+h*K;
%f(t(i)+0.5*h,y(:,i)+0.5*K*h);
end
%Gráfica conjunta de las tres funciones
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
legend('Número de presas','Número de depredadores de la primera especie','Número de depredadores de la segunda especie','Location','Best')
title('Ecosistema')
xlabel('Tiempo')
ylabel('Millones de individuos')
hold off
%Dibujo las gráficas de las trayectorias
figure
plot(y(1,:),y(2,:))
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
figure
plot(y(1,:),y(3,:))
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
figure
plot(y(2,:),y(3,:))
title('Trayectoria')
xlabel('Número de depredadores de primera especie')
ylabel('Número de depredadores de segunda especie')
2.1 . Para [math]T=[0,100][/math] años
A continuación, se adjuntan las gráficas obtenidas para el modelo con un periodo de [math]T=[0,100][/math] años:
Como se puede comprobar en la gráfica que indica la variación del número de individuos de cada especie respecto al tiempo, el ecosistema comienza con una disminución de la población de las presas, lo que provoca que el número de individuos de las dos especies depredadoras también decrezca. Debido a esta disminución de los depredadores el número de presas aumenta hasta alcanzar su máximo: 15.0969 millones de individuos. Las especies depredadoras evolucionan de forma muy parecida durante los primeros 30 años,pero a partir de los 31 años las poblaciones de ambas especies se igualan. Desde este momento la segunda especie depredadora, se reduce hasta alcanzar su mínimo que es 1900 individuos.
- Análisis trayectoria Presa-Depredador de la Primera especie: en esta gráfica se puede observar que la relación entre la presa y el depredador de primera especie es una interacción cíclica. Inicialmente, cuando disminuyen las presas, disminuye también el número de depredadores. Se puede ver como después el número de depredadores crece de una manera menos pronunciada que el número de presas, y cuando las presas descienden, todavía los depredadores siguen aumentando hasta un punto en el que empieza su decrecimiento. Siempre que se llega al decrecimiento brusco de depredadores empieza un crecimiento intenso de presas. También cabe destacar que debido a la disminución del depredador de la segunda especie, la gráfica cambia de patrón y tiende a formar una figura elíptica que se mantendrá constante a lo largo del tiempo, como podremos comprobar para T=300, ésto ocurrirá cuando la segunda especie depredadora se reduzca hasta extinguirse.
- Análisis trayectoria Presa-Depredador de la Segunda especie: en la gráfica se comprueba que la interacción entre la presa y el depredador de la segunda especie es también cíclica inicialmente, el número de presas es mayor que el de depredadores y una disminución del número de presas irá asociada a una disminución del número de depredadores. Sin embargo llega un momento en el que el número de depredadores disminuye y al haber menos depredadores, el número de presas volverá a crecer, pegándose la gráfica a la linea horizontal.
-Análisis trayectoria Depredador de la Primera especie-Depredador de la Segunda especie: en esta gráfica se comprueba que inicialmente ambas poblaciones se distribuyen de manera similar, siendo algo superior el número de individuos de la segunda especie depredadora. Inicialmente su variación es parecida, es decir, la disminución de presas provoca que se disminuya el mismo número de predadores de ambas especies, sin embargo, el número de depredadores de la segunda especie comienza a disminuir hasta que se iguala con los de la primera especie y pasado ese límite se reduce hasta casi extinguirse lo cual hace que la gráfica en los instantes finales se encuentre pegado al eje de abscisas.
2.2 . Para [math]T=[0,300][/math] años
A continuación, se adjuntan las gráficas obtenidas para el modelo con un periodo de [math]T=[0,300][/math] años:
Como consecuencia de la desaparición de la segunda especie depredadora, mencionado en el apartado anterior a los 149 años desaparece el último miembro de la segunda especie depredadora lo que supondrá la extinción de la especie , y como se puede observar en la gráfica adjunta, la primera especie depredadora aumenta progresivamente debido a que el número de individuos de las presas también aumenta y a que no tiene especie competidora de alimento. La evolución de ambas especies es simultánea, es decir, el aumento de la presa conlleva un aumento de la especie depredadora mientras que la disminución del número de presas disminuye también el número de depredadores, tal y como describe el modelo Volterra-Lotka.
El ecosistema no es estable ya que para considerarse un ecosistema estable los valores de las especies deberían mantenerse muy similares en todo momento y como hemos analizado previamente una de las especies se extingue y las otras dos no alcanzan nunca valores estables.
Comparando las gráficas para [math]T=[0,300][/math] años frente a las de [math]T=[0,100][/math], podemos observar que para [math]T=[0,300][/math] son la continuación de las gráficas de [math]T=[0,100][/math].
- Análisis trayectoria Presa-Depredador de la Primera especie: es la única gráfica en la que se puede apreciar a simple vista una distinción frente a aquella para [math]T=[0,100][/math] años. A partir de los 100 años la relación entre estas especies se convierte en una espiral que crece con el tiempo ya que los valores aumentan continuamente.
- Análisis trayectoria Presa-Depredador de la Segunda especie: es una de las gráficas que se mantiene similar a la de [math]T=[0,100][/math] años, puesto que, como ya se ha mencionado con anterioridad, la población de depredadores de la segunda especie se extingue, haciendo que la función se confunda con el eje de abscisas.
-Análisis trayectoria Depredador de la Primera especie-Depredador de la Segunda especie: al igual que para la trayectoria anterior, se puede decir que es muy similar a la trayectoria con [math]T=[0,100][/math] años. Una vez se extingue el segundo depredador la función se encuentra pegada al eje de abcisas.
3 . Resolución del problema mediante el método de Euler
Una vez resuelto el problema mediante el método de Euler modificado, se ha planteado su resolución empleando Euler. Para ello, se ha seleccionado un periodo de tiempo de [math]T=[0,300][/math] años y se ha tratado de resolver el problema para dos pasos distintos, por lo que el programa solicita este dato. Una vez incluido, se obtienen tres gráficas, siendo la primera la variación del número de individuos de las tres especies respecto al tiempo y las dos restantes la trayectoria de las presas frente a los depredadores de la primera y segunda especie. Además, el programa calcula los valores máximos y mínimos de cada población y el instante en el que se dan estos máximos y mínimos.
%Datos del problema
ne= 3;
t0=0;
tn=300;
y0=[3.5;1;0.2];
ftexto='[0.5*y(1)-0.25*y(1)*y(2)-0.3*y(1)*y(3);-0.4*y(2)+0.07*y(1)*y(2)-0.05*y(2)*y(3);-0.38*y(3)+0.045*y(1)*y(3)-0.05*y(2)*y(3)]';
f=inline(ftexto,'t','y');
%Así se convierte el texto en una función
%Discretización
h=input('Introduce el tamaño de paso: '); %Que será 1 o 0.1
N=round((tn-t0)/h);
%Creación del vector t
t=t0:h:tn;
%Se crea el vector y de ceros
y=zeros(ne,N+1);
%Se adjudica el primer valor de y
y(:,1)=y0;
%Aplicación del esquema numérico
for i=1:N;
y(:,i+1)=y(:,i)+h*f(t(i),y(:,i));
end
%Se realiza la gráfica conjunta de las tres funciones
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
legend('Número de presas','Número de depredadores de la primera especie','Número de depredadores de la segunda especie','Location','Best')
xlabel('Tiempo')
ylabel('Millones de individuos')
hold off
%Para calcular las trayectorias se utiliza:
figure
plot(y(1,:),y(2,:))
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
figure
plot(y(1,:),y(3,:))
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
figure
%Máximos y mínimos
maximopresa=max(y(1,:))
maximopredador1=max(y(2,:))
maximopredador2=max(y(3,:))
minimopresa=min(y(1,:))
minimopredador1=min(y(2,:))
minimopredador2=min(y(3,:))
posmaxpresa=ceil(find(y==max(y(1,:)))/3);
posmaxpred1=ceil(find(y==max(y(2,:)))/3);
posmaxpred2=ceil(find(y==max(y(3,:)))/3);
posminpresa=ceil(find(y==min(y(1,:)))/3);
posminpred1=ceil(find(y==min(y(2,:)))/3);
posminpred2=ceil(find(y==min(y(3,:)))/3);
tmaxpresa=t(posmaxpresa)
tminpresa=t(posminpresa)
tmaxpred1=t(posmaxpred1)
tminpred1=t(posminpred1)
tmaxpred2=t(posmaxpred2)
tminpred2=t(posminpred2)
3.1 . Para un paso de [math]h=1[/math]
El tamaño de paso [math]h=1[/math] no es adecuado para el método de Euler, pues es demasiado elevado y produce un error muy grande. Es por ello que, en las gráficas adjuntas de variación del número de individuos respecto al tiempo y trayectoria de las presas frente al resto de especies, los valores que se alcanzan carecen de sentido.
3.2 . Para un paso de [math]h=0.1[/math]
Como se ha mencionado anteriormente, el programa obtiene los valores máximos y mínimos de cada especie, así como el momento en el que se producen. Los datos obtenidos para el paso [math]h=0.1[/math] son los siguientes:
-Máximo número de presas: 103.5183 millones, que se produce en un tiempo de 238.8 años.
-Mínimo número de presas: 1.9024e-10 millones, que se produce en un tiempo de 246.5 años.
-Máximo número de depredadores de la primera especie: 33.4513 millones, que se produce en un tiempo de 239.5 años.
-Mínimo número de depredadores de la primera especie: 2.8502e-09 millones, que se produce en un tiempo de 298.4 años.
-Máximo número de depredadores de la segunda especie: 0.2000 millones, que se produce en un tiempo de 0 años.
-Mínimo número de depredadores de la segunda especie: 2.2093e-34 millones, que se produce en un tiempo de 299.2 años.
Además, las gráficas mostradas por el programa son las siguientes:
- Análisis variación del número de individuos respecto al tiempo: se puede comprobar que el ecosistema sería inestable debido a que los valores que alcanzan tanto las presas como los depredadores de ambas especies, son tremendamente pequeños, menores de 10^-9.
- Análisis trayectoria Presa-Depredador de la Primera especie: lo primero que se puede observar con claridad es que cuando las presas alcanzan un máximo, luego descienden de forma brusca y además, es cuando éstas alcanzan un mínimo cuando los valores del número de individuos de los depredadores de la primera especie son máximos. Además, si estudiamos los valores calculados anteriormente por el programa, se puede ver que los tiempos para los que el número de presas es mínimo y el número de depredadores de la primera especie máximo, prácticamente coinciden, siendo 246.5 y 239.5 años, respectivamente.
- Análisis trayectoria Presa-Depredador de la Segunda especie: se puede observar en la gráfica que el número de depredadores de la segunda especie desciende en todo momento hasta su extinción, excepto en dos momentos, que coinciden con el mayor crecimiento del número de presas. Aunque el programa obtiene como mínimo para el número de depredadores de la segunda especie 2.2093e-34 millones y no cero, es un valor tan pequeño que podemos afirmar que la especie se extinguirá.
4 . Resolución del problema mediante el método de Runge-Kutta
Por último, se ha planteado la resolución del problema empleando el método de Runge-Kutta para un periodo de tiempo de [math]T=[0,500][/math] años. Para este programa, se ha solicitado el número inicial de individuos de cada población, obteniendo una gráfica con la variación de cada especie a lo largo del tiempo, así como los valores máximos y mínimos de cada población y el instante en el que se producen.
%Datos del problema
ne=3;
t0=0;
tn=500;
y0=input('Introduce el número de individuos inicial para cada población como vector columna: ');%será [3.5;0.001;0.0002] o [3.5;0.00001;0.2]
ftexto='[0.4*y(1)-0.3*y(1)*y(2)-0.4*y(1)*y(3);-0.37*y(2)+0.05*y(1)*y(2)-0.1*y(2)*y(3);-0.28*y(3)+0.7*y(1)*y(3)-0.1*y(2)*y(3)]';
f=inline(ftexto,'t','y');
%Así convierto el texto en una función
%Discretización
h=0.1;
N=round((tn-t0)/h);
%Creación del vector t
t=t0:h:tn;
%Creación el vector y de ceros
y=zeros(ne,N+1);
%Se adjudica el primer valor de y
y(:,1)=y0;
%Aplicación del esquema numérico
for i=1:N;
K1=f(t(i),y(:,i));
K2=f(t(i)+0.5*h,y(:,i)+0.5*K1*h);
K3=f(t(i)+0.5*h,y(:,i)+0.5*K2*h);
K4=f(t(i)+h,y(:,i)+K3*h);
y(:,i+1)=y(:,i)+h/6*(K1+2*K2+2*K3+K4);
end
%Gráfica conjunta de las tres funciones
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
legend('Número de presas','Número de depredadores de la primera especie','Número de depredadores de la segunda especie','Location','Best')
title('Ecosistema')
xlabel('Tiempo')
ylabel('Millones de individuos')
hold off
%Máximos y mínimos
maximopresa=max(y(1,:))
maximopredador1=max(y(2,:))
maximopredador2=max(y(3,:))
minimopresa=min(y(1,:))
minimopredador1=min(y(2,:))
minimopredador2=min(y(3,:))
posmaxpresa=ceil(find(y==max(y(1,:)))/3);
posminpresa=ceil(find(y==min(y(1,:)))/3);
posmaxpred1=ceil(find(y==max(y(2,:)))/3);
posminpred1=ceil(find(y==min(y(2,:)))/3);
posmaxpred2=ceil(find(y==max(y(3,:)))/3);
posminpred2=ceil(find(y==min(y(3,:)))/3);
tmaxpresa=t(posmaxpresa)
tminpresa=t(posminpresa)
tmaxpred1=t(posmaxpred1)
tminpred1=t(posminpred1)
tmaxpred2=t(posmaxpred2)
tminpred2=t(posminpred2)
%Otra opcion sería:
%[M1,pt]=max(y(1,:));
%tmax1=t(pt)
4.1 . Para [math]p_{0}=3,5 millones[/math] de presas, [math]d_{0}=0,001 millones[/math] de predadores de la primera especie, [math]e_{0}=0,0002 millones[/math] de predadores de la segunda especie
Como se ha mencionado anteriormente, el programa obtiene los valores máximos y mínimos de cada especie, así como el momento en el que se producen. Los datos obtenidos para los valores iniciales [math]p_{0}=3,5 millones[/math] de presas, [math]d_{0}=0,001 millones[/math] de predadores de la primera especie, [math]e_{0}=0,0002 millones[/math] de predadores de la segunda especie son los siguientes:
-Máximo número de presas: 8.1307 millones, que se produce en un tiempo de 2.3 años.
-Mínimo número de presas: 1.2206e-08 millones, que se produce en un tiempo de 13.2 años.
-Máximo número de depredadores de la primera especie: 1.0000e-03 millones, que se produce en un tiempo de 0 años.
-Mínimo número de depredadores de la primera especie: 1.9540e-102 millones, que se produce en un tiempo de 500 años.
-Máximo número de depredadores de la segunda especie: 15.1285 millones, que se produce en un tiempo de 68.7 años.
-Mínimo número de depredadores de la segunda especie: 4.0629e-06 millones, que se produce en un tiempo de 60 años.
La gráfica de las variaciones de las especies en relación con el paso del tiempo, es la que se muestra a continuación:
Estudiando la gráfica y teniendo en cuenta además los valores máximos y mínimos de cada población, así como el instante en el que se producen, se puede ver que tanto las presas como los depredadores de la segunda especie crecerán al mismo ritmo, hasta que las primeras se extingan, momento en el que los depredadores comienzan también a desaparecer. Por otro lado, los depredadores de la primera especie prácticamente no aparecen en el ecosistema, y los pocos que hay se extinguen rápidamente. Según los valores alcanzados, en una situación real se extinguirían prácticamente las tres especies, por lo tanto sería un ecosistema inestable.
