Logística con umbral (Grupo 10-C)

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

Trabajo realizado por estudiantes
Título	Logística con umbral (Grupo 10-C)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Íñigo Díez García (1051) Pablo Molinero Brito (1172) Diego Navarro Gozalo (1049) Dariusz Adam Pabian (1187) Javier Santander Gimeno (1223) Pablo Vázquez Melgarejo (1090)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

\left \{ \begin{array}{c} y'= -ry(1-\frac{y}{M1})(1-\frac{y}{M2}), t>0 \\ y=y0 \end{array} \right \}

2 Problema de valor inicial

Nos encontramos con el siguiente problema de valor inicial. Lo resolveremos con los métodos de Euler, Heun y Runge-Kutta de orden 4. Los datos proporcionados son los siguientes: el intervalo I=[0,100], los saltos de paso son h=1, h=0.1, h=0.01 (con ello vamos a ver cómo cambia la precisión de la solución en función de los pasos de salto escogidos), y las constantes r, M1, M2, cuyos valores son 0.04, 30 y 100, respectivamente. El valor inicial de la función y es igual a 60.

%Modelo logístico 

%DATOS DEL PROBLEMA
t0=0;
tN=100;
y0=60;
h=input('introduce numero de paso: ');
M1=30;
M2=100;
r=0.04;

%Calculamos número de subintervalos
N=round((tN-t0)/h);

%Definimos la variable independiente
t=t0:h:tN;

%Creamos los vectores solución y,w,z 
y=zeros(1,N+1); %Euler
w=zeros(1,N+1); %Heun
z=zeros(1,N+1); %Runge-kutta

%Damos el valor inicial a la primera posicion del vector
y(1)=y0;
w(1)=y0;
z(1)=y0;

%Bucle
for i=1:N
 %Euler
 y(i+1)=y(i)+h*(-r*y(i)*(1-(y(i)/M1))*(1-(y(i)/M2))); 
 %Heun
 k1=(-r*w(i)*(1-(w(i)/M1))*(1-(w(i)/M2))); 
 k2=(-r*(w(i)+h*k1)*(1-((w(i)+h*k1)/M1))*(1-((w(i)+h*k1)/M2)));
 w(i+1)=w(i)+(h/2)*(k1+k2); 
 % Runge-Kutta
 A1=(-r*z(i)*(1-(z(i)/M1))*(1-(z(i)/M2)));
 A2=(-r*(z(i)+1/2*A1*h)*(1-((z(i)+1/2*A1*h)/M1))*(1-(z(i)+1/2*A1*h)/M2)); %z(i)+1/2*A1*h)*(1-(z(i)+1/2*A1*h))
 A3=(-r*(z(i)+1/2*A2*h)*(1-((z(i)+1/2*A2*h)/M1))*(1-(z(i)+1/2*A2*h)/M2));
 A4=(-r*(z(i)+A3*h)*(1-((z(i)+A3*h)/M1))*(1-(z(i)+A3*h)/M2));
 z(i+1)= z(i)+h/6*(A1+2*A2+2*A3+A4);
end

%Obtenemos la tabla de resultados
[t',y',w',z']
hold on

%Gráfico
plot(t,y,'r')
plot(t,w,'g')
plot(t,z)
legend('Euler','Heun','Runge-Kutta','Location','best'); 
hold off

%Máximo error entre métodos 
e1=max(abs(y-w))   %Euler y Heun
e2=max(abs(z-y))   %Runge-Kutta y Euler
e3=max(abs(w-z))   %Heun y Runge-Kutta

2.1 Interpretación numérica

Se observa que los tres métodos empleados en la resolución del sistema se aproximan entre sí, cuanto más se disminuye el número de salto de paso. Todos los métodos se consideran correctos, porque proporcionan información precisa de la solución. Para mejor visualización y compresión del ejercicio, presentamos tres gráficas. En cada una de ellas, aparecen las tres soluciones (cada una de ellas corresponde a un método). La primera gráfica corresponde a un salto de paso h=1; la segunda a un h=0.1, y la tercera a un h=0.01.

2.1.1 Observacion h=1

Grafica para h=1

Grafica para h=1 zoom

Vemos que los tres metodos se asemejan mucho a simple vista como se muestra en la grafica, para ver la precision real hemos calculado los errores entre los metodos, siendo e1 el error cometido entre Euler y Heun , e2 el error entre Runge-Kutta y Euler y e3 el error entre Heun y Runge-Kutta.

                                        e1 = 0.1860 ;  e2 = 0.1789  ;  e3 = 0.0074

2.1.2 Observacion h=0.1

Grafica para h=0.1

Grafica para h=0.1 zoom

Grafica para h=0.1 zoom entre Heun y Runge-Kutta

Vemos que los errores esta vez quedan mas pequeños ya que ha disminuido el intervalo, hemos dibujado dos graficas con zoom para apreciar mejor las diferencias entre los metodos.

                             e1 = 0.0177 ;  e2 = 0.0177 ;  e3 = 7.1680e-05

2.1.3 Observacion h=0.01

Grafica para h=0.01

Vemos como con este tamaño de paso tan pequeño los errores son muy bajos .No hemos hecho el zoom como en el apartado anterior ya que no se apreciarian mas diferencias.

                                       e1 = 0.0018 ;  e2 = 0.0018 ;  e3 = 7.1475e-07

2.2 Interpretacion desde el punto de vista poblacional

2.2.1 Apartado 3

Grafica para h=1

Vemos que con poblacion inicial moderada 60 el creicimiento de la especie es posible, y este se produce de manera logaritmica hasta estabilizarse la poblacion al alcanzar un valor cercano a 100.

2.2.2 Apartado 4

Grafica para heun con poblacion inicial 120

En este caso la poblacion inicial es mucho mayor y esto hace que sea imposible su superviviencia dado que no es posible el acceso a recursos para toda la poblacion. El medio esta superpoblado y por ello la poblacion se ve reducida de manera drastica hasta estabilizarse al aproximarse a valores cercanos a 100, como en el caso anterior en que la poblacion inicial es menor y esta crece.Podemos concluir que la capacidad de carga de la especie en el medio es un valor cercano a 100.

%Modelo logístico ejercicio 8. apartado 4
clear all
clf
%DATOS DEL PROBLEMA
t0=0;
tN=100;
y0=120;
h=0.1;
M1=30;
M2=100;
r=0.04;
%Calculamos número de subintervalos
N=round((tN-t0)/h);
% Definimos la variable independiente
t=t0:h:tN;

w=zeros(1,N+1); %Heun
w(1)=y0;

for i=1:N
 k1=(-r*w(i)*(1-(w(i)/M1))*(1-(w(i)/M2))); %Heun
 k2=(-r*(w(i)+h*k1)*(1-((w(i)+h*k1)/M1))*(1-((w(i)+h*k1)/M2))); %Heun t(i)+h+w(i)+k1*h
 w(i+1)=w(i)+(h/2)*(k1+k2); %Heun
end

hold on
%gráfico
plot(t,w)
legend('Heun','Location','best'); 
hold off

2.2.3 Apartado 5

Grafica para poblacion inicial 20

Con la poblacion inicial en 20 vemos como la especie no es capaz de sobrevivir en el medio cayendo a valores cercanos a 0.

%Modelo logístico ejercicio 8. apartado 5
clear all
clf
%DATOS DEL PROBLEMA
t0=0;
tN=100;
y0=20;
h=0.1;
M1=30;
M2=100;
r=0.04;
%Calculamos número de subintervalos
N=round((tN-t0)/h);
% Definimos la variable independiente
t=t0:h:tN;

w=zeros(1,N+1); %Heun
w(1)=y0;

for i=1:N
 k1=(-r*w(i)*(1-(w(i)/M1))*(1-(w(i)/M2))); %Heun
 k2=(-r*(w(i)+h*k1)*(1-((w(i)+h*k1)/M1))*(1-((w(i)+h*k1)/M2))); 
 w(i+1)=w(i)+(h/2)*(k1+k2);
end

hold on
%gráfico
plot(t,w)
legend('Heun','Location','best'); 
hold off

3 Sistema poblacional no lineal

3.1 Introduccion

3.2 Resolucion numerica

3.3 Casos Particulares

3.3.1 Caso 1

Para paso h=0.1 el metodo de Euler es inestable y no nos sirve para resolver el problema. Heun si es valido la grafica muestra nos permite observar el desarrollo de dos poblacion una de ellas un depredador y otra una presa. El aumento y descenco de ellas se ve relacionado en tanto que al aumentar mucho los depredadores las presas reducen mucho su numnero, esto causa que a largo plazo las presas tenga que disminuir su numero para sobrevivir. De manera inversa ocurre con las presas cuyo numero aumenta cuando la cantidad de presas se reduce. Las graficas nor permiten ver que aunque la cantidad de presas o depredadores se vea muy reducida nunca llega al 0 lo cual permite que se recuperen.

Al utilizar un tamaño de paso igual a 0.0001 el metodo de Euler se vuelve preciso y por tanto valido para resolver el problema, al haber reducido tanto el tamaño de paso el trabajo que tiene que realizar el ordenador es muy grande por lo que el metodo aunque sea valido no es practico.