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]x'(t)=ax(t)[/math]

[math]y'(t)=-by(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]x'(t)=ax(t)-[término de interacción][/math]

[math]y'(t)=-by(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]x'(t)=ax(t)-cx(t)y(t)[/math]

[math]y'(t)=-by(t)+dx(t)y(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]X=b/d[/math]:[math]Y=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 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:

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

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