Logística con Umbral (15-A)

Trabajo realizado por estudiantes
Título	Logística con Umbral (Grupo 15-A)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	• Javier Colorado Martínez (1278) • Álvaro Llera Fernández (501) • Francisco Javier Alcaraz de Amuriza (567) • Daniel Ballesteros Gálvez (1376) • Alberto Rodríguez Fernández (357) • Francisco Javier Barral González (632)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Nuestro Problema de Valor Inicial (P.V.I.) viene definido por la siguiente ecuación: [math]\left\{\begin{matrix}y’=-ry(1-{y\over M_1})(1-{y\over M_2}), t\gt0\\\ltbr/\gt\\y=y_0\end{matrix}\right.[/math] Donde \(M_1=30\), \(M_2=100\), \(r=0.04\) e \(y_0=60\). La ecuación a estudiar está definida en el intervalo \(I=[0,100]\).

2 Resolución numérica

A continuación, se va a proceder a modelizar dicha ecuación mediante los métodos de Euler, Heun y Runge-Kutta de orden 4. El procedimiento por el que se ha optado ha consistido en introducir en MATLAB cada uno de los métodos con los diferentes pasos, a saber, \(h=1\), \(h=0.1\) y \(h=0.01\).

2.1 Euler

%MÉTODO DE EULER
hold on
%Datos conocidos e intervalo de estudio.
t0=0; tN=100; r=0.04; M1=30; M2=100; y0=60;
%-------------------------------------------------------
%CASO 1: paso de 1
h1=1;
N1=(tN-t0)/h1;
t1=t0:h1:tN;
y1=zeros(1,N1+1);
y1(1)=y0;
for n=1:N1
  y1(n+1)=y1(n)+h1*(-r*y1(n)*(1-y1(n)/M1)*(1-y1(n)/M2));
end
%-------------------------------------------------------
%CASO 2: paso de 0.1
h2=0.1;
N2=(tN-t0)/h2;
t2=t0:h2:tN;
y2=zeros(1,N2+1);
y2(1)=y0;
for n=1:N2
  y2(n+1)=y2(n)+h2*(-r*y2(n)*(1-y2(n)/M1)*(1-y2(n)/M2));
end
%------------------------------------------------------
%CASO 3: PASO DE 0.01
h3=0.01;
N3=(tN-t0)/h3;
t3=t0:h3:tN;
y3=zeros(1,N3+1);
y3(1)=y0;
for n=1:N3
  y3(n+1)=y3(n)+h3*-r*y3(n)*(1-y3(n)/M1)*(1-y3(n)/M2);
end
%-----------------------------------------------------
%Representación en gráficas.
plot(t1,y1,'r') %Caso 1.
plot(t2,y2,'g') %Caso 2.
plot(t3,y3,'b') %Caso 3.
title('Método de Euler con todos los casos')
xlabel('Tiempo')
ylabel('Evolución de la y')
legend('Caso 1 Euler','Caso 2 Euler','Caso 3 Euler','Location','best')
hold off

Como podemos apreciar en las gráficas

2.2 Heun

%MÉTODO DE HEUN
hold on
%Datos conocidos e intervalo de estudio.
t0=0; tN=100; r=0.04; M1=30; M2=100; y0=60;
%--------------------------------------------------------------
%CASO 1: paso de 1
h1=1;
N1=(tN-t0)/h1;
t1=t0:h1:tN;
y1=zeros(1,N1+1);
y1(1)=y0;
for n=1:N1
  k1=-r*y1(n)*(1-y1(n)/M1)*(1-y1(n)/M2);
  k2=-r*(y1(n)+k1*h1)*(1-(y1(n)+k1*h1)/M1)*(1-(y1(n)+k1*h1)/M2);
  y1(n+1)=y1(n)+h1/2*(k1+k2);
end
%---------------------------------------------------------------
%CASO 2: PASO DE 0.1
h2=0.1;
N2=(tN-t0)/h2;
t2=t0:h2:tN;
y2=zeros(1,N2+1);
y2(1)=y0;
for n=1:N2
  k1=-r*y2(n)*(1-y2(n)/M1)*(1-y2(n)/M2);
  k2=-r*(y2(n)+k1*h2)*(1-(y2(n)+k1*h2)/M1)*(1-(y2(n)+k1*h2)/M2);
  y2(n+1)=y2(n)+h2/2*(k1+k2);
end
%---------------------------------------------------------------
%CASO 3: PASO DE 0.01
h3=0.01;
N3=(tN-t0)/h3;
t3=t0:h3:tN;
y3=zeros(1,N3+1);
y3(1)=y0;
for n=1:N3
  k1=-r*y3(n)*(1-y3(n)/M1)*(1-y3(n)/M2);
  k2=-r*(y3(n)+k1*h3)*(1-(y3(n)+k1*h3)/M1)*(1-(y3(n)+k1*h3)/M2);
  y3(n+1)=y3(n)+h3/2*(k1+k2);
end
%---------------------------------------------------------------
%Representación en gráficas.
plot(t1,y1,'r') %Caso 1.
plot(t2,y2,'g') %Caso 2.
plot(t3,y3,'b') %Caso 3.
title('Método de Heun con todos los casos')
xlabel('Tiempo')
ylabel('Evolución de la y')
legend('Caso 1 Heun','Caso 2 Heun','Caso 3 Heun','Location','best')
hold off

Como podemos apreciar

2.3 Runge-Kutta de orden 4

%MÉTODO DE RUNGE-KUTTA DE ORDEN 4
hold on
%Datos conocidos e intervalo de estudio.
t0=0; tN=100; r=0.04; M1=30; M2=100; y0=60;
%--------------------------------------------------------------------------
%CASO 1: paso de 1
h1=1;
N1=(tN-t0)/h1;
t1=t0:h1:tN;
y1=zeros(1,N1+1);
y1(1)=y0;
for n=1:N1
  k1=-r*y1(n)*(1-y1(n)/M1)*(1-y1(n)/M2);
  k2=-r*(y1(n)+1/2*k1*h1)*(1-(y1(n)+1/2*k1*h1)/M1)*(1-(y1(n)+1/2*k1*h1)/M2);
  k3=-r*(y1(n)+1/2*k2*h1)*(1-(y1(n)+1/2*k2*h1)/M1)*(1-(y1(n)+1/2*k2*h1)/M2);
  k4=-r*(y1(n)+k3*h1)*(1-(y1(n)+k3*h1)/M1)*(1-(y1(n)+k3*h1)/M2);
  y1(n+1)=y1(n)+h1/6*(k1+2*k2+2*k3+k4);
end
%--------------------------------------------------------------------------
%CASO 2: paso de 0.1
h2=0.1;
N2=(tN-t0)/h2;
t2=t0:h2:tN;
y2=zeros(1,N2+1);
y2(1)=y0;
for n=1:N2
  k1=-r*y2(n)*(1-y2(n)/M1)*(1-y2(n)/M2);
  k2=-r*(y2(n)+1/2*k1*h2)*(1-(y2(n)+1/2*k1*h2)/M1)*(1-(y2(n)+1/2*k1*h2)/M2);
  k3=-r*(y2(n)+1/2*k2*h2)*(1-(y2(n)+1/2*k2*h2)/M1)*(1-(y2(n)+1/2*k2*h2)/M2);
  k4=-r*(y2(n)+k3*h2)*(1-(y2(n)+k3*h2)/M1)*(1-(y2(n)+k3*h2)/M2);
  y2(n+1)=y2(n)+h2/6*(k1+2*k2+2*k3+k4);
end
%--------------------------------------------------------------------------
%CASO 3: paso de 0.01
h3=0.01;
N3=(tN-t0)/h3;
t3=t0:h3:tN;
y3=zeros(1,N3+1);
y3(1)=y0;
for n=1:N3
  k1=-r*y3(n)*(1-y3(n)/M1)*(1-y3(n)/M2);
  k2=-r*(y3(n)+1/2*k1*h3)*(1-(y3(n)+1/2*k1*h3)/M1)*(1-(y3(n)+1/2*k1*h3)/M2);
  k3=-r*(y3(n)+1/2*k2*h3)*(1-(y3(n)+1/2*k2*h3)/M1)*(1-(y3(n)+1/2*k2*h3)/M2);
  k4=-r*(y3(n)+k3*h3)*(1-(y3(n)+k3*h3)/M1)*(1-(y3(n)+k3*h3)/M2);
  y3(n+1)=y3(n)+h3/6*(k1+2*k2+2*k3+k4);
end
%--------------------------------------------------------------------------
%Representación en gráficas.
plot(t1,y1,'r') %Caso 1.
plot(t2,y2,'g') %Caso 2.
plot(t3,y3,'b') %Caso 3.
title('Método de Runge-Kutta con todos los casos')
xlabel('Tiempo')
ylabel('Evolución de la y')
legend('Caso 1 RK4','Caso 2 RK4','Caso 3 RK4','Location','best')
hold off

Como podemos en la gráfica

2.4 Interpretaciones

De todas las gráficas representadas anteriormente, se puede apreciar que aquellas que corresponden a un paso de \(h=1\) (el caso 1, considerado como un paso amplio), se alejan de la solución exacta, a diferencia de lo que sucede con pasos menores (casos 2 (\(h=0.1\)) y 3 (\(h=0.01\))). A su vez, la solución más cercana a la exacta es la obtenida mediante el método de Runge-Kutta de orden 4.

2.5 Dinámica de poblaciones

Desde un punto de vista de la dinámica de poblaciones, dado que se trata de un método logístico, se puede apreciar que la población no crece indefinidamente, sino que llega a un punto en el que el ecosistema no admite más individuos de la especie en cuestión.

3 Método de Heun para diferentes poblaciones iniciales

%MÉTODO DE HEUN CON DISTINTOS VALORES INICIALES
clf; clear;
% Datos Iniciales.
t0=0; tN=100; r=0.04; M1=30; M2=100; 
y0=120; z0=20; 
%Paso e intervalo de estudio.
h=0.1;
N=(tN-t0)/h;
t=t0:h:tN;
%--------------------------------------------------------------------------
%Heun con valor inicial 120.
y1=zeros(1,N+1);
y1(1)=y0;
for n=1:N
  k1=-r*y1(n)*(1-y1(n)/M1)*(1-y1(n)/M2);
  k2=-r*(y1(n)+k1*h)*(1-(y1(n)+k1*h)/M1)*(1-(y1(n)+k1*h)/M2);
  y1(n+1)=y1(n)+h/2*(k1+k2);
end
%--------------------------------------------------------------------------
%Heun con valor incial 20.
z1=zeros(1,N+1);
z1(1)=z0;
for n=1:N
  K1=-r*z1(n)*(1-z1(n)/M1)*(1-z1(n)/M2);
  K2=-r*(z1(n)+K1*h)*(1-(z1(n)+K1*h)/M1)*(1-(z1(n)+k1*h)/M2);
  z1(n+1)=z1(n)+h/2*(K1+K2);
end
%--------------------------------------------------------------------------
%Comparación de gráficas.
hold on
plot(t,y1,'r')
plot(t,z1,'g')
title('Comparativa del método de Heun')
xlabel('Tiempo')
ylabel('Evolución de la población')
legend('Heun 120','Heun 20','Location','best')
hold off

3.1 Interpretación

Según lo que se aprecia en la gráfica, las poblaciones siempre disminiuyen hasta establizarse en un valor determinado. En la primera de las poblaciones, la cual es mayor (en este caso, 120 individuos), ésta se estabiliza con una mayor rapidez de la que lo hace en el segundo caso, con una población de 20, donde tiende a extinguirse hasta llegar a 0.

4 Modelos de competencia

[math]\left\{\begin{matrix}x’=a_1x-b_1x^2-c_1xy, \\\ltbr/\gt\\ y’=a_2x-b_2x^2-c_2xy \end{matrix}\right.[/math]

%%SISTEMA NO LINEAL POR EL MÉTODO DE EULER
% x'=a1*x-b1*x^2-c1*x*y
% y'=a2*y-b2*y^2-c2*x*y
clear all; clf
%Introducción de los coeficientes, instante final y pasos de los diferentes casos.
a1=input('a1: '); a2=input('a2: ');
b1=input('b1: '); b2=input('b2: ');
c1=input('c1: '); c2=input('c2: ');
h=input('Paso: ');
tN=input('Limite final: ');
%Números de subintervalos.
t0=0; 
N=round((tN-t0)/h);
t=[t0:h:tN];
x=zeros(1,N+1);
y=zeros(1,N+1);
%Valores iniciales..
x0=2; y0=7;
x(1)=x0;
y(1)=y0;
%Bucle planteado con el método de Euler.
for i=1:N
  x(i+1)=x(i)+h*(a1*x(i)-b1*(x(i))^2-c1*x(i).*y(i));
  y(i+1)=y(i)+h*(a2*y(i)-b2*(y(i))^2-c2*x(i).*y(i));
end
%Para conocer los valores máximos y mínimos de las funciones.
[Mx, pMx]=max(x); [My, pMy]=max(y);
[mx, pmy]=min(x); [mx, pmy]=min(y);
%Representación de las gráficas de ambas poblaciones.
hold on
plot(t,x)
plot(t,y, 'r')
legend('Animal 1', 'Animal 2', 'Location' , 'best')
xlabel('Tiempo');
ylabel('Población');
title('Euler')
hold off

%SISTEMA NO LINEAL POR EL MÉTODO DE HEUN
% x'=a1*x-b1*x^2-c1*x*y
% y'=a2*y-b2*y^2-c2*x*y
clear all; clf
%Introducción de los coeficientes, instante final y pasos de los diferentes casos.
a1=input('a1: '); a2=input('a2: '); 
b1=input('b1: '); b2=input('b2: ');
c1=input('c1: '); c2=input('c2: ');
h=input('Paso: ');
tN=input('Limite final: ');
%Números de subintervalos.
t0=0;
N=round((tN-t0)/h);
t=[t0:h:tN];
x=zeros(1,N+1);
y=zeros(1,N+1);
%Valores iniciales..
x0=2; y0=7;
x(1)=x0;
y(1)=y0;
%Bucle planteado con el método de Heun.
for i=1:N
  k1x=a1*x(i)-b1*x(i)^2-c1*x(i)*y(i);
  k1y=a2*y(i)-b2*y(i)^2-c2*x(i)*y(i);
  k2x=a1*(x(i)+k1x*h)-b1*(x(i)+k1x*h)^2-c1*(x(i)+k1x*h)*y(i);
  k2y=a2*(y(i)+k1y*h)-b2*(y(i)+k1y*h)^2-c2*x(i)*(y(i)+k1y*h);
  
  x(i+1)=x(i)+h/2*(k1x+k2x);
  y(i+1)=y(i)+h/2*(k1y+k2y);
end
%Para conocer los valores máximos y mísimos de las funciones.
[Mx, pMx]=max(x); [My, pMy]=max(y);
[mx, pmx]=min(x); [mx, pmy]=min(y);
%Representación de las gráficas de ambas poblaciones.
hold on
plot(t,x)
plot(t, y,'r')
legend('Animal 1', 'Animal 2', 'Location', 'best')
xlabel('Tiempo');
ylabel('Población');
title('Heun')
hold off

4.1 Situaciones posibles

4.1.1 Parasitismo

Los coeficientes no se corresponden con la realidad porque el animal 1 renace pero de 0 a 15 sería un parasitismo. El problema de Euler es que es menos aproximado que cualquiera de los otros métodos. Hay un máximo de 50.773 animal 1 en 3.1. Hay un máximo de 28.612 animal 2 en 2. No sólo desaparece sino que renace y todo. El crecimiento del animal 1 se ve afectado positivamente por el crecimiento del animal 2 pero el crecimiento del animal 2 se ve afectado negativamente por el crecimiento del animal 1 Ecosistema no estable ya que el animal 2 desaparece y se puede apreciar que el animal 1 va por el mismo camino.

asa

4.1.2 Neutralismo

c1 y c2 son 0 por lo que no hay interacciones entre ellas. Neutralismo porque se elimina el término que relaciona las dos especies. No existen máximos y mínimos salvo en los extremos del intervalo de estudio (teorema de Weierstrass). No tiende a desaparecer ninguna especie. Tienden a ser estables las dos especies. Al tener parámetros c igual a 0 las especies no interactúan entre ellas quedando ecuaciones diferenciales de primer grado independientes. Así que el crecimiento de cada una de ellas no afecta sobre la otra. El ecosistema es estable ya que las dos especies tienden a un número.

asa

4.1.3 Simbiosis o cooperación

Simbiosis ya que ambas especies empiezan a crecer al final del dominio de tiempo dado. Hay un mínimo de 3.366 animal 2 en 9. No tiende a desaparecer ninguna especie. En nuestro intervalo de estudio aparece un crecimiento, si se establece un intervalo mayor, se ve que las gráficas tienden a un número. Se afectan positivamente porque los coeficientes son negativos y además en la gráfica se ve que las dos especies crecen al final del intervalo. Es un ecosistema estable porque tienden ambas a un número.

asa

4.1.4 Competición

Es una competición ya que hemos comparado con una gráfica suponiendo que no hay interacción y viendo que ambos crecían más que en la gráfica pedida. Pero gráficamente con solo la referencia de la gráfica pedida parece un amensalismo ya que el animal 1 perjudica el crecimiento del animal 2 sin verse el 1 afectado. El animal 2 tiene un máximo de 8.576 en 1.9. Tienden a estabilizarse las dos especies. Se afectan negativamente entre sí. Es un ecosistema estable ya que las dos especies tienden a un número.

asa

4.1.5 Amensalismo

Es un amensalismo ya que el animal 1 crece sin verse afectado por el animal 2 mientras que el animal 2 muere estrepitosa y vergonzosamente entre terribles dolores y sufrimientos. Teorema de Weierstrass. En animal 2 desaparece mientras que el animal 1 tiende a estabilizarse en 25. Al animal 1 no le afecta (c_1=0) mientras que al animal 2 le afecta negativamente extinguiéndose. Es un sistema estable ya que ambos animales tienden a un número (aunque uno de esos números sea 0.

asa

4.1.6 Comensalismo

Es un comensalismo debido a que el crecimiento del animal 1 no se ve afectado mientras que el animal 2 se ve beneficiado a costa del primero pero sin perjudicarlo ni beneficiarlo. Teorema de Weierstrass. Tienden a estabilizarse ambas especies. En el animal 1 no afecta ya que el parámetro es 0 quedando una ecuación diferencial de primer grado independiente, el animal 2 se ve beneficiado por el crecimiento del animal 1. Es un ecosistema estable ya que ambos animales tienden a un número.

asa

4.2 Comparación de especies

Hacemos uso del siguiente código de MATLAB:

%SISTEMA NO LINEAL POR EL MÉTODO DE HEUN
% x'=a1*x-b1*x^2-c1*x*y
% y'=a2*y-b2*y^2-c2*x*y
clear all; clf
%Introducción de los coeficientes, instante final y pasos de los diferentes casos.
a1=input('a1: '); a2=input('a2: '); 
b1=input('b1: '); b2=input('b2: ');
c1=input('c1: '); c2=input('c2: ');
h=input('Paso: ');
tN=input('Limite final: ');
%Números de subintervalos.
t0=0;
N=round((tN-t0)/h);
t=[t0:h:tN];
x=zeros(1,N+1);
y=zeros(1,N+1);
%Valores iniciales..
x0=2; y0=7;
x(1)=x0;
y(1)=y0;
%Bucle planteado con el método de Heun.
for i=1:N
  k1x=a1*x(i)-b1*x(i)^2-c1*x(i)*y(i);
  k1y=a2*y(i)-b2*y(i)^2-c2*x(i)*y(i);
  k2x=a1*(x(i)+k1x*h)-b1*(x(i)+k1x*h)^2-c1*(x(i)+k1x*h)*y(i);
  k2y=a2*(y(i)+k1y*h)-b2*(y(i)+k1y*h)^2-c2*x(i)*(y(i)+k1y*h);
  
  x(i+1)=x(i)+h/2*(k1x+k2x);
  y(i+1)=y(i)+h/2*(k1y+k2y);
end
%Para conocer los valores máximos y mísimos de las funciones.
[Mx, pMx]=max(x); [My, pMy]=max(y);
[mx, pmx]=min(x); [mx, pmy]=min(y);
%Representación de las gráficas de ambas poblaciones.

subplot(1,2,1)
hold on
plot(t,x)
plot(t, y,'r')
xlabel('Tiempo');
ylabel('Animales');
legend('Animal 1', 'Animal 2', 'Location', 'best')
title('Heun')
subplot(1,2,2)
plot(x,y,'k')
xlabel('Animal 1');
ylabel('Animal 2');
title('Comparación de especies')
hold off

4.2.1 Parasitismo

En cuanto el animal 2 aumenta su población, el animal 1 comienza a aumentar su población a su costa, se trata de un parasitismo. Cuando el animal 2 está casi extinto, el animal 1 empieza también a desaparecer.

4.2.2 Neutralismo

Aunque pudiera parecer un amensalismo ya que cuando una especie aumenta mientras que la otra desciende, se trata de un neutralismo ya que ninguna de las dos especies tiende a desaparecer.

4.2.3 Simbiosis o cooperación

Se trata de una simbiosis ya que en la parte final de la gráfica ambas especies tienden a crecer.

4.2.4 Competición

A primera vista se podría sospechar que se trata de un parasitismo. No obstante, si fuera así, al verse reducido el número de individuos de la especie 2, el animal 1 (supuesto parásito) debería tender a disminuir también su población al quedarse sin huésped. De ahí deducimos que no es un caso de parasitismo, porque al verse disminuida la especie 2, la especie 1 no tiende a disminuir su población, sino que ambas se estabilizan. Podría tratarse de un caso de competición, porque el crecimiento de ambas especies llega a un máximo, a partir del cual el animal 1 sigue en tendencia positiva de crecimiento, mientras que el animal 2 decrece, aunque no se ha llegado a un consenso claro en esta interpretación.

4.2.5 Amensalismo

En este caso se puede apreciar un amensalismo, ya que la especie 2 se ve reducida a medida que aumenta el número de individuos de la especie 1. Como observación, en la gráfica se ve que cuantos más animales 1 hay, más rápido decrece la población del animal 2, y además, cuantos menos individuos hay de la especie 2, más despacio decrece ésta.

4.2.6 Comensalismo

Vemos que una de las especies ve incrementada su población en proporción al número de individuos de la segunda especie, mientras ésta no se ve afectada por la primera.