Diferencia entre revisiones de «Depredador-Presa ( grupo 12B)»

Revisión del 02:54 5 mar 2014

El contenido de este artículo nos muestra la resolución del modelo matemático de A. Lotka y V. Volterra. Es un sistema formado por ecuaciones diferenciales de primer orden no lineales, que modeliza una lucha constante por la supervivencia de dos especies competidoras que viven en un mismo hábitat siendo de esa manera una la depredadora y la otra su presa. Las ecuaciones fueron propuestas de forma independiente por Alfred Lotka en 1925 y Vito Volterra en 1926. El modelo de Volterra-Lotka se conoce también como modelo depredador- presa.

1 Interpretación

1.1 Hipótesis iniciales

1. La especie depredadora se alimenta solo de la especie presa, mientras que éstas siempre tienen una cantidad suficiente de alimento a su disposición.
2. x1 y x2 son dos especies de presas que conviven en un mismo ecosistema junto con una población x3 de depredadores.
3. En el instante inicial del estudio, la especie x1 tiene una población de p0 individuos, x2 de q0, y x3 de d0.
4. Las iteraciones entre las presas y el depredador es la misma.
5. El sistema de ecuaciones, modeliza las velocidades de crecimiento de las tres poblaciones, teniendo en cuenta únicamente las interacciones entre las tres especies, (no se tienen en cuenta posibles enfermedades o variaciones demográficas repentinas de las distintas especies del ecosistema por causa de desastres naturales).

En un hábitat en el que hay abundante vegetación dejamos evolucionar dos especies una depredadora y otra presa. Si nos centramos en la primera ecuación del sistema, vemos que la velocidad de crecimiento de la especie x1, a partir de ahora G (gazelle), La tasa de natalidad de éstas, debe seguir la ley malthusiana o exponencial sin la presencia de los depredadores, esto es, la tasa de natalidad de la presa A1 x1, donde A1 es una constante positiva, que define la tasa de crecimiento de la población de gacelas en ausencia de depredadores, x3, a partir de ahora L (lion); al existir leones, se añade un término más a la ecuación, definido por otra constante A2, que define la tasa de mortalidad de la presa y en que afecta la interacción entre las gacelas y los leones (las gacelas son comidas por los leones) a la población de gacelas.

[math]\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\ [/math]

Además de esto, y cambiando levemente la ecuación clásica de Lotka-Volterra, contamos con otra especie de tipo presa, x2, a partir de ahora Z (zebra), cuya ecuación (variables y constantes) es análoga a la de las gacelas, siendo su depredador natural el león. En estas dos ecuaciones se observa que no hay competencia entre cebras y gacelas, ya que si ésta existiera, deberíamos añadir a ambas ecuaciones un término más con A3 y B3, para modelizar dicha competencia.

[math]\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\[/math]

Centrándonos en la tercera ecuación, que rige la velocidad de crecimiento de la población de leones, se puede observar que en ausencia de presas, la población de leones decrecería hasta llegar a cero (extinción de la especie), ya que su tasa de crecimiento (C1) está multiplicada por “-1” también siguiendo un modelo malthusiano. El sentido que tiene esto en la naturaleza, en ausencia de presas los leones no podrían comer y la población decrecería, mientras que si existen presas suficientes, la población de leones aumentará. El tiempo negativo no tiene sentido, por tanto, es obvio que t debe ser mayor que cero para nuestras ecuaciones.

[math]\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\[/math]

Combinando estas tres ecuaciones, se obtiene el sistema de Volterra-Lotka para la dinámica de las poblaciones de las dos especies competidoras, que se resolverá como un problema de valor inicial:

[math] \left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0},x_2(0)=q_{0},x_3(0)=d_{0}\end{matrix}\right. [/math]

2 Método de Euler modificado

2.1 Programación método numérico

Utilizando el método de Euler modificado con una longitud de paso de [math]h=0.1[/math], en un intervalo de tiempo tЄ[0,100], con una población inicial de [math]p_0=2[/math] millones de presas de un tipo, [math]q_0=1.4[/math] millones de presas de otro tipo y [math]d_0=1[/math] millon de depredadores y teniendo en cuenta que los valores que adoptan los parámetros son

[math] \left\{\begin{matrix}\ A_1=0.35, A_2=0.6\\ B_1=0.3, B_2=0.5\\ C_1=0.37, C_2=0.04, C_3=0.035\end{matrix}\right. [/math]

El programa resultante es el que se muestra a continuación

clear all
clc
hold off
%Definiciones
t0=0; tn=100; % Extremos del intervalo temporal a estudiar
y0=[2;1.4;1]; %Valores iniciales
h=0.1; N=(tn-t0)/h; %Número y amplitud de subintervalos
y=y0; %Asignacion de valores iniciales
G(1)=y(1); %Iniciamos el vector de población G(gacelas)
Z(1)=y(2);%Iniciamos el vector de población Z(cebras)
L(1)=y(3);%Iniciamos el vector de población L(leones)
A=[0.35 0 0;0 0.3 0;0 0 -0.37]; %coeficientes de la parte lineal del sistema
%Algoritmo
for i=1:N
    k1=(A*y+y(1)*y(3)*[-0.6;0;0.04]+y(2)*y(3)*[0;-0.5;0.035]);
    k2=(A*(y+h*k1)+(y(1)+h*k1(1))*(y(3)+h*k1(3))*[-0.6;0;0.04]+(y(2)+h*k1(2))*(y(3)+h*k1(3))*[0;-0.5;0.035]);
    y=y+(h/2)*(k1+k2);  %Aplicación del método de Euler modificado 
    G(i+1)=y(1);
    Z(i+1)=y(2);
    L(i+1)=y(3);
end
t=t0:h:tn; %Vector de subintervalos de tiempo
%Dibujo de gráficas
figure(1)
hold on
plot(t,G,'g') 
plot(t,Z,'b')
plot(t,L,'r') 
title('Población de gacelas(verde), cebras(azul) y leones(rojo) en función del tiempo');
xlabel('t(años)');
ylabel('población (millones de individuos)')
hold off
figure(2)
subplot(3,1,1)
plot(G,L,'b')
title('Población de leones en función de la de gacelas');
xlabel('G(milllones de gacelas)');
ylabel('L(millones de leones)');
subplot(3,1,2)
plot(Z,L,'y')
title('Población de leones en función de la de cebras');
xlabel('Z(milllones de cebras)');
ylabel('L(millones de leones)');
subplot(3,1,3)
plot(Z,G,'g')
title('Población de gacelas en función de la de cebras');
xlabel('Z(milllones de cebras)');
ylabel('G(millones de gacelas)');

Utilizando el método de Euler modificado con [math]h=0.1[/math], en un intervalo de tiempo tЄ[0,300] en el programa de matlab únicamente hay que sustituir el valor final del tiempo

2.2 Gráficas

2.2.1 Evolución e iteración de especies en tє[0,100]

En el estudio de la evolución de especies en el tiempo hemos obtenido las siguientes gráficas

En este gráfico se observan las variaciones poblacionales de las tres especies: leones como depredador y gacelas y cebras como presas a lo largo del tiempo. Podemos interpretar que la población de leones empieza a aumentar cuando la población de presas es máxima. Este crecimiento es sostenido hasta que el nivel de presas llega a un punto crítico, donde deja de haber suficientes presas para sostener la población de leones, lo que genera una disminución en el número de dichos depredadores, que a su vez da lugar a un aumento en el número de presas.

Por otra parte, se puede ver que las velocidades de crecimiento de las presas es mucho mayor que la de leones. Si esto no fuese así, conduciría a la extinción de las tres especies debido a un aumento desmesurado de la población de leones cuando todavía no hay presas suficientes.

Además de esto, otro aspecto a tener en cuenta es que a pesar de la existencia de dos especies distintas de presas, los predadores tienen predilección por una de ellas, en particular las gacelas, lo que provoca una disminución progresiva en los niveles máximos de éstas y por consiguiente un aumento de cebras en el tiempo. Esta situación se prolonga en el tiempo, lo que implica que transcurrido tiempo suficiente, las gacelas se extinguirán, quedando las cebras como única presa.

La iteración de poblaciones tє[0,100] será

Las dos primeras gráficas, en forma de espiral, muestran la evolución de leones en función de las gacelas y las cebras respectivamente. Los puntos de equilibrio en los que el sistema sería estable serían los orígenes (0,0) que los despreciamos, y los puntos donde se situarían los centros de las espirales. Como vemos la curva no converge al mismo por lo que no podríamos asegurar que ni el ecosistema formado por leones y gacelas ni el formado por leones y cebras sean estables.

La tercera gráfica indica la evolución de cebras en función de las gacelas, los puntos donde la gráfica crece son los puntos de crecimiento de la población total de presas. Si el intervalo de tiempo fuera lo suficientemente amplio, se podría comprobar que según aumenta el número de cebras, disminuye el número de gacelas, esto se ve muy claramente en los puntos de máxima y mínima población de presas.

2.2.2 Evolución e iteración de especies en tє[0,300]

En nuestro estudio de la evolución de las especies el resultado es el siguiente

Estas gráficas son iguales a las comentadas con anterioridad con la salvedad de que el intervalo de tiempo más amplio, motivo por el cual se observan mejor los crecimientos, decrecimientos y relaciones de las poblaciones, tanto en el tiempo como entre si.Los leones mantiene su evolución oscilatoria, las gacelas al alcanzar los 300 años casi se extinguen y las cebras alcanza su población máxima al final del periodo. Esto es así debido a que los leones se centran más en un tipo de presa haciendo que la otra aumente considerablemente su población.

La iteración de las poblaciones por tanto:

Las gráficas Depredador-Presa se muestran más alargadas que antes debido a que las poblaciones máximas adquieren mayor variación por ser mayor el tiempo estudiado.

3 Método de Euler

3.1 Programación método numérico

Utilizando el método de Euler con [math]h=0.1[/math], en un intervalo de tiempo tЄ[0,300], con una población inicial de [math]p_0=0.8[/math] millones de presas de un tipo, [math]q_0=2.4[/math] millones de presas de otro tipo y [math]d_0=0.2[/math] millon de depredadores y teniendo en cuenta que los valores que adoptan los parámetros son

[math] \left\{\begin{matrix}\ A_1=0.4, A_2=0.7\\ B_1=0.2, B_2=0.4\\ C_1=0.37, C_2=0.04, C_3=0.03\end{matrix}\right. [/math]

El codigo matlab que introduciremos será el siguiente:

clear all
clc
hold off
%Definiciones
t0=0; tn=300; % Extremos del intervalo temporal a estudiar
y0=[0.8;2.4;0.2]; %Valores iniciales
h=1; N=(tn-t0)/h; %Número y amplitud de subintervalos
y=y0; %Asignacion de valores iniciales
G(1)=y(1); %Iniciamos el vector de población G(gacelas)
Z(1)=y(2);%Iniciamos el vector de población Z(cebras)
L(1)=y(3);%Iniciamos el vector de población L(leones)
A=[0.35 0 0;0 0.3 0;0 0 -0.37]; %coeficientes de la parte lineal del sistema
%Algoritmo
for i=1:N
    y=y+h*(A*y+y(1)*y(3)*[-0.6;0;0.04]+y(2)*y(3)*[0;-0.5;0.035]);  %Aplicación del método de Euler 
    G(i+1)=y(1);
    Z(i+1)=y(2);
    L(i+1)=y(3);
end
t=t0:h:tn; %Vector de subintervalos de tiempo
%Dibujo de gráficas
figure(1)
hold on
plot(t,G,'g') 
plot(t,Z,'b')
plot(t,L,'r') 
title('Población de gacelas(verde), cebras(azul) y leones(rojo) en función del tiempo');
xlabel('t(años)');
ylabel('población (millones de individuos)');
hold off
figure(2)
subplot(2,1,1)
plot(G,L,'b')
title('Población de leones en función de la de gacelas');
xlabel('G(milllones de gacelas)');
ylabel('L(millones de leones)');
subplot(2,1,2)
plot(Z,L,'y')
title('Población de leones en función de la de cebras');
xlabel('Z(milllones de cebras)');
ylabel('L(millones de leones)');

3.2 Gráficas

3.2.1 Evolución de especies en el tiempo

Para una longitud de paso de h=0.1 tendremos la siguiente gráfica de evolución en el tiempo de las especies :

Por el contrario, para h=1 será:

3.2.2 Iteración de especies y trayectoria descrita

Para h=1 obtenemos el siguiente resultado

Para h=0,1 tendremos

3.3 Comparación de resultados

El método de Euler, a diferentes longitudes de paso h de subintervalo, realiza una aproximación por la recta tangente a determinadas funciones. En soluciones oscilatorias o cíclicas, como es nuestro caso, no obtenemos buenos resultados para longitudes de paso grandes como h = 1. Debido a las diferencias entre máximos y mínimos de la función, Euler no proporciona una aproximación adecuada para h excesivos, ya que el valor de dicha aproximación se aleja demasiado de la solución real acumulándose los errores cuanto mayor sea nuestra h.También es de destacar que para h = 1, el método no es capaz ni siquiera de aproximar la solución mínimamente, obteniéndose como solución rectas cuyos valores máximos son del orden de 10¹⁷⁰ millones de individuos.

Por el contrario tomamos una longitud de paso más pequeña, como h = 0,1, el método de Euler alcanza un nivel mayor de convergencia cuanto menor es el valor de h, pero esto puede suponer un problema, debido a que si tomamos un h demasiado pequeña puede que nuestro programa no responda, necesitando encontrar un valor de h adecuado para que aporte una buena precisión y pueda ser soportado por el programa.

De este método de integración, se puede decir que, si bien con un valor pequeño de h, la solución se puede aproximar a la real, no se obtiene ni de lejos la precisión que proporcionan otros métodos de integración para valores de h diez o cien veces mayores.

4 Método del trapecio

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. Al no ser posible su resolución, se descarta dicho método para comparar resultados.

5 Runge Kutta

5.1 Programación método numérico

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 [math]h=0.1[/math],con una población inicial de [math]p_0=3.5[/math] millones de presas de un tipo, [math]d_0=0.2[/math] millones de presas de otro tipo y [math]d_0=0.4[/math] millones de depredadores y teniendo en cuenta que los valores que adoptan los parámetros son los siguientes:

[math] 
\left\{\begin{matrix}\ A_1=0.4, A_2=0.7\\ B_1=0.2, B_2=0.4\\ C_1=0.37, C_2=0.04, C_3=0.035\end{matrix}\right.

[/math]

clear all
clc
hold off
%Valores iniciales

    t0=0; tN=500; %Extremos del intervalo a estudiar

    y0=[3.5;0.2;0.4]; %Valores de población iniciales

    h=0.1; N=(tN-t0)/h;%Nº y  amplitud de subintervalos

    y=y0;

    G(1)=y(1);%Valor inicial de G (población de gacelas)
    Z(1)=y(1);%Valor inicial de Z (población de cebras)
    L(1)=y(2);% Valor inicial de L (población de leones)

   A=[0.4 0 0;0 0.2 0;0 0 -0.37]; %Parámetros

    %Algoritmo

    for n=1:N

        k1=(A*y+y(1)*y(3)*[-0.7;0;0.04]+y(2)*y(3)*[0;-0.4;0.03]);

        k2=(A*(y+(h/2)*k1)+(y(1)+(h/2)*k1(1))*(y(3)+(h/2)*k1(3))*[-0.7;0;0.04]+(y(2)+(h/2)*k1(2))*(y(3)+(h/2)*k1(3))*[0;-0.4;0.03]);

        k3=(A*(y+(h/2)*k2)+(y(1)+(h/2)*k2(1))*(y(3)+(h/2)*k2(3))*[-0.7;0;0.04]+(y(2)+(h/2)*k2(2))*(y(3)+(h/2)*k2(3))*[0;-0.4;0.03]);

        k4=(A*(y+h*k3)+(y(1)+h*k3(1))*(y(3)+h*k3(3))*[-0.7;0;0.04]+(y(2)+h*k3(2))*(y(3)+h*k3(3))*[0;-0.4;0.03]);

        y=y+h/6*(k1+2*k2+2*k3+k4);

        G(n+1)=y(1);
        Z(n+1)=y(2);
        L(n+1)=y(3);

    end

    %Dibujo

    t=t0:h:tN;

    figure(1)

    hold on
   plot(t,G,'g') 
   plot(t,Z,'b')
   plot(t,L,'r') 
   title('Población de gacelas(verde), cebras(azul) y leones(rojo) en función del tiempo');
   xlabel('t(años)');
   ylabel('población (millones de individuos)')
   hold off

A diferencia de los otros métodos numéricos realizados con anterioridad, con el método de Runge-Kutta se obtienen soluciones más precisas y permite buenas aproximaciones con valores de h elevados, debido a que este método no depende de tangentes que se desvían de la solución y por tanto no habrá acumulaciones de errores que puedan hacer nuestra solución inservible.

5.2 Gráficas

Obtenemos como resultado la siguiente gráfica con la evolución en el tiempo del ecosistema:

En este gráfico podemos observar que debido a las interacciones entre las distintas especies las cebras se extinguen rápidamente. A raíz de esto, se obtiene un valor cuasi-constante a lo largo del tiempo de las poblaciones máximas y mínimas tanto de gacelas, como de leones.

Como podemos comprobar en el siguiente gráfico las tres especies : gacelas, cebras y leones van cambiando conforme avanza el tiempo. Las gacelas, representadas en color verde, las cebras representadas en color azul y los depredadores representados en color rojo. A medida que va aumentando la población de presas habrá más comida para los depredadores ,así que cuando la pobláción de presas es máxima, la población de depredadores aumentará. En nuestro caso, la población inicial de cebras es menos que la de leones por lo que acaba extinguiéndose. Por el contrario, en el caso de las gacelas, tienen una población mayor que la de los leones, por ello aunque éstos se alimentes de estas presas, la velocidad de reproducción y crecimiento de las gacelas es mayor y por ello nunca se extinguirán, dando lugar a ciclos de crecimiento y crecimiento siendo en sus máximos y mínimos constantes. La población de depredadores sigue aumentando hasta que el nivel de presas llega a un mínimo, entonces no habrá suficientes presas para alimentar a los depredadores, por lo que el numero de predadores disminuirá, y la cantidad de presas comenzará a subir.

6 Conclusiones

Con este ejemplo, y los resultados obtenidos por los distintos métodos de integración numérica, se puede deducir que mientras que euler no consigue aproximarse a la solución y no nos da estabilidad, runge-kutta y euler modificado si que lo hacen, dando un resultado casi exacto. Esto se debe al número de intervalos h que el método adopta (implícitamente) dentro de cada intervalo h de tiempo definido por nosotros: Euler únicamente utiliza el h que nosotros hemos definido, mientras que el método euler modificado utiliza 2h. Aquí es donde radica el porqué de la gráfica recta obtenida por el método de euler con h=1

Runge-kutta merece mención aparte ya que, aunque el método utilizado para integrar nuestras ecuaciones ha sido el de orden 4 (4h), dependiendo de la precisión que queramos obtener, podemos crear programas con el orden que queramos, lo que significa que podemos obtener la precisión que queramos. Aún así los órdenes más utilizados están entre el 2 y el 4, puesto que a mayor orden, más complicados serán los cálculos del ordenador, y más tiempo le llevará arrojar los resultados, y aunque éstos den resultados más exactos, la diferencia entre las precisiones de un "orden 4" y otro mayor son muy pequeñas.