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

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.

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.

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.