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;
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')
%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.
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
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,00001 millones[/math] de predadores de la primera especie, [math]e_{0}=0,2 millones[/math] de predadores de la segunda especie son los siguientes:
-Máximo número de presas: 4.0175 millones, que se produce en un tiempo de 77.4 años.
-Mínimo número de presas: 1.7460e-04 millones, que se produce en un tiempo de 9.9 años.
-Máximo número de depredadores de la primera especie: 1.0000e-05 millones, que se produce en un tiempo de 0 años.
-Mínimo número de depredadores de la primera especie: 2.3015e-103 millones, que se produce en un tiempo de 500 años.
-Máximo número de depredadores de la segunda especie: 7.7652 millones, que se produce en un tiempo de 232.4 años.
-Mínimo número de depredadores de la segunda especie: 0.0033 millones, que se produce en un tiempo de 32.8 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:
Observando la gráfica, así como los valores máximos y mínimos, se puede concluir que en este caso el ecosistema sí sería estable, puesto que al contrario que en el apartado anterior, ni las presas ni los depredadores de segunda especie se extinguirían. Esto se debe a que en el apartado anterior, tanto el número de depredadores de la primera especie como el de la segunda era muy pequeño, haciendo que las interacciones entre ellos y su tasa de mortalidad fuesen motivo suficiente para extinguir las especies. Sin embargo, en este apartado es únicamente el valor inicial de depredadores de primera especie el pequeño, haciendo que se extingan pronto, pero manteniendo un ecosistema que sigue el modelo depredador-presa de Lotka-Volterra. Es por ello que ambas poblaciones, las presas y los depredadores de segunda especie crecen de forma similar y al morir muchas presas, desciende también el número de depredadores de segunda especie a gran velocidad. Las presas se reducen hasta los 174 individuos, y los depredadores de segunda especie hasta los 3300 individuos.
5 . Comparación de resultados
Se puede ver con facilidad que a pesar de estar tratando el mismo problema, las soluciones varían en función del método que se utilice; no a nivel de resultados numéricos, pues los datos iniciales varían para cada apartado, sino en lo que a exactitud se refiere. Por ello, se va a realizar una comparación entre los siguientes métodos:
5.1 . Euler frente a Euler Modificado
Para comparar los resultados, se ha decidido estudiar los apartados 2.2 y 3.2, es decir, los estudios para un periodo de [math]T=[0,300][/math] años y un paso de [math]h=0.1[/math], aunque sin equiparar el número de individuos inicial de depredadores de la segunda especie, que para Euler Modificado es de [math]d0=1.2 millones[/math] y para Euler de [math]d0=0.2 millones[/math].
Se puede observar que para Euler, el número inicial de depredadores de la segunda especie es tan pequeño que se extingue en un periodo de tiempo muy corto, sin llegar a crecer en ningún momento, convirtiendo el modelo de dos depredadores y una presa en un modelo depredador-presa. Aunque para el programa realizado con Euler Modificado la segunda especie depredadora también se acaba extinguiendo, ésta lo hace en un periodo de tiempo mucho más largo. Por el mismo motivo de la diferencia entre los valores iniciales de los depredadores de la segunda especie, el número de presas para Euler Modificado únicamente alcanza un valor máximo de aproximadamente 15 millones, mientras que para Euler llega a ser de 103.5183 millones, al igual que el número de individuos de la primera especie, que para Euler Modificado no llega a 3 millones, mientras que para Euler alcanza los 33.4513 millones de individuos. En resumen, si la población de depredadores de segunda especie se extingue rápidamente, esto permitirá al resto de especies a crecer a mayor velocidad, puesto que no habrá competencia entre depredadores y las presas contarán con menos individuos que las extingan.
Para ver si la diferencia de métodos ha influido en la variación de los resultados, se ha decidido comparar ambos con un programa para el que los dos tengan los mismos datos iniciales.
Para comparar estos métodos, se ha decidido generar un programa que permitiese calcular ambos a la vez y en el que se obtuviese como solución las gráficas de las trayectorias de cada población, comparadas, para un periodo de [math]T=[0,300][/math] años y un paso de [math]h=0.1[/math]. Se ha generado un programa nuevo dado que para cada apartado anterior los datos del problema eran distintos. Además, el programa compara la evolución de las poblaciones para ambos métodos y obtiene un vector "Variacion" que indica la diferencia entre métodos para cada posición:
%Datos del problema
ne=3;
t0=0;
tn=300;
y0=[3.5;1;1.2];
ftexto='[0.4*yem(1)-0.3*yem(1)*yem(2)-0.35*yem(1)*yem(3);-0.3*yem(2)+0.05*yem(1)*yem(2)-0.1*yem(2)*yem(3);-0.28*yem(3)+0.045*yem(1)*yem(3)-0.1*yem(2)*yem(3)]';
f=inline(ftexto,'t','yem');
gtexto='[0.4*ye(1)-0.3*ye(1)*ye(2)-0.35*ye(1)*ye(3);-0.3*ye(2)+0.05*ye(1)*ye(2)-0.1*ye(2)*ye(3);-0.28*ye(3)+0.045*ye(1)*ye(3)-0.1*ye(2)*ye(3)]';
g=inline(gtexto,'t','ye');
%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
yem=zeros(ne,N+1);
%Se adjudia el primer valor de y
yem(:,1)=y0;
%Se crea el vector y de ceros
ye=zeros(ne,N+1);
%Se adjudica el primer valor de y
ye(:,1)=y0;
%Aplicación del esquema numérico
for i=1:N;
K=f(t(i),yem(:,i));
yem(:,i+1)=yem(:,i)+h*K;
ye(:,i+1)=ye(:,i)+h*g(t(i),ye(:,i));
end
%Dibujo las gráficas de las trayectorias
figure
subplot(1,2,1)
plot(yem(1,:),yem(2,:),'b')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
subplot(1,2,2)
plot(ye(1,:),ye(2,:),'r')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
figure
subplot(1,2,1)
plot(yem(1,:),yem(3,:),'b')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
subplot(1,2,2)
plot(ye(1,:),ye(3,:),'r')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
figure
subplot(1,2,1)
plot(yem(2,:),yem(3,:),'b')
title('Trayectoria')
xlabel('Número de depredadores de primera especie')
ylabel('Número de depredadores de segunda especie')
subplot(1,2,2)
plot(ye(2,:),ye(3,:),'r')
title('Trayectoria')
xlabel('Número de depredadores de primera especie')
ylabel('Número de depredadores de segunda especie')
figure
subplot(1,3,1)
plot(yem(1,:),ye(1,:),'b')
title('Comparación')
xlabel('Número de presas para Euler Modificado')
ylabel('Número de presas para Euler')
subplot(1,3,2)
plot(yem(2,:),ye(2,:),'b')
title('Comparación')
xlabel('Número de depredadores de primera especie para Euler Modificado')
ylabel('Número de depredadores de primera especie para Euler')
subplot(1,3,3)
plot(yem(3,:),ye(3,:),'b')
title('Comparación')
xlabel('Número de depredadores de segunda especie para Euler Modificado')
ylabel('Número de depredadores de segunda especie para Euler')
%Para obtener la variación entre métodos se utilizará:
Variacion=ye-yem
Así, las gráficas de comparación entre las trayectorias para cada método son las siguientes:
Como se puede comprobar, los dos métodos calculan las trayectorias de la misma manera y la diferencia entre ambos es tan pequeña que es inapreciable. Ésto se percibir también si se observa el último gráfico que devuelve el programa, en el que se compara el número de individuos para cada población en ambos métodos. Por último, se ve que el vector "Variación" devuelve ceros como resultado.
Por tanto se puede afirmar que ambos métodos son igual de exactos, lo que se confirma teniendo en cuenta que ambos métodos son de orden similar y por lo que se deduce que la diferencia entre resultados es debida únicamente a la variación a la hora de establecer los datos iniciales del problema.
5.2 . Euler frente a Runge-Kutta
Si comparamos el método de Runge-Kutta con el de Euler, podemos observar que el tamaño de paso tiene una gran relevancia ya que, dependiendo de que método se aplique, se podrá hallar una solución válida para ese paso o no. Con Euler no se podía aplicar el paso [math]h=1[/math], como ya vimos anteriormente, y en el primer caso de Runge-Kutta, con unos valores iniciales de [math][p0,d0,e0]=[3.5,0.001,0.0002] millones[/math], tampoco, puesto que como se puede ver en el gráfico los valores se disparan y no serían aceptables para un estudio, aunque sí se obtendría una solución coherente para ese mismo caso con unos valores iniciales de [math][p0,e0,d0]=[3.5,0.00001,0.2] millones[/math]. Para el paso [math]h=0.1[/math] funcionaría en todos los casos.
FALTA PONER MÁS DIFERENCIAS
Una vez obtenidas las diferencias entre ambos apartados, se pretende estudiar la influencia de que los métodos empleados sean distintos. Aunque en el apartado anterior, el hecho de que el cálculo se hubiese realizado con distintos métodos no influía, esta vez se está comparando el método de Euler con el de Runge-Kutta 4, siendo el primero de orden uno y el segundo de orden cuatro, por lo que probablemente si tenga relevancia y deba comprobarse. Para ello se ha utilizado un programa similar al del apartado anterior, pero esta vez para el método Runge-Kutta 4 en lugar de Euler Modificado.
%Datos del problema
ne=3;
t0=0;
tn=300;
y0=[3.5;1;1.2];
ftexto='[0.4*ye(1)-0.3*ye(1)*ye(2)-0.35*ye(1)*ye(3);-0.3*ye(2)+0.05*ye(1)*ye(2)-0.1*ye(2)*ye(3);-0.28*ye(3)+0.045*ye(1)*ye(3)-0.1*ye(2)*ye(3)]';
f=inline(ftexto,'t','ye');
gtexto='[0.4*yrk(1)-0.3*yrk(1)*yrk(2)-0.35*yrk(1)*yrk(3);-0.3*yrk(2)+0.05*yrk(1)*yrk(2)-0.1*yrk(2)*yrk(3);-0.28*yrk(3)+0.045*yrk(1)*yrk(3)-0.1*yrk(2)*yrk(3)]';
g=inline(gtexto,'t','yrk');
%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;
%Se crea el vector y de ceros
ye=zeros(ne,N+1);
%Se adjudica el primer valor de y
ye(:,1)=y0;
%Creación el vector y de ceros
yrk=zeros(ne,N+1);
%Se adjudica el primer valor de y
yrk(:,1)=y0;
%Aplicación del esquema numérico
for i=1:N;
ye(:,i+1)=ye(:,i)+h*g(t(i),ye(:,i));
K1=g(t(i),yrk(:,i));
K2=g(t(i)+0.5*h,yrk(:,i)+0.5*K1*h);
K3=g(t(i)+0.5*h,yrk(:,i)+0.5*K2*h);
K4=g(t(i)+h,yrk(:,i)+K3*h);
yrk(:,i+1)=yrk(:,i)+h/6*(K1+2*K2+2*K3+K4);
end
%Dibujo las gráficas de las trayectorias
figure
subplot(1,2,1)
plot(ye(1,:),ye(2,:),'b')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
subplot(1,2,2)
plot(yrk(1,:),yrk(2,:),'r')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de primera especie')
figure
subplot(1,2,1)
plot(ye(1,:),ye(3,:),'b')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
subplot(1,2,2)
plot(yrk(1,:),yrk(3,:),'r')
title('Trayectoria')
xlabel('Número de presas')
ylabel('Número de depredadores de segunda especie')
figure
subplot(1,2,1)
plot(ye(2,:),ye(3,:),'b')
title('Trayectoria')
xlabel('Número de depredadores de primera especie')
ylabel('Número de depredadores de segunda especie')
subplot(1,2,2)
plot(yrk(2,:),yrk(3,:),'r')
title('Trayectoria')
xlabel('Número de depredadores de primera especie')
ylabel('Número de depredadores de segunda especie')
figure
subplot(1,3,1)
plot(ye(1,:),yrk(1,:),'b')
title('Comparación')
xlabel('Número de presas para Euler Modificado')
ylabel('Número de presas para Euler')
subplot(1,3,2)
plot(ye(2,:),yrk(2,:),'b')
title('Comparación')
xlabel('Número de depredadores de primera especie para Euler Modificado')
ylabel('Número de depredadores de primera especie para Euler')
subplot(1,3,3)
plot(ye(3,:),yrk(3,:),'b')
title('Comparación')
xlabel('Número de depredadores de segunda especie para Euler Modificado')
ylabel('Número de depredadores de segunda especie para Euler')
%Para obtener la variación entre métodos se utilizará:
Variacion=yrk-ye
Así, las gráficas de comparación entre las trayectorias para cada método son las siguientes:
Como se puede comprobar, esta vez las trayectorias obtenidas por cada método, a pesar de haber utilizado los mismos datos iniciales, no son iguales. Como se anticipaba anteriormente, esto es debido a que el orden de los métodos utilizados es distinto.
Aunque no es necesario, puesto que las gráficas de las trayectorias son lo suficientemente determinantes, se ha obtenido una gráfica que relaciona el número de individuos de la misma especie para ambos métodos. Como era de esperar, ésta relación no sigue una recta en la que ambos valores sean iguales a lo largo del tiempo. Por último, esta vez el vector "Variacion" no es un vector nulo, sino que presenta la diferencia de número de individuos para las distintas especies entre ambos métodos, en la misma posición.
