Revisión del 14:16 27 feb 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 natalidad de la presa;
• DXpXd es la tasa de natalidad de la presa;

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

clear all
t0=0; tN=100;
y0=[3.5;1;1.2];
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN;
A1=0.4; A2=0.3; A3=0.35; B1=0.3;
B2=0.05; D=0.1; C1=0.28;  C2=0.045;
yy=y0;
y(:,1)=yy;
A=[A1 0 0;0 -A2 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

Error = abs( exact solution - numerical approximation)

Numerical approximation of the solution

Error = abs( exact solution - numerical approximation)

En la primera gráfica podemos apreciar la extinción de una de las especies depredadoras 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 incial 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.

3 Explicación

3.1 Solución mediante el método Euler modificado

clear all
t0=0;tN=300;
h=1;
N=(tN-t0)/h;
t=t0:h:tN;
 
%VALORES
A1=0.5; A2=0.25; A3=0.3; B1=0.4;
B2=0.07; D=0.05; C1=0.38;  C2=0.045;
 
%DATOS
y0=[3.5;1;0.02];
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)];
 
 
    K2=[A1*(yy(1)+h*K1(1))-A2*(yy(1)+h*K1(1))*(yy(2)+h*K1(2))-A3*(yy(1)+h*K1(1))*(yy(3)+h*K1(3));
        -B1*(yy(2)+h*K1(2))+B2*(yy(1)+h*K1(1))*(yy(2)+h*K1(2))-D*(yy(2)+h*K1(2))*(yy(3)+h*K1(3));
        -C1*(yy(3)+h*K1(3))+C2*(yy(1)+h*K1(1))*(yy(3)+h*K1(3))-D*(yy(3)+h*K1(3))*(yy(3)+h*K1(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')

Numerical approximation of the solution

Error = abs( exact solution - numerical approximation)

Numerical approximation of the solution

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

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)];
 
 
    K2=[A1*(yy(1)+(h/2)*K1(1))-A2*(yy(1)+(h/2)*K1(1))*(yy(2)+(h/2)*K1(1))-A3*(yy(1)+(h/2)*K1(1))*(yy(3)+(h/2)*K1(1));
        -B1*(yy(2)+(h/2)*K1(2))+B2*(yy(1)+(h/2)*K1(2))*(yy(2)+(h/2)*K1(2))-D*(yy(2)+(h/2)*K1(2))*(yy(3)+(h/2)*K1(2));
        -C1*(yy(3)+(h/2)*K1(3))+C2*(yy(1)+(h/2)*K1(3))*(yy(3)+(h/2)*K1(3))-D*(yy(3)+(h/2)*K1(3))*(yy(3)+(h/2)*K1(3))];
    
    K3=[A1*(yy(1)+(h/2)*K2(1))-A2*(yy(1)+(h/2)*K2(1))*(yy(2)+(h/2)*K2(1))-A3*(yy(1)+(h/2)*K2(1))*(yy(3)+(h/2)*K2(1));
        -B1*(yy(2)+(h/2)*K2(2))+B2*(yy(1)+(h/2)*K2(2))*(yy(2)+(h/2)*K2(2))-D*(yy(2)+(h/2)*K2(2))*(yy(3)+(h/2)*K2(2));
        -C1*(yy(3)+(h/2)*K2(3))+C2*(yy(1)+(h/2)*K2(3))*(yy(3)+(h/2)*K2(3))-D*(yy(3)+(h/2)*K2(3))*(yy(3)+(h/2)*K2(3))];
    
    K4=[A1*(yy(1)+h*K3(1))-A2*(yy(1)+h*K3(1))*(yy(2)+h*K3(1))-A3*(yy(1)+h*K3(1))*(yy(3)+h*K3(1));
        -B1*(yy(2)+h*K3(2))+B2*(yy(1)+h*K3(2))*(yy(2)+h*K3(2))-D*(yy(2)+h*K3(2))*(yy(3)+h*K3(2));
        -C1*(yy(3)+h*K3(3))+C2*(yy(1)+h*K3(3))*(yy(3)+h*K3(3))-D*(yy(3)+h*K3(3))*(yy(3)+h*K3(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