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

En este apartado vamos a resolver el PVI mediante el método de Euler modificado. Como datos tenemoscuna 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 tomamos un intervalo de tiempo t[0,100]. El programa se describe a continuación:

%Datos del problema
t0=0;
tN=100;
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
y0=[2 1.4 1]';
%Matrices de coeficientes del sistema
A=[0.35 0 0;0 0.3 0;0 0 -0.37];
B=[-0.6 0 0.04]';
C=[0 -0.5 0.035]';
t=t0:h:tN;
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,'y')
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

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.