Modelo Predador-Presa de Lokta-Volterra

Las ecuaciones de Lokta-Volterra tratan de modelizar la evolución de dos poblaciones (depredador y presa) en el tiempo.

1 Planteamiento

El modelo se basa en los siguientes supuestos:

• El crecimiento de la población R(t) (Rabbits, en adelante simplemente R) de presas en ausencia de predadores es proporcional a la población R. Siendo a la diferencia entre las tasas de natalidad y mortalidad de la presa::

[math] \frac{\mathrm{d} R}{\mathrm{d} t}=aR [/math]

• El crecimiento de la población R se ve afectado por la acción de los depredadores. Esta acción es proporcional a la cantidad de interacciones FR entre la población R de presas y la población F(t) (Foxes, en adelante simplemente F) de predadores. El factor de proporcionalidad c indica el grado de efectividad del proceso.:

[math] \frac{\mathrm{d} R}{\mathrm{d} t}=aR-bRT [/math]

• Este mismo factor afecta de forma positiva al crecimiento de la población F de predadores, que se ven beneficiados por la caza. El factor de proporcionalidad d indica el grado en que afecta el éxito en la caza al crecimiento de la población de predadores.:

[math] \frac{\mathrm{d} F}{\mathrm{d} t}=dRT [/math]

• Por último, la competitividad asociada al crecimiento excesivo de la población de predadores F resulta perjudicial para el crecimiento de esta población. El grado en que esto afecta al crecimiento de la población se representa a través del parámetro c.:

[math] \frac{\mathrm{d} F}{\mathrm{d} t}=-cF+dRT [/math]

Por tanto, la expresión final del modelo se puede expresar a través del problema de condiciones iniciales o de Cauchy::

[math] \left\{\begin{matrix}\frac{\mathrm{d} R}{\mathrm{d} t}=aR-bRT\\\frac{\mathrm{d} F}{\mathrm{d} t}=-cF+dRT\\R(t_{0})=R_{0}, F(t_{0})=F_{0}\end{matrix}\right. [/math]

Siendo [math]R_{0},F_{0}[/math] el tamaño de las poblaciones R y F en [math]t=0[/math]

2 Resolución Numérica

Puesto que no se trata de una ecuación no lineal (aunque es posible linealizarla en un entorno de los puntos de equilibrio), vamos a obtener sus soluciones a través de métodos iterativos.

2.1 Resolución por el Método de Euler

Para una población inicial de Presas [math]R_{0}=3000[/math] y de Predadores [math]F_{0}=1000[/math], en un intervalo [math][/math] y para los valores [math]a=0.4,b=0.37,c=0.3,d=0.05[/math] de las constantes, el sistema se resuelve con el siguiente código en MATLAB®:

clear all
N=1000; %Número de iteraciones.
t0=0; tN=100; h=(tN-t0)/N;  %Tiempo inicial, final y espacio entre iteraciones sucesivas.
a=0.4;b=0.37;c=0.3;d=0.05;  %Valor de las constantes del modelo
R(1)=3; F(1)=1;  %Condiciones iniciales

for i = 1:N
  %Cálculo del próximo valor de R, F por el método de Euler.
  R(i+1) = R(i) + h*R(i)*(a-c*F(i));
  F(i+1) = F(i) + h*F(i)*(-b+d*R(i));
end;
T=t0:h:tN;

%Gráfica de la solución.
hold on
figure(1)
plot(T,R,'r')
plot(T,F,'g')
figure(2)
plot(R,F)
hold off

Gráficamente, las soluciones son:

centro

centro

2.2 Resolución por el Método de Runge-Kutta

Para una población inicial de Presas [math]R_{0}=3000[/math] y de Predadores [math]F_{0}=1000[/math], en un intervalo [math][/math] y para los valores [math]a=0.4,b=0.37,c=0.3,d=0.05[/math] de las constantes, el sistema se resuelve con el siguiente código en MATLAB®:

clear all
N=10000; %Número de iteraciones.
t0=0; tN=100; h=(tN-t0)/N; %Tiempo inicial, final y espacio entre iteraciones sucesivas.
a=0.4;b=0.37;c=0.3;d=0.05; %Valor de las constantes del modelo.
R(1)=3; F(1)=1; %Condiciones iniciales.

for i = 1:N       
  %Cálculo de las Ctes k1-k4.    
   k1=[a*R(i)-c*F(i)*R(i), -b*F(i)+d*F(i)*R(i)];    
   k2=[a*(R(i)+h*k1(1)/2)-c*(F(i)+h*k1(2)/2)*(R(i)+h*k1(1)/2), -b*(F(i)+h*k1(2)/2)+d*(F(i)+h*k1(2)/2)*(R(i)+h*k1(1)/2)];    
   k3=[a*(R(i)+h*k2(1)/2)-c*(F(i)+h*k2(2)/2)*(R(i)+h*k2(1)/2), -b*(F(i)+h*k2(2)/2)+d*(F(i)+h*k2(2)/2)*(R(i)+h*k2(1)/2)];    
   k4=[a*(R(i)+h*k3(1))-c*(F(i)+h*k3(2))*(R(i)+h*k3(1)), -b*(F(i)+h*k3(2))+d*(F(i)+h*k3(2))*(R(i)+h*k3(1))];    
  %Cálculo del próximo valor de R, F.    
   R(i+1)=R(i)+h*(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6;    
   F(i+1)=F(i)+h*(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6;
end;
T=t0:h:tN; 

%Gráfica de la solución.
hold on
figure(1)
plot(T,R,'r')
plot(T,F,'g')
figure(2)
plot(R,F)
hold off

Gráficamente, las soluciones son:

centro

centro