Modelo predador-presa de A. Lotka y V. Volterra

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. Grupo 6B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Paula De Santos Muñoz, María del Mar García Reinaldos, Elisabet Sánchez López, Ana Santos Martín, Isabel Sastre Furones, Ángel Antonio Villa Figueroa
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El modelo predador-presa fue propuesto de forma independiente a comienzos del siglo pasado por A. Lotka y V. Volterra. Es un modelo mátematico para la dinámica de poblaciones de especies competidoras (predadores y presas), que consiste en un sistema de ecuaciones diferenciales de primer orden no lineales correspondientes al cambio de las poblaciones de presa y predador por unidad de tiempo.

1 Interpretación del problema de valor inicial

El problema planteado trabaja la existencia de dos especies de presas distintas que comparten una población de depredadores. Llamaremos:

• x₁: especie presa 1
• x₂: especie presa 2
• x₃: población depredadora

Nuestro trabajo se basa en un estudio del crecimiento de las presas y del depredador en el tiempo, observando las interacciones entre ellas y teniendo en cuenta que:

1. En condiciones iniciales, x₁ cuenta con una población de p₀, x₂ con una población de q₀ y finalmente x₃ con una población d₀.
2. Las poblaciones evolucionan interactuando entre ellas, sin tener en cuenta agentes externos.
3. Las presas cuenta con una cantidad suficiente de comida para ser alimentada.
4. La alimentación de los depredadores depende completamente de la presa.

Con estos datos podemos definir las ecuaciones de nuestro sistema. La primera ecuación nos define la evolución de la presa en el tiempo cuando está interactuando con la especie depredadora donde:

• A₁x₁ es la tasa de natalidad de la presa 1, siguiendo una ley exponencial.
• A₂x₁x₃ es la tasa de mortalidad de la presa 1 afectada por los depredadores.

[math]\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\[/math]

De forma análoga a la primera ecuación, podemos definir la segunda.

[math]\frac{\mathrm{d} x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\[/math]

En la tercera ecuación definimos la evolución de x₃, definiendo la constante C₁ multiplicada por (-1), lo que nos hace ver la extinción de los depredadores ante la inexistencia de presas y el crecimiento de los mismos gracias a la alimentación aportada por las presas.

[math]\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\[/math]

Con estás tres ecuaciones queda definido el sistema de ecuaciones de Lotka y Volterra a resolver como un problema de valor inicial.

[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]

2 Resolución por el método de Euler

2.1 Programación del método numérico

El programa utilizado para resolver el problema de valor inicial mediante el método de Euler en un intervalo de tiempo tє[0,300], para una longitud de paso h=0.1, constantes A₁=0.35, A₂=0.6, B₁=0.3, B₂=0.5, C₁=0.37, C₂=0.04 y C₃=0.035, y valores iniciales de p₀=0.8 millones de presas de un tipo, q₀=2.4 millones de presas de otro tipo y d₀=0.2 millones de predadores es el siguiente:

%Datos del problema
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;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=0.8;
q0=2.4;
d0=0.2;
y0=[p0 q0 d0]';
%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 y matriz de soluciones aproximadas
t=t0:h:tf;
y=zeros(3,N+1);
%Inicialización
y(:,1)=y0;
yy=y0;
for n=1:N
    yy=yy+h*(A*yy+B*yy(1)*yy(3)+C*yy(2)*yy(3));
    y(:,n+1)=yy;
end
figure(1)
plot(t,y)
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
figure(2)
subplot(1,2,1)
plot(y(1,:),y(3,:))
title('Población del primer tipo de presa en  función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y(2,:),y(3,:))
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight

2.2 Gráficas

2.2.1 Considerando la evolución de las especies en el tiempo

De la aplicación del método numérico se obtienen las siguientes gráficas, de las que sacaremos conclusiones posteriormente.

En esta primera gráfica, obtenida como resultado de utilizar como longitud de paso h=0.1, se muestran las tres funciones equivalentes a la evolución en el tiempo de los dos tipos de presas y del predador. Cada curva presenta unos máximos periódicos desplazados los unos respecto de los otros, lo cual se debe al desfase con el que evolucionan las distintas especies: primero aumentan las poblaciones de ambas presas (creciendo más aquella de la que inicialmente hay más ejemplares) y después comienza a aumentar la población de la especie predadora, pues puede alimentarse fácilmente. Ésta última alcanzará su primer máximo cuando las presas disminuyen, tras lo cual disminuye también el número de predadores al escasear el alimento. Este proceso se dará periódicamente a lo largo de los 300 años. Es notable la diferencia en la evolución de ambos tipos de presas: mientras que uno tiende a aumentar, el otro disminuye.

Si ahora utilizamos como longitud de paso h=1, obtenemos esta segunda gráfica de la cual podemos concluir que, debido precisamente al aumento de h, el método de Euler pierde su validez.

2.2.2 Considerando la interacción entre especies

Dado que hemos comprobado anteriormente que el método de Euler con longitud de paso h=1 no es efectivo, evaluaremos la interacción entre especies utilizando h=0.1 como longitud de paso. Las gráficas resultantes son las siguientes:

2.3 Interpretación de resultados

Al aumentar la longitud paso de h=0.1 a h=1, lo que observamos es una gran pérdida de precisión en los resultados, ya que el valor de la aproximación calculada por el método de Euler se aleja demasiado de la solución exacta del problema de valor inicial. Esto se debe a que este método numérico se basa en aproximaciones mediante la tangente a una determinada función. En nuestro caso, las diferencias entre los máximos y mínimos de las funciones que describen la evolución de las especies son considerables (son funciones oscilatorias) y se producen en intervalos de tiempo cortos, por lo que, al utilizar una longitud de paso grande (como puede ser h=1), la inexactitud de la aproximación calculada resulta inaceptable. Por el contrario, si disminuimos la longitud de paso, el valor de la aproximación se considera válido. Esto se explica porque el método de Euler alcanza mayor nivel de convergencia cuanto menor es el valor de h; sin embargo, hay que tener en cuenta que, si disminuimos mucho este valor, el programa puede llegar a fallar.

3 Resolución por el método de Euler modificado

3.1 Programación del método numérico

A continuación, se resuelve el problema de valor inicial mediante el método de Euler modificado. Como datos tenemos una longitud de paso h=0.1, constantes A₁=0.35, A₂=0.6, B₁=0.3, B₂=0.5, C₁=0.37, C₂=0.04 y C₃=0.035, y valores iniciales de p₀=2 millones de presas de un tipo, q₀=1.4 millones de presas de otro tipo y d₀=1 millones de predadores. Primero tomaremos un intervalo de tiempo tє[0,100]. El programa se describe a continuación:

%Datos del problema
t0=0;
tf=100;
A1=0.35; A2=0.6;
B1=0.3; B2=0.5;
C1=0.37; C2=0.04; C3=0.035;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=2;
q0=1.4;
d0=1;
y0=[p0 q0 d0]';
%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;
%Inicialización
y=y0;
for n=1:N
    k1=(A*y+y(1)*y(3)*B+y(2)*y(3)*C);
    k2=(A*(y+h*k1)+(y(1)+h*k1(1))*(y(3)+h*k1(3))*B+(y(2)+h*k1(2))*(y(3)+h*k1(3))*C);
    y=y+h*(k1+k2)/2;
    y1(n+1)=y(1);
    y2(n+1)=y(2);
    y3(n+1)=y(3);
end
%Gráfica poblaciones-tiempo
figure(1)
hold on
plot(t,y1,'b')
plot(t,y2,'g')
plot(t,y3,'r')
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
hold off
%Gráficas presa-predador
figure(2)
subplot(1,2,1)
plot(y1,y3,'b')
title('Población del primer tipo de presa en  función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y2,y3,'r')
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight

Hacemos lo mismo para el segundo intervalo de tiempo tє[0,300]. El único cambio que sufre el programa es tf=300.

3.2 Gráficas

3.2.1 Considerando la evolución de las especies en el tiempo

De la aplicación del método numérico para un tiempo tє[0,100], tenemos la siguiente gráfica:

De ella podemos concluir que, durante los primeros 100 años, la evolución de ambos tipos de presas sigue direcciones opuestas: mientras la curva que describe la evolución de la presa de la que hay más ejemplares inicialmente tiende a que sus máximos sean cada vez más pequeños, los máximos de la otra son cada vez mayores. Se observa que, al finalizar el período, el número de ejemplares de ambas presas es similar. Por otro lado, la evolución de la cantidad de predadores en el tiempo es similar a la obtenida al aplicar el método de Euler, aumentando el número de ejemplares tras producirse un pico en las cantidades de presas. Cabe destacar, sin embargo, que con el método de Euler modificado obtenemos gráficas más suaves, lo cual permite aumentar un poco la longitud de paso h sin acumular demasiado error.

Si ahora aumentamos nuestro intervalo de tiempo hasta tє[0,300], podemos comprobar en la gráfica siguiente que la evolución del número de presas sigue la tendencia indicada durante los primeros 100 años: la especie predadora mantiene su evolución oscilatoria con máximos constantes, las presas cuyos máximos disminuían progresivamente casi se han extinguido transcurridos los 300 años y las otras alcanzan sus cifras máximas al final del período. Esto se debe a que los predadores se centran en un determinado tipo de presa, con lo cual la otra se encuentra menos amenazada.

3.2.2 Considerando la interacción entre especies

Para evaluar la interacción entre especies tendremos en cuenta las gráficas obtenidas para un tiempo tє[0,300], ya que se observa mejor la variación de las poblaciones de presas respecto de la de predadores.

Estas gráficas nos permiten concluir que los puntos en los que el ecosistema sería estable son los orígenes (0,0), los cuales despreciamos (pues si lo interpretamos en términos de poblaciones, no hay especies interactuando), y los puntos donde se encuentran los centros de las espirales. Como podemos comprobar, ninguna de las dos curvas converge al mismo centro, así que no podemos asegurar que los ecosistemas sean estables.

4 Resolución por el método trapezoidal

El método trapezoidal es un método implícito que, en nuestro caso, requiere la resolución de una ecuación no lineal. Al no ser posible su resolución, se descarta el método y no se tendrá en cuenta a la hora de comparar resultados.

5 Resolución por el método de Runge-Kutta

5.1 Programación del método numérico

A continuación, se resuelve el problema de valor inicial mediante el método de Runge- Kutta de cuarto orden. Como datos tenemos una longitud de paso h=0.1, constantes A₁=0.4, A₂=0.7, B₁=0.2, B₂=0.4, C₁=0.37, C₂=0.04 y C₃=0.03, y valores iniciales de p₀=3.5 millones de presas de un tipo, q₀=0.2 millones de presas de otro tipo y d₀=0.4 millones de predadores. En este caso tomaremos un intervalo de tiempo tє[0,500]. El programa de obtención del método numérico se especifica a continuación:

t0=0; tN=500; %Extremos del intervalo
y0=[3.5;0.2;0.4]; %Valores de población iniciales
y=y0;
P1(1)=y(1); %Valor inicial de P1=presas tipo I
P2(1)=y(1); %Valor inicial de P2=presas tipo II
D(1)=y(1); %Valor inicial de D=depredadores
h=0.1; %Amplitud de los subintervalos
N=(tN-t0)/h; %Número de subintervalos
A=[ 0.4 0 0;0 0.2 0;0 0 -0.37]; %Parámetros
%Algoritmo
for k=1:N
  k1=(A*y+y(1)*y(3)*[-0.7;0;0.04]+y(2)*y(3)*[0;-0.4;0.03]);
  k2=(A*(y+(h/2)*k1)+(y(1)+(h/2)*k1(1))*(y(3)+(h/2)*k1(3))*[-0.7;0;0.04]+...
  (y(2)+(h/2)*k1(2))*(y(3)+(h/2)*k1(3))*[0;-0.4;0.03]);
  k3=(A*(y+(h/2)*k2)+(y(1)+(h/2)*k2(1))*(y(3)+(h/2)*k2(3))*[-0.7;0;0.04]+...
  (y(2)+(h/2)*k2(2))*(y(3)+(h/2)*k2(3))*[0;-0.4;0.03]);
  k4=(A*(y+h*k3)+(y(1)+h*k3(1))*(y(3)+h*k3(3))*[-0.7;0;0.04]+...
  (y(2)+h*k3(2))*(y(3)+h*k3(3))*[0;-0.4;0.03]);
  y=y+h/6*(k1+2*k2+2*k3+k4);
  P1(k+1)=y(1);
  P2(k+1)=y(2);
  D(k+1)=y(3);
end

t=t0:h:tN;

%CÁLCULO DE INSTANTES MÁXIMOS Y MÍNIMOS

[maximopresaI,x1]= max(P1)
t1=t(x1)
[minimopresaI,x2]= min(P1)
t2=t(x2)
[maximopresaII,x3]= max(P2)
t3=t(x3)
[minimopresaII,x4]= min(P2)
t4=t(x4)
[maximopredadores,x5]= max(D)
t5=t(x5)
[minimopredadores,x6]= min(D)
t6=t(x6)


%GRÁFICA PRESA TIPO I/PRESA TIPO II/DEPREDADORES

hold on
figure(1)
plot(t,P1,'c')
plot(t,P2,'r')
plot(t,D,'g')
title(' Población de presas (tipos I y II) y depredadores')
xlabel('t(años)')
ylabel('población(millones de individuos)')
hold off

5.2 Gráficas: Evolución de las especies en el tiempo

5.3 Interpretación de resultados

En este gráfico se observan las variaciones poblacionales a lo largo del tiempo de las tres especies; la presa tipo I en color cian, la presa de tipo II en color rojo y los predadores se muestran en verde. Interpretamos que la población de predadores empieza a aumentar cuando la población de presas es máxima, es decir cuando tienen mayor cantidad de comida. Este crecimiento se mantiene hasta que el nivel de presas alcanza un punto crítico, donde deja de haber suficientes presas para sostener la población de predadores, lo que genera una disminución en el número de dichos depredadores, que a su vez da lugar a un aumento en el número de presas.

Podemos observar dos comportamientos distintos según el tipo de presa al que nos referimos:

• En la presa de tipo II observamos a simple vista que la población inicial es muy reducida, incluso menor que la de depredadores. Esto provoca que a pesar de que este tipo de presa no sea el único del que se alimenta el depredador, esta especie termina por extinguirse, como bien se aprecia en la gráfica.
• En la presa de tipo I, vemos que a parte de que su población inicial es muchísimo mayor que la de depredadores, su velocidad de crecimiento también es mayor, y esto asegura la supervivencia de ambas especies.

Como resultado conluimos que tanto la cantidad de depredadores como la de presas del tipo que no se han extinguido mantienen sus ciclos de crecemiento y decrecimiento por lo que ambas se mantienen constantes en lo que se refiere a máximos y mínimos.

Dichos máximos y mínimos se adjuntan en la siguiente tabla, con los correspondientes instantes en los que tienen lugar:

6 Conclusiones

Tras analizar los resultados obtenidos aplicando los diferentes métodos numéricos, llegamos a las siguientes conclusiones:

• El método de Euler es un método de orden 1, por lo que nos ofrecerá peores resultados que los métodos de órdenes superiores, como son el método de Euler modificado (Runge-Kutta de orden 2) y el Runge-Kutta de orden 4. Esto se comprueba al comparar los distintos resultados obtenidos.
• En comparación con los otros métodos, es preferible el método de Runge-Kutta, pues aporta soluciones más precisas. Esto se debe a que permite buenas aproximaciones aumentando el valor de la longitud de paso h, ya que no se acumula el error del cálculo mediante tangentes que se desvían de la solución exacta.