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

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]