Revisión actual del 18:21 6 mar 2014

Las ecuaciones de Lotka-Volterra o ecuaciones 'predador-presa', son ecuaciones diferenciales ordinarias de primer orden no lineales que se usan para el modelado de dos o más poblaciones que interactúan, presas y depredadores.Las ecuaciones fueron propuestas de forma independiente por Alfred J. Lotka y Vito Volterra en 1900.

Tales ecuaciones se definen como:

[math]\frac{dXp}{dt} = AXp-BXdXp[/math]

[math]\frac{dXd}{dt} = -CXd+DXpXd[/math]

donde

• Xd es el número de los predadores ;
• Xp es el número de sus presas ;
• dXp/dt y dXd/dt representa el crecimiento de las dos poblaciones en el tiempo;
• t representa el tiempo;
• AXp es la tasa de natalidad de la presa;
• BXdXp es la tasa de mortalidad de la presa;
• CXd es la tasa de mortalidad de los depredadores;
• DXpXd es la tasa de natalidad de los depredadores;

1 Explicación de las ecuaciones

1.1 Presa

[math]\frac{dXp}{dt} =AXp-BXdXp [/math]

En el modelo se asume que las presas tienen una cantidad de comida ilimitada, y se reproducen exponencialmente o siguiendo la ley malthusiana. Este crecimiento exponencial está representado en la ecuación por el término AXp. El término de la ecuación BXdXp representa la iteración entre presas y depredadores. Podemos interpretar la ecuación como que el cambio del número de presas por unidad de tiempo viene determinado por su propio crecimiento menos la tasa de encuentros con predadores.

1.2 Depredador

[math]\frac{dXd}{dt} =-CXd+DXpXd [/math]

El término de la ecuación DXpXd representa que la natalidad de los depredadores depende de las relaciones con las presas ya que su alimentación depende completamente de estas. Debido a que la población de presas iría desapareciendo a una razón proporcional a la población presente si no hay alimentos, la tasa de mortalidad de esta especie es -CXd.

De esta manera, interpretamos la ecuación como el crecimiento de los depredadores por la caza de presas menos la muerte de éstos.

Las ecuaciones anteriores hacen referencia al modelo de Lotka-Volterra, pero nuestro caso es mas complejo, puesto que también hay que tener en cuenta la influencia de las especies depredadoras entre sí.

[math]\frac{dX1}{dt} = A1X1-A2X1X2-A3X1X3 (1)[/math]

[math]\frac{dX2}{dt} = -B1X2+B2X1X2-DX2X3 (2)[/math]

[math]\frac{dX3}{dt} = -C1X3+C2X1X3-DX2X3 (3)[/math]

La ecuación (1) representa la variación de la población de las presas. Ésta, será proporcional a la tasa de natalidad y la población existente, pero habrá que restar la interacción con las especies depredadoras X2 y X3, con la tasa de mortalidad respectiva. La (2) y la (3) representan la variación de las especies depredadoras que serán proporcionales a los encuentros que haya con la presa, pero tendrá su propia tasa de mortalidad a la que habrá que restar también los encuentros entre las dos especies depredadoras.

2 Solución mediante el método Euler

En un intervalo de tiempo [0,100], particularizando los parámetros con los valores indicados y con una población inicial de 3.5 millones de presas ,1 millón de un tipo de depredador y 1.2 millones del segundo tipo de depredador, el programa en MATLAB es el siguiente:

clear all

    t0=0; tN=300;

    y0=[3.5;1;0.2];

    h=1;

    N=floor((tN-t0)/h);

    t=t0:h:tN;

    A1=0.5; A2=0.25; A3=0.3; B1=0.4;

    B2=0.07; D=0.05; C1=0.38;  C2=0.045;

    yy=y0;

    y(:,1)=yy;

    A=[A1 0 0;0 -B1 0;0 0 -C1]; %coeficientes de la parte lineal del sistema

     

    for n=1:N

      yy=yy+h*(A*yy+yy(1)*yy(2)*[-A2;B2;0]+yy(1)*yy(3)*[-A3;0;C2]+yy(2)*yy(3)*[0;-D;-D]);

     

      y(:,n+1)=yy;

    end

     

    figure(1)

    hold on

    plot(t,y(1,:),'r')

    plot(t,y(2,:),'g')

    plot(t,y(3,:))

    hold off

     

    figure(2)

    plot(y(1,:),y(2,:),'b') %Relación entre ambas poblaciones

    figure(3)

    plot(y(1,:),y(3,:),'k')

    figure(4)

    plot(y(2,:),y(3,:),'m')

Población de depredadores y presa en función del tiempo

Población del depredador 1 en función de la población de la presa

Población del depredador 2 en función de la población de la presa

Población del depredador 2 en función de la población del depredador 1

Se puede apreciar que dado que el paso "h" es excesivamente grande (en este caso, 1) el método de Euler no proporciona una aproximación adecuada a la solución real, ya que las 3 poblaciones, especialmente la de la presa, sufren una variación relativamente rápida respecto al tiempo. De esta forma, el paso de una unidad no refleja claramente los importantes cambios de pendiente que se producen en la gráfica. El fallo es tal que se da la incoherencia de que aparezcan poblaciones negativas. Probemos entonces con un paso 10 veces menor.

h=0.1:

Población de las 3 especies del ecosistema a lo largo del tiempo

Población del depredador 1 en función de la población de la presa

Población del depredador 2 en función de la población de la presa

Población del depredador 2 en función de la población del depredador 1

En la primera gráfica podemos apreciar la extinción de una de las especies depredadoras (en color azul), algo lógico teniendo en cuenta que se encuentran en un ecosistema en el que conviven dos especies depredadoras y una única presa. Este hecho puede ser debido a que hay un menor número inicial de individuos, tiene una mayor tasa de mortalidad y una menor interacción con las presas, lo que incita a su extinción. Se observa un cierto desfase en las representaciones de las poblaciones, ya que para que los depredadores puedan aumentar su número, debe haber antes un crecimiento en el número de presas. Una vez que la segunda especie depredadora se ha extinguido, se puede apreciar una cierta estabilidad en la correlación entre el número de individuos de las especies supervivientes.

3 Explicación

3.1 Solución mediante el método Euler modificado

A continuación se muestra la resolución numérica mediante el método de Euler modificado , es decir, el método de Runge Kutta de orden 2, para un intervalo de tiempo [0,100].

clear all
t0=0;tN=100;
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN;
 
%VALORES
A1=0.4; A2=0.3; A3=0.35; B1=0.3;
B2=0.05; D=0.1; C1=0.28;  C2=0.045;
 
%DATOS
y0=[3.5;1;1.2];
yy=y0;
y(:,1)=yy;
 
%bucle
 
for n=1:N
    
    
    K1=[A1*yy(1)-A2*yy(1)*yy(2)-A3*yy(1)*yy(3);
        -B1*yy(2)+B2*yy(1)*yy(2)-D*yy(2)*yy(3);
        -C1*yy(3)+C2*yy(1)*yy(3)-D*yy(2)*yy(3)];
 
    yp= yy+h*K1;
    K2=[A1*yp(1)-A2*yp(1)*yp(2)-A3*yp(1)*yp(3);
        -B1*yp(2)+B2*yp(1)*yp(2)-D*yp(2)*yp(3);
        -C1*yp(3)+C2*yp(1)*yp(3)-D*yp(2)*yp(3)];
    
     yy=yy+(h/2)*(K1+K2);
    
    
     y(:,n+1)=yy;
    
end
 
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'g')
plot(t,y(3,:))
hold off
figure(2)
plot(y(1,:),y(2,:),'b') %Relación entre ambas poblaciones
figure(3)
plot(y(1,:),y(3,:),'k')
figure(4)

plot(y(2,:),y(3,:),'m')

Población de depredadores y presa en función del tiempo

Población del depredador 1 en función de la población de la presa

Población del depredador 2 en función de la población de la presa

Población del depredador 2 en función de la población del depredador 1

En el transcurso del tiempo se puede percibir como la población, tanto de presas como de depredadores del primer tipo tienden a estabilizarse, siempre dentro de una estructura cíclica, y aumentando esta estabilidad con la extinción de la especie depredadora 2. Previo a la extinción se observa una clara recesión en el número de presas ya que ambos depredadores conviven durante un cierto periodo de tiempo. La extinción del depredador 2 se aprecia claramente en la 3º gráfica, que muestra la relación de su población con la de la presa. En un principio se nota la estructura cíclica inicial pero pronto comienzan a decrecer hasta alcanzar el 0 tras un número escaso de ciclos (5 aproximadamente).

La resolución numérica para el intervalo de tiempo [0,300]:

clear all
t0=0;tN=300;
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN;
 
%VALORES
A1=0.4; A2=0.3; A3=0.35; B1=0.3;
B2=0.05; D=0.1; C1=0.28;  C2=0.045;
 
%DATOS
y0=[3.5;1;1.2];
yy=y0;
y(:,1)=yy;
 
%bucle
 
for n=1:N
    
    
    K1=[A1*yy(1)-A2*yy(1)*yy(2)-A3*yy(1)*yy(3);
        -B1*yy(2)+B2*yy(1)*yy(2)-D*yy(2)*yy(3);
        -C1*yy(3)+C2*yy(1)*yy(3)-D*yy(2)*yy(3)];
 
    yp= yy+h*K1;
    K2=[A1*yp(1)-A2*yp(1)*yp(2)-A3*yp(1)*yp(3);
        -B1*yp(2)+B2*yp(1)*yp(2)-D*yp(2)*yp(3);
        -C1*yp(3)+C2*yp(1)*yp(3)-D*yp(2)*yp(3)];
    
     yy=yy+(h/2)*(K1+K2);
    
    
     y(:,n+1)=yy;
    
end
 
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'g')
plot(t,y(3,:))
hold off
figure(2)
plot(y(1,:),y(2,:),'b') %Relación entre ambas poblaciones
figure(3)
plot(y(1,:),y(3,:),'k')
figure(4)

plot(y(2,:),y(3,:),'m')

Población de depredadores y presa en función del tiempo

Población del depredador 1 en función de la población de la presa

Población del depredador 2 en función dela población de la presa

Población del depredador 2 en función de la población del depredador 1

Al aumentar el intervalo temporal se observa con mayor claridad la estabilización del ecosistema tras la extinción

3.2 Solución mediante el método Runge-Kutta

Por último, la resolución numérica mediante el método de Runge Kutta de 4ºorden ,en un intervalo de tiempo [0,500] y con valores iniciales 3.5 millones de presas,0.001 millones de depredadores 1 y 0.0002 millones de depredadores 2, es el siguiente:

clear all
t0=0;tN=500;
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN;
 
%VALORES
A1=0.4; A2=0.3; A3=0.4; B1=0.37;
B2=0.05; D=0.1; C1=0.28;  C2=0.07;
 
%DATOS
y0=[3.5;0.001;0.0002];
yy=y0;
y(:,1)=yy;
 
%bucle
 
for n=1:N
    
    
    K1=[A1*yy(1)-A2*yy(1)*yy(2)-A3*yy(1)*yy(3);
        -B1*yy(2)+B2*yy(1)*yy(2)-D*yy(2)*yy(3);
        -C1*yy(3)+C2*yy(1)*yy(3)-D*yy(2)*yy(3)];
 
    yp= yy+h*K1/2;
    K2=[A1*yp(1)-A2*yp(1)*yp(2)-A3*yp(1)*yp(3);
        -B1*yp(2)+B2*yp(1)*yp(2)-D*yp(2)*yp(3);
        -C1*yp(3)+C2*yp(1)*yp(3)-D*yp(2)*yp(3)];
    
     yp= yy+h*K2/2;
    K3=[A1*yp(1)-A2*yp(1)*yp(2)-A3*yp(1)*yp(3);
        -B1*yp(2)+B2*yp(1)*yp(2)-D*yp(2)*yp(3);
        -C1*yp(3)+C2*yp(1)*yp(3)-D*yp(2)*yp(3)];
    
     yp= yy+h*K3;
    K4=[A1*yp(1)-A2*yp(1)*yp(2)-A3*yp(1)*yp(3);
        -B1*yp(2)+B2*yp(1)*yp(2)-D*yp(2)*yp(3);
        -C1*yp(3)+C2*yp(1)*yp(3)-D*yp(2)*yp(3)];
    
     yy=yy+(h/6)*(K1+2*K2+2*K3+K4);
    
    
     y(:,n+1)=yy;
    
end
 
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'g')
plot(t,y(3,:))
hold off
figure(2)
plot(y(1,:),y(2,:),'b') %Relación entre ambas poblaciones
figure(3)
plot(y(1,:),y(3,:),'k')
figure(4)

plot(y(2,:),y(3,:),'m')

Población de las 3 especies del ecosistema en función del tiempo

Población del depredador 1 en función de la de la presa

Población del depredador 2 en función de la de la presa

Población del depredador 2 en función de la del depredador 1

Al ser Runge Kutta un método más preciso, se nota con mayor claridad la estabilización de nuestro sistema y cómo, tras la desaparición de una especie, las otras dos oscilan en torno a unos valores coincidiendo sus máximos y mínimos tras un cierto desfase. También podemos observar cómo los mínimos valores de ambas poblaciones están muy próximos al 0, lo que supondría la extinción de las mismas. En caso de extinguirse ambos depredadores, se produciría un aumento descontrolado de la población de presas. Por el contrario, la extinción de las presas provocaría, una vez transcurrido ese desfase temporal, la extinción de los depredadores supervivientes.

Cambiando las poblaciones iniciales de los depredadores a 0,00001 millones y 0,2 millones respectivamente, obtenemos las gráficas siguientes:

Población de las 3 especies del ecosistema en función del tiempo

Población del depredador 1 en función de la de la presa

Población del depredador 2 en función de la de la presa

Población del depredador 2 en función de la del depredador 1

La variación de las poblaciones iniciales provoca, al haber muchos menos individuos de la población de depredadores 1 (en verde), que éstos se extingan mucho más rápidamente que en el caso anterior. El ecosistema adquiere entonces la estabilidad mucho antes y da lugar a un número mayor de ciclos de crecimiento y decrecimiento de poblaciones de las especies supervivientes. Además, al tener la especie depredadora superviviente una población mayor en el inicio del estudio, en ningún momento se encontrará tan próxima a la desaparición total y absoluta.

3.3 Comparación de métodos

Aproximación por el método de euler

Aproximación por el método de euler modificado

Aproximación por el método de ruge-kutta