Diferencia entre revisiones de «Logística con umbral (Grupo 10-C)»

De MateWiki
Saltar a: navegación, buscar
Línea 23: Línea 23:
 
clear all
 
clear all
 
clf
 
clf
 +
 
%DATOS DEL PROBLEMA
 
%DATOS DEL PROBLEMA
 
t0=0;
 
t0=0;
Línea 31: Línea 32:
 
M2=100;
 
M2=100;
 
r=0.04;
 
r=0.04;
 +
 
%Calculamos número de subintervalos
 
%Calculamos número de subintervalos
 
N=round((tN-t0)/h);
 
N=round((tN-t0)/h);
 
%Definimos la variable independiente
 
%Definimos la variable independiente
 
t=t0:h:tN;
 
t=t0:h:tN;
 +
 
%Ahora vamos a guardar los valores de la solución aproximada en el
 
%Ahora vamos a guardar los valores de la solución aproximada en el
 
%vector y
 
%vector y
y=zeros(1,N+1); %euler
+
y=zeros(1,N+1); %Euler
 
w=zeros(1,N+1); %Heun
 
w=zeros(1,N+1); %Heun
 
z=zeros(1,N+1); %Runge-kutta
 
z=zeros(1,N+1); %Runge-kutta
Línea 44: Línea 47:
 
z(1)=y0;
 
z(1)=y0;
 
for i=1:N
 
for i=1:N
  y(i+1)=y(i)+h*(-r*y(i)*(1-(y(i)/M1))*(1-(y(i)/M2))); %euler
+
  y(i+1)=y(i)+h*(-r*y(i)*(1-(y(i)/M1))*(1-(y(i)/M2))); %Euler
 
  k1=(-r*w(i)*(1-(w(i)/M1))*(1-(w(i)/M2))); %Heun
 
  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
 
  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
 
  w(i+1)=w(i)+(h/2)*(k1+k2); %Heun
  % Runge-kutta
+
  % Runge-Kutta
 
  A1=(-r*z(i)*(1-(z(i)/M1))*(1-(z(i)/M2)));
 
  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))
 
  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))
Línea 55: Línea 58:
 
  z(i+1)= z(i)+h/6*(A1+2*A2+2*A3+A4);
 
  z(i+1)= z(i)+h/6*(A1+2*A2+2*A3+A4);
 
end
 
end
 +
 
%Obtenemos la tabla de resultados
 
%Obtenemos la tabla de resultados
 
[t',y',w',z']
 
[t',y',w',z']
 
hold on
 
hold on
 +
 
%Gráfico
 
%Gráfico
 
plot(t,y,'r')
 
plot(t,y,'r')
Línea 67: Línea 72:
 
%Máximo error entre métodos  
 
%Máximo error entre métodos  
 
e1=max(abs(y-w))  %Euler y Heun
 
e1=max(abs(y-w))  %Euler y Heun
e2=max(abs(z-y))   %Runge-Kutta y Euler
+
e2=max(abs(z-y))   %Runge-Kutta y Euler
 
e3=max(abs(w-z))  %Heun y Runge-Kutta
 
e3=max(abs(w-z))  %Heun y Runge-Kutta
 
}}
 
}}

Revisión del 11:51 5 mar 2015

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

2 Problema de valor inicial

Nos encontramos con el siguiente problema de valor inicial. Empezamos por resolverlo con los métodos de Euler, Heun y Runge-Kutta de orden 4. Los datos proporcionados son los siguientes: el intervalo I=[0,100], el salto de paso es 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 ejercicio 8
%Euler
clear all
clf

%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;

%Ahora vamos a guardar los valores de la solución aproximada en el
%vector y
y=zeros(1,N+1); %Euler
w=zeros(1,N+1); %Heun
z=zeros(1,N+1); %Runge-kutta
y(1)=y0;
w(1)=y0;
z(1)=y0;
for i=1:N
 y(i+1)=y(i)+h*(-r*y(i)*(1-(y(i)/M1))*(1-(y(i)/M2))); %Euler
 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
 % 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'); % lo último es para que la leyenda
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


3 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 e 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.

Grafica para h=1
Grafica para h=1 zoom


Vemos que