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

Revisión actual del 17:04 25 jul 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-bRF [/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}=dRF [/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+dRF [/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-bRF\\\frac{\mathrm{d} F}{\mathrm{d} t}=-cF+dRF\\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 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] (a escala 1:1000), 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] (a escala 1:1000), 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

2.3 Resolución por el Método del Trapecio

Este procedimiento se dejó de lado, puesto que se trata de un método implícito y no es sencillo despejar el valor del próximo valor de [math](R,F)[/math] para procesarlo.

3 Interpretación de los resultados

La solución del sistema de Lokta-Volterra nos permite estudiar la evolución de la población de las dos especies para diferentes condiciones iniciales. Por ejemplo, con los mismos parámetros usados anteriormente, pero variando la población inicial de Presas observamos:

• 

  Diagrama Población-Tiempo para una [math]R_0[/math] de 1500

• 

  Diagrama [math]R-F[/math] para una [math]R_0[/math] de 1500

• 

  Diagrama Población-Tiempo para una [math]R_0[/math] de 1000

• 

  Diagrama [math]R-F[/math] para una [math]R_0[/math] de 1000

• 

  Diagrama Población-Tiempo para una [math]R_0[/math] de 250

• 

  Diagrama [math]R-F[/math] para una [math]R_0[/math] de 250

Una de las aplicaciones de este estudio es tratar de predecir la extinción de las especies. Por ejemplo, si suponemos que para una población de Predadores inferior a 50 individuos la extinción se vuelve inevitable debido a otros factores, podríamos determinar que el depredador correría riesgo de extinción para la ya estudiada población inicial de presas de 250, ya que entonces, según nuestra aproximación, la población de predadores llegaría a descender hasta la cifra de 52 individuos.

También se puede utilizar este modelo para calcular la población de predadores y presas en un instante [math]t[/math] anterior a [math]t_0[/math]. Por ejemplo, para calcular la población de presas y predadores en [math]t=0[/math], sabiendo que en [math]t=100[/math], [math]R=2200,F=320[/math], podemos utilizar esta pequeña variación de la anterior aproximación por el método de Runge-Kutta (y análogamente por el método de Euler), obteniendo unos valores de [math]R_0= 1435,F_0=490[/math] en [math] t=0[/math]:

clear all
N=1000; %Número de iteraciones.
t0=100; tN=0; 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(N+1)=2.2; F(N+1)=0.32; %Condiciones iniciales


for i =  N+1:-1:2  
%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; %Vector de valores temporales asociados a cada valor de R,F.


%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