Modelo Depredador-Presa de Lotka-Volterra (grupo 16)

Archivo:LotkaVolterra.jpg

izquierda: Alfred James Lotka; derecha: Vito Volterra

El modelo depredador-presa de Lotka-Volterra es un sistema formado por un par de ecuaciones diferenciales de primer orden no lineales que modeliza el crecimiento de dos poblaciones que interactúan (depredador y presa).

1 Introducción

1.1 Vito Volterra

Vito Volterra (1860-1940), físico y matemático italiano fue catedrático de la universidad de Roma y senador. Su oposición al fascismo y su origen judío significaron la expulsión de su cátedra y de las sociedades científicas italianas. Exiliado en Francia hasta 1939, impartió cursos en distintos países, entre ellos España. Volterra desarrolló la solución a las ecuaciones integrales de límites variables que lleva su nombre y tras la primera guerra mundial, en la que se alistó en el cuerpo de ingenieros, se interesó por la aplicación matemática en la biología, extendiendo y desarrollando la obra del matemático belga Pierre François Verhulst, uno de los “padres” de la ecuación logística.

1.2 Alfred James Lotka

Alfred James Lotka (1880-1949), químico, demógrafo y matemático norteamericano de origen ucraniano escribió un libro de Biología teórica y varios artículos sobre procesos oscilantes en Química, en donde de manera independiente a Volterra trabajó con la misma ecuación logística de Verhulst pero con el fin de describir una reacción química en la cual las concentraciones oscilan y estableció el modelo que hoy se conoce con el nombre de ambos Lotka-Volterra y que representa aún la base de los estudios teóricos acerca de la dinámica de poblaciones y otros modelos matemáticos en campos tan diversos como la economía, interdependencia compleja, sostenibilidad, tratamiento de plagas, etc.

1.3 Explicación de las ecuaciones

Para iniciar su investigación matemática establecieron los siguientes parámetros:

• La especie depredadora se alimenta solo de la especie presa, mientras que ésta se nutre de un recurso que se encuentra en el hábitat en grandes cantidades.
• Las dos poblaciones eran homogéneas, es decir, los parámetros de edad y sexo no cuentan.
• Las características son las mismas en todo el hábitat.
• La probabilidad de interacción entre ambas especies es la misma.

Por lo que se concluye que solo existen dos variables: el tamaño poblacional de la especie depredadora y el de la especie presa y que éstos dependen exclusivamente del tiempo.

Por tanto, si no existiesen depredadores, la población de presas crecería malthusianamente; mientras que si no hubiese presas, la especie depredadora decrecería, también siguiendo un modelo malthusiano, es decir:

[math]R'(t)=aR(t)[/math]

[math]F'(t)=-bF(t)[/math]

Ahora bien, como la interacción beneficia a la especie depredadora y perjudica a la presa, él supuso que seria necesario modificar a los depredadores en un término que diera cuenta del perjuicio para una y del beneficio para la otra, lo que tendría que ser:

[math]R'(t)=aR(t)-[término de interacción][/math]

[math]F'(t)=-bF(t)+[término de interacción][/math]

Luego el problema radica en encontrar una forma analítica para cada término que aparece entre corchetes y, basándose en el argumento de que cuantos más encuentros por unidad de tiempo haya entre presas y depredadores, mayor ha de ser el perjuicio de unos y el beneficio de otros; luego matemáticamente queda:

[math]R'(t)=aR(t)-cR(t)F(t)[/math]

[math]F'(t)=-bF(t)+dR(t)F(t)[/math]

1.3.1 Significado de los parámetros

• a: es la tasa instantánea de aumento de conejos en ausencia de zorros
• b: es la tasa instantánea de disminución de zorros en el caso de ausencia de conejos.
• c: mide la susceptibilidad de los conejos a ser cazados.
• d: mide la capacidad de depredación de los zorros

En resumen, los parámetros a y b representan la razón de nacimiento y muerte entre ambas especies, y las constantes c y d reflejan su interacción, beneficiosa para los zorros(d), perjudicial para los conejos(b). Volterra, basándose en esto, formuló la ley que posteriormente se llamaría "Ley de la periodicidad de Volterra", la cual dice que “el cambio de los tamaños poblacionales de ambas especies son periódicos y el periodo depende solamente de a, b, c y d y del tamaño inicial de las dos especies.”

Posteriormente, formuló otras dos leyes relacionadas con su modelo presa-depredador, las cuales son:

• Su Ley de Conservación de los Promedios. Según ésta ley los promedios de los tamaños poblacionales de la especie presa(x) y de la depredadora (y) son independientes de su tamaño inicial. De la que destacamos las siguientes relaciones:

[math]R=b/d[/math]:[math]F=a/c[/math]

• Ley de la perturbación de los promedios.En relación con la anterior, dice que si las poblaciones de ambas especies son destruidas a una razón proporcional a su tamaño poblacional, el promedio de las presas aumenta mientras que el de los depredadores disminuye

Ésta ultima ley es conocida como "El principio de Volterra" y tiene importantes implicaciones en el uso d insecticidas que destruyen tanto a los insectos depredadores como a sus insectos presa.

2 Resolución mediante el método de Euler

El programa utilizado para la resolución mediante el método de Euler en un intervalo de tiempo tЄ[0,100], teniendo en cuenta que los valores que adoptan los parámetros son a = 0,4, b = 0,37, c = 0,3 y d = 0,05, con una población inicial de 3000 conejos y 1000 zorros, es el siguiente:

clear all
%Definiciones
t0=0; tn=100; % Extremos del intervalo temporal a estudiar
y0=[3;1]; %Valores iniciales
N=10000; h=(tn-t0)/N; %Número y amplitud de subintervalos
y=y0; %Asignacion de valores iniciales
R(1)=y(1); %Iniciamos el vector de población R(conejos)
F(1)=y(2); %Iniciamos el vector de población F(zorros)
A=[0.4 0;0 -0.37]; %coeficientes de la parte lineal del sistema
%Algoritmo
for i=1:N
    y=y+h*(A*y+y(1)*y(2)*[-0.3;0.05]);  %Aplicación del método de Euler 
    R(i+1)=y(1);
    F(i+1)=y(2);
end
t=t0:h:tn; %Vector de subintervalos de tiempo
%Dibujo de gráficas
figure(1)
hold on
plot(t,R,'g') %Crecimiento de población R(conejos) en función del tiempo
plot(t,F,'r') %Crecimiento de población F(zorros) en función del tiempo
title('Población de conejos(verde) y zorros(rojo) en función del tiempo');
xlabel('t(años)');
hold off
figure(2)
plot(R,F,'b') %Relación entre ambas poblaciones
title('Población de zorros en función de la de conejos');
xlabel('R(miles de conejos)');
ylabel('F(miles de zorros)');

Obtendremos como resultado las siguientes gráficas:

En esta gráfica se muestran dos funciones que cada una guarda relación (conejos y zorros) respecto del tiempo. Cada curva presenta unos máximos periódicos pero éstos están desplazados los unos respecto de los otros. Este fenómeno es debido a que existe un desfase en el tiempo entre la evolución de ambas especies, ya que primero ha de crecer la especie presa para que la depredadora tenga alimento y pueda llegar a su máximo, haciéndolo más tarde que la especie anterior. Una vez que la especie depredadora alcanza su máximo es cuando la especie presa empieza a descender y por consiguiente el alimento empieza a escasear, lo que produce una disminución de la población de la especie depredadora; aunque dicho suceso se producirá con más lentitud respecto de la presa. Esto se observa en la gráfica, puesto que la pendiente verde es mucho más acentuada que la roja, por lo tanto la especie presa alcanza antes el mínimo.

Esta gráfica en forma de espiral muestra la evolución de los zorros en función de los conejos. Los puntos de equilibrio en los que el sistema sería estable serían (0,0) que lo despreciamos, y el punto donde se situaría el centro de la espiral. Como vemos la curva no converge al mismo por lo que no podríamos asegurar que el ecosistema formado por conejos y zorros fuera estable.

2.1 Variación de poblaciones iniciales

A continuación, pasamos a resolver la ecuación diferencial variando la población inicial de conejos a 1500, 1000 y 250 respectivamente, manteniendo fijos el resto de parámetros para analizar así el comportamiento de ambas poblaciones en función de esta magnitud. Para ello, basta con sustituir en el programa anterior la primera componente del vector y0. Se obtendrán, en cada caso, las siguientes gráficas:

Para R_o=1500:

Para R_o=1000:

Para R_o=250:

Se observa en estas gráficas, por tanto, que cuanto menor es la población de conejos en el momento inicial, menor es también la de zorros al principio de cada ciclo y más tiempo tardan en empezar a crecer ambas, esto se debe a que la población de conejos es la base del ecosistema conejos-zorros, pues los últimos necesitan a los primeros para poder alimentarse y con ello sobrevivir y aumentar en población.

Sin embargo, también se observa que al disminuir la población de conejos inicial aumenta la población máxima que llega a haber de conejos (máxima abscisa en la gráfica azul) y también la máxima población que llega a existir de zorros en el ecosistema (máxima ordenada en la gráfica azul), esto es causado porque al comenzar con un número más bajo de conejos, el número de zorros que puede sobrevivir es menor por la alimentación y que el número de zorros sea menor provoca, a su vez, que puedan llegar a existir más conejos, pues baja el índicede depredación.

Siguiendo el ciclo, llega un punto donde la población de conejos comienza a ser excesiva en función a la de zorros (máximo en la población de conejos) y esto hace que la de zorros aumente, pues disponen de mucha más comida. Este aumento en la población de zorros, conlleva, por tanto, una disminución en la de conejos, la cual llega un momento en que se convierte en insuficiente para alimentar a los zorros (máxima población de zorros), por lo que el número de zorros disminuye muy acusadamente, volviendo al comienzo del ciclo.

Parece obvio que a medida que disminuye la población inicial de conejos, y manteniendo fijos el resto de parámetros de las ecuaciones, los límites en los que cambia el signo de crecimiento de las poblaciones sean mayores, ya que los distintos periodos que sufre el ecosistema se desarrollan de una forma más lenta, lo cual puede ser observado en la primera gráfica, que muestra que el número de repeticiones cíclicas durante el mismo periodo de tiempo es menor al modificar el parámetro R.

Extinción de población dezorros: Si suponemos, pues, que si lapoblación de zorros disminuye por debajo de los 50 elementos la extinción de laespecie es inevitable debido a factores externos, podemos observar en lasgráficas antes expuestas que la extinción de los zorros es prácticamenteinevitable para una población inicial de 250 conejos, ya que encontramos unarelación inicial entre zorros y conejos (parámetros a y b) de ¼, es decir, quepor cada conejo existen 4 zorros, lo que da lugar a un gran decrecimientoinicial en la población de depredadores por falta de alimento, pues el aumentode población de zorros por comer conejos no es lo suficientemente grande comopara compensar la falta de comida, lo que hace inviable el mantenimiento de lapoblación.

3 Resolución mediante el método Runge-Kutta

A continuación se muestra el programa usado para la resolución numérica de nuestro problema mediante el método de Runge-Kutta de 4º orden para 1000 intervalos.

%Valores iniciales
t0=0; tN=100; %Extremos del intervalo a estudiar
y0=[3,1]'; %Valores de población iniciales
N=1000; h=(tN-t0)/N;%Nº y  amplitud de subintervalos
y=y0;
R(1)=y(1);%Valor inicial de R (población de conejos)
F(1)=y(2);% Valor inicial de F (población de zorros)
a=0.4; b=0.37; c=0.3; d=0.05; %Parámetros
%Algoritmo
for n=1:N
    k1=[a*y(1)-c*y(1)*y(2);-b*y(2)+d*y(2)*y(1)];
    k2=[a*(y(1)+1/2*h*k1(1))-c*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));-b*(y(2)+1/2*h*k1(2))+d*(y(2)+1/2*h*k1(2))*(y(1)+1/2*h*k1(1))];
    k3=[a*(y(1)+1/2*h*k2(1))-c*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));-b*(y(2)+1/2*h*k2(2))+d*(y(2)+1/2*h*k2(2))*(y(1)+1/2*h*k2(1))];
    k4=[a*(y(1)+h*k3(1))-c*(y(1)+h*k3(1))*(y(2)+h*k3(2));-b*(y(2)+h*k3(2))+d*(y(2)+h*k3(2))*(y(1)+h*k3(1))];
    y=y+h/6*(k1+2*k2+2*k3+k4);
    R(n+1)=y(1);
    F(n+1)=y(2);
end
%Dibujo
t=t0:h:tN;
figure(1)
hold on
plot(t,R,'g')%Población de conejos en función del tiempo
plot(t,F,'r')%Población de zorros en función del tiempo
title('Población de conejos(verde) y zorros(rojo) en función del tiempo');
ylabel('t(años)');
hold off
figure(2)
plot(R,F,'b')%Población de zorros en función de la de conejos
title('Población de zorros en función de la de conejos');
xlabel('R(miles de conejos)');
ylabel('F(miles de zorros)');

Obtenemos como resultado las siguientes gráficas:

3.1 Ventajas

En comparación con otros métodos de resolución numérica, el de Runge-Kutta tiene ciertas propiedades que pueden hacerlo preferible frente a los otros.

• Si comparamos con el método de Euler tomando diferentes amplitudes (h) de subintervalo en ambos métodos observamos que el método de Runge-Kutta aporta soluciones más precisas.
  

Esto puede deberse a que con el método de Euler la aproximación se realiza por la recta tangente y por ello, en soluciones cíclicas u oscilatorias, como es nuestro caso, no dan muy buenos resultados a no ser que h sea muy pequeña con respecto al intervalo de estudio. Esto se explica porque el método de Euler alcanza un mayor nivel de convergencia cuanto menor es el valor de h.

Destacar al respecto que podría surgir el problema de que al escoger un valor de h demasiado pequeño el programa no responda, por lo que es necesario encontrar un valor adecuado de h que aporte una buena precisión y que pueda ser soportado.

Sin embargo, esto no ocurre con el método de Runge-Kutta, que permite buenas aproximaciones con mayores valores de h, ya que no depende de tangentes que se desvían de la solución acumulando error que puede llegar a hacer inservible la solución final. Así lo vemos en las siguientes gráficas para las que se ha tomado un valor de h=1:

Con el método de Euler:

Con Runge-Kutta:

• El empleo del método de los trapecios en la resolución de nuestro sistema resulta inviable, ya que se trata de un sistema no lineal.

4 Caso especial: Problema inverso

Supongamos que en vez de tener las condiciones iniciales (t=0) y querer conocer las condiciones dadas a un tiempo determinado, se nos plantease el problema inverso, es decir, dadas las condiciones en t=100 debemos hallar las condiciones iniciales (t=0), en este caso nuestro método sigue siendo útil si le aplicamos los cambios adecuados. Necesitamos hacer que nuestra función avance hacia atrás en el tiempo, en otras palabras, tomar el eje t negativo en dirección derecha. Para ello basta con cambiar de signo todos los términos que dependan de t para que, al invertir el sentido del eje, hagan efecto espejo y nos muestren nuestra función avanzando hacia atrás. Habiendo realizado esto, simplemente hace falta pedirle al programa que nos dé los últimos valores de ambas poblaciones.

clear all
t0=0; tN=100;
y0=[2.2,0.32]';
N=1000; h=(tN-t0)/N;
y=y0;
R(1)=y(1);
F(1)=y(2);
a=-0.4; b=-0.37; c=-0.3; d=-0.05;
for n=1:N
   k1=[a*y(1)-c*y(1)*y(2);-b*y(2)+d*y(2)*y(1)];
   k2=[a*(y(1)+1/2*h*k1(1))-c*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));-b*(y(2)+1/2*h*k1(2))+d*(y(2)+1/2*h*k1(2))*(y(1)+1/2*h*k1(1))];
   k3=[a*(y(1)+1/2*h*k2(1))-c*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));-b*(y(2)+1/2*h*k2(2))+d*(y(2)+1/2*h*k2(2))*(y(1)+1/2*h*k2(1))];
   k4=[a*(y(1)+h*k3(1))-c*(y(1)+h*k3(1))*(y(2)+h*k3(2));-b*(y(2)+h*k3(2))+d*(y(2)+h*k3(2))*(y(1)+h*k3(1))];
   y=y+h/6*(k1+2*k2+2*k3+k4);
   R(n+1)=y(1);
   F(n+1)=y(2);
end
t=t0:h:tN;
figure(1)
hold on
plot(t,R,'g')
plot(t,F,'r')
title('Población de conejos(verde) y zorros(rojo) en función del tiempo');
xlabel('t(años)');
hold off
figure(2)
plot(R,F,'b')
title('Población de zorros en función de la de conejos');
xlabel('R(miles de conejos)');
ylabel('F(miles de zorros)'); 
R(101)
F(101)

Los valores obtenidos serán 1.410,6 para la población de conejos y 521,6 para la de zorros.

Se puede observar en la gráfica 1 que efectivamente no hemos hecho otra cosa que darle la vuelta al eje x (tiempo), ya que es la misma gráfica que en el problema directo (con distintas condiciones iniciales) dada la vuelta.