Diferencia entre revisiones de «Logística con umbral A7»

clear all
%DATOS DEL PROBLEMA
t0=0;tN=100;    %tiempos inicial y final
y0=60;          %condición inicial
%El tamaño de paso lo haremos con h=1,0.1,0.01
h=input('Introduce tamaño de paso:');      %Tamaño de paso
N=(tN-t0)/h;            %discretización
r=0.04;m1=30;m2=100;  %constantes de la ecuación
t=t0:h:tN;          %vector de tiempos
%Vector y su condición inicial en cada metodo
y=zeros(1,N+1); y(1)=y0;   %Runge-Kutta
z=zeros(1,N+1);z(1)=y0;    %Euler
he=zeros(1,N+1);he(1)=y0;  %Heun
for i=1:N
 %Método Runge-Kutta:
 K1=-r*y(i)*(1-y(i)/m1)*(1-y(i)/m2);
 K2=-r*(y(i)+1/2*K1*h)*(1-((y(i)+1/2*K1*h)/m1))*(1-((y(i)+1/2*K1*h)/m2));
 K3=-r*(y(i)+1/2*K2*h)*(1-((y(i)+1/2*K2*h)/m1))*(1-((y(i)+1/2*K2*h)/m2));
 K4=-r*(y(i)+K3*h)*(1-((y(i)+K3*h)/m1))*(1-((y(i)+K3*h)/m2));
 y(i+1)= y(i)+h/6*(K1+2*K2+2*K3+K4);
 %Método Euler:
 z(i+1)=z(i)+h*(-r*(z(i)+1/2*K1*h)*(1-(z(i)+1/2*K1*h)/m1)*(1-(z(i)+1/2*K1*h)/m2));
 %Método de Heun:
 Kh1=-r*he(i)*(1-he(i)/m1)*(1-he(i)/m2);
 Kh2=-r*(he(i)+Kh1*h)*(1-((he(i)+Kh1*h)/m1))*(1-((he(i)+Kh1*h)/m2));
 he(i+1)= he(i)+h/2*(Kh1+Kh2);      
end
%Dibujamos las distintas gráficas:
subplot(3,1,1);
plot(t,y,'b','linewidth',2)
legend('RUNGE-KUTTA','location','best')
subplot(3,1,2);
plot(t,z,'r','linewidth',2)
legend('EULER','location','best')
subplot(3,1,3);
plot(t,he,'g','linewidth',2)
legend('HEUN','location','best')
%Cálculo errores
C1=max(abs(y-z))    %error entre Runge-Kutta y Euler
C2=max(abs(z-he))   %error entre Euler y Heun
C3=max(abs(he-y))   %error entre Heun y Runge-Kutta

clear all
%DATOS DEL PROBLEMA
t0=0;tN=100;    %tiempos inicial y final
y0=120; y02=20;       %condición inicial
%El tamaño de paso lo haremos con h=1,0.1,0.01
h=0.1;                %Tamaño de paso
N=(tN-t0)/h;          %discretización
r=0.04;m1=30;m2=100;  %constantes de la ecuación
t=t0:h:tN;            %vector de tiempos
%Vector y su condición inicial en cada y0
he=zeros(1,N+1);he(1)=y0;  
he2=zeros(1,N+1);he2(1)=y02;
for i=1:N
    %Con y0=120
 K1=-r*he(i)*(1-he(i)/m1)*(1-he(i)/m2);
 K2=-r*(he(i)+K1*h)*(1-((he(i)+K1*h)/m1))*(1-((he(i)+K1*h)/m2));
 he(i+1)= he(i)+h/2*(K1+K2);     
    %Con y0=20
 Kh1=-r*he2(i)*(1-he2(i)/m1)*(1-he2(i)/m2);
 Kh2=-r*(he2(i)+Kh1*h)*(1-((he2(i)+Kh1*h)/m1))*(1-((he2(i)+Kh1*h)/m2));
 he2(i+1)= he2(i)+h/2*(Kh1+Kh2);
end
%dibujamos las graficas:
hold on
plot(t,he,'g','linewidth',2)
plot(t,he2,'b','linewidth',2)
legend('y0=120','y0=20','location','best')
hold off
xlabel('tiempo')
ylabel('población')