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

Revisión actual del 09:43 7 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 interaciones 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 en Matlab

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 y relación entre 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, éstos empiezan a pasar hambre 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. Este proceso desencadena un nuevo incremento de la población de depredadores y, conforme va pasando el tiempo, se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las especies.

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 relación de especies en 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 y relación entre 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 relación entre las especies 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 en Matlab

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.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 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 y trayectoria descrita h = 0.1

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

En esta primera gráfica podemos observar en primer lugar abunda la especie presa, gacelas y cebras, siendo las cebras (representadas en color azul) las que más aumentan en población debido al mayor número de ejemplares. En segundo lugar la especie depredadora, leones (representadas en color rojo), tiene alimento, se fortalece y aumenta la población de los depredadores, alcanzando su primer máximo como consecuencia de la disminución de las presas. Cada máximo de cada función oscilatoria están desfasadas unas con respecto a otras debido a la distinta evolución de cada especie.

Como en el caso de Euler modificado, debido a que las curvas no convergen al centro de las espirales no podemos asegurar que estos ecosistemas sean estables.

Los resultados obtenidos de calcular los instantes en los que se dan los máximos y mínimos poblacionales de cada especie, así como los valores de pico de dichas poblaciones vienen recogidos en la siguiente tabla:

3.2.2 Evolución de especies en el tiempo y trayectoria descrita h = 1

Se puede apreciar que dado que el paso "h" es excesivamente grande, las 3 poblaciones, especialmente las de las presas, 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.

Los resultados obtenidos de calcular los instantes en los que se dan los máximos y mínimos poblacionales de cada especie, así como los valores de pico de dichas poblaciones vienen recogidos en la siguiente tabla:

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 de las funciones, con la salvedad de que si otros métodos de integracion subdividen cada paso, minimizando errores, euler no lo hace por lo que cada error que comete se va acumulando, ocasionando la gran imprecisión que este método tiene para resolver sistemas. 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. Otra incoherencia es el hecho de que el modelo presente una población negativa en los leones, suceso que en la realidad no es posible.

Por el contrario, si tomamos una longitud de paso más pequeña como h = 0,1, la solución por el método de Euler alcanza un nivel mayor de convergencia cuanto menor es el valor de h (aunque no vamos a conseguir con este método la precisión obtenida con otros), pero esto puede suponer un problema, debido a que si tomamos un h demasiado pequeña puede que nuestro programa no responda por el elevado número de operaciones a realizar, lo que provoca que este método no sea útil para resolver sistemas de ecuaciones no lineales.

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 en Matlab

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.03\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(2);%Valor inicial de Z (población de cebras)
    L(1)=y(3);% 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:

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.A raíz de esto, se obtiene un valor cuasi-constante a lo largo del tiempo de los valores máximos y mínimos de las poblaciones, tanto de gacelas como de leones. 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.

Los resultados obtenidos de calcular los instantes en los que se dan los máximos y mínimos poblacionales de cada especie, así como los valores de pico de dichas poblaciones vienen recogidos en la siguiente tabla:

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 el método de 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 divisiones que el método adopta (implícitamente) dentro de la longitud de paso h definida por nosotros (orden 1=h original; orden 2=h/2; orden 3=h/3, etc.):

Euler únicamente utiliza el h que nosotros hemos definido, mientras que el método Euler modificado utiliza una longitud h/2. 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 (longitud del paso de curva adoptado por el método=h/4), 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.