Ver el código de «Modelo Predador-Presa de Lokta-Volterra»

{{Trabajo|Modelo Depredador-Presa de Lotka-Volterra|[[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:Trabajos 2012-13|2012-13]]}}
Las ecuaciones de Lokta-Volterra tratan de modelizar la evolución de dos poblaciones (depredador y presa) en el tiempo.
== 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>

== 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.
=== 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®:


{{matlab|codigo=
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: 

[[Archivo:Euler1.jpg|marco|centro]][[Archivo:Euler2.jpg|marco|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.

=== 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®:
{{matlab|codigo=
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:
[[Archivo:RungeKutta1.jpg|marco|centro]]
[[Archivo:RungeKutta2.jpg|marco|centro]]

=== 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.

== 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:

<gallery>
Archivo:1500R1.jpg|Diagrama Población-Tiempo para una <math>R_0</math> de 1500
Archivo:1500R2.jpg|Diagrama <math>R-F</math> para una <math>R_0</math> de 1500
</gallery>
<gallery>
Archivo:1000R1.jpg|Diagrama Población-Tiempo para una <math>R_0</math> de 1000
Archivo:1000R2.jpg|Diagrama <math>R-F</math> para una <math>R_0</math> de 1000
</gallery>
<gallery>
Archivo:250R1.jpg|Diagrama Población-Tiempo para una <math>R_0</math> de 250
Archivo:250R2.jpg|Diagrama <math>R-F</math> para una <math>R_0</math> de 250
</gallery>

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>:
{{matlab|codigo=
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
}}



[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:Trabajos 2012-13]]