Diferencia entre revisiones de «Modelo Predador-Presa de Lokta-Volterra»

Revisión del 21:58 23 feb 2013

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

Estas soluciones se han realizado tomando un número de iteraciones relativamente pequeño. Por eso, el sistema de la solución no es cíclico. Para valores mayores de N, es decir, trabajando con intervalos h más pequeño, la solución del sistema por el método de Euler se aproxima más a la solución real. Una mejor aproximación a esta solución real se consigue utilizando el Método de Runge-Kutta. El aspecto de la solución por el método de Euler para valores mayores de N se asemeja paulatinamente a la solución por el método de Runge-Kutta que se detalla a continuación.

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