Diferencia entre revisiones de «Modelo predador-presa (Grupo 16B)»
De MateWiki
(→Resolución del modelo mediante el método de Euler modificado) |
|||
| Línea 4: | Línea 4: | ||
==Interpretación== | ==Interpretación== | ||
==Resolución del modelo mediante el método de Euler modificado== | ==Resolución del modelo mediante el método de Euler modificado== | ||
| + | {{matlab|codigo= | ||
| + | % Datos del problema | ||
| + | t0=0; | ||
| + | tf=300; | ||
| + | %Coeficientes | ||
| + | A1=0.4; A2=0.3; A3=0.35; | ||
| + | B1=0.3; B2=0.05; | ||
| + | D=0.1; | ||
| + | C1=0.28; C2=0.045; | ||
| + | %Condiciones iniciales | ||
| + | p0=3.5; | ||
| + | d0=1; | ||
| + | e0=1.2; | ||
| + | Y0=[p0 d0 e0]'; | ||
| + | % Datos discretización | ||
| + | h=0.1; | ||
| + | N=(tf-t0)/h; | ||
| + | % Vectores de tiempo y soluciones app | ||
| + | t=t0:h:tf; | ||
| + | % Inicialización | ||
| + | Y(:,1)=Y0; | ||
| + | YY=Y0; | ||
| + | for n=1:N | ||
| + | % Calculamos k1=f(t_n,y_n) | ||
| + | k1=[A1*YY(1)-A2*YY(1)*YY(2)-A3*YY(1)*YY(3);-B1*YY(2)+B2*YY(1)*YY(2)-D*YY(2)*YY(3);-C1*YY(3)+C2*YY(1)*YY(3)-D*YY(2)*YY(3)]; | ||
| + | % Calculamos k2=f(t_n+h,y_n+h*k1) | ||
| + | k2=[A1*(YY(1)+h*k1(1))-A2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-A3*(YY(1)+h*k1(1))*(YY(3)+h*k1(3));-B1*(YY(2)+h*k1(2))+B2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3));-C1*(YY(3)+h*k1(3))+C2*(YY(1)+h*k1(1))*(YY(3)+h*k1(3))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3))]; | ||
| + | YY=YY+(h/2)*(k1+k2); | ||
| + | Y(:,n+1)=YY; | ||
| + | end | ||
| + | figure (1) | ||
| + | hold on | ||
| + | plot(t,Y(1,:),'r') | ||
| + | plot(t,Y(2,:),'b') | ||
| + | plot(t,Y(3,:),'g') | ||
| + | hold off | ||
| + | figure (2) | ||
| + | plot(Y(1,:),Y(2,:)) | ||
| + | legend('Trayectoria presas, primer tipo depredadores','Location','best') | ||
| + | figure (3) | ||
| + | plot(Y(1,:),Y(3,:)) | ||
| + | legend('Trayectoria presas, segundo tipo depredadores','Location','best') | ||
| + | figure(4) | ||
| + | plot(Y(2,:),Y(3,:)) | ||
| + | legend('Trayectoria primer tipo depredadores, segundo tipo depredadores','Location','best') | ||
| + | }} | ||
==Resolución del modelo mediante el método de Euler== | ==Resolución del modelo mediante el método de Euler== | ||
Revisión del 23:08 3 mar 2014
 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 | Modelo predador-presa. Grupo 16 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2013-14 |
| Autores | Javier Díez Olaya 121 Javier Lozano Aragoneses 248 Enrique Martínez Mur 271 Begoña Bigeriego Alvarez 637 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Interpretación
2 Resolución del modelo mediante el método de Euler modificado
% Datos del problema
t0=0;
tf=300;
%Coeficientes
A1=0.4; A2=0.3; A3=0.35;
B1=0.3; B2=0.05;
D=0.1;
C1=0.28; C2=0.045;
%Condiciones iniciales
p0=3.5;
d0=1;
e0=1.2;
Y0=[p0 d0 e0]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización
Y(:,1)=Y0;
YY=Y0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*YY(1)-A2*YY(1)*YY(2)-A3*YY(1)*YY(3);-B1*YY(2)+B2*YY(1)*YY(2)-D*YY(2)*YY(3);-C1*YY(3)+C2*YY(1)*YY(3)-D*YY(2)*YY(3)];
% Calculamos k2=f(t_n+h,y_n+h*k1)
k2=[A1*(YY(1)+h*k1(1))-A2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-A3*(YY(1)+h*k1(1))*(YY(3)+h*k1(3));-B1*(YY(2)+h*k1(2))+B2*(YY(1)+h*k1(1))*(YY(2)+h*k1(2))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3));-C1*(YY(3)+h*k1(3))+C2*(YY(1)+h*k1(1))*(YY(3)+h*k1(3))-D*(YY(2)+h*k1(2))*(YY(3)+h*k1(3))];
YY=YY+(h/2)*(k1+k2);
Y(:,n+1)=YY;
end
figure (1)
hold on
plot(t,Y(1,:),'r')
plot(t,Y(2,:),'b')
plot(t,Y(3,:),'g')
hold off
figure (2)
plot(Y(1,:),Y(2,:))
legend('Trayectoria presas, primer tipo depredadores','Location','best')
figure (3)
plot(Y(1,:),Y(3,:))
legend('Trayectoria presas, segundo tipo depredadores','Location','best')
figure(4)
plot(Y(2,:),Y(3,:))
legend('Trayectoria primer tipo depredadores, segundo tipo depredadores','Location','best')
3 Resolución del modelo mediante el método de Euler
t0=0;
tf=300;
%Variables
A1=0.5; A2=0.25; A3=0.3;
B1=0.4; B2=0.07;
D=0.05;
C1=0.38; C2=0.045;
%Condiciones iniciales
p0=3.5;
d0=1;
e0=0.2;
X0=[p0 d0 e0]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización
X(:,1)=X0;
XX=X0;
% Iteraciones
for n=1:N
XX=XX+h*[A1*XX(1)-A2*XX(1)*XX(2)-A3*XX(1)*XX(3);-B1*XX(2)+B2*XX(1)*XX(2)-D*XX(2)*XX(3);-C1*XX(3)+C2*XX(1)*XX(3)-D*XX(2)*XX(3)];
X(:,1+n)=XX;
end
figure (5)
hold on
plot(t,X(1,:),'r')
plot(t,X(2,:),'b')
plot(t,X(3,:),'k')
hold off
pm=max(X(1,:));
d1m=max(X(2,:));
d2m=max(X(3,:));
figure (6)
plot(X(1,:),X(2,:))
legend('Trayectoria presas, primer tipo depredadores','Location','best')
figure (7)
plot(X(1,:),X(3,:))
legend('Trayectoria presas, segundo tipo depredadores','Location','best')
4 Resolución del modelo mediante el método de Runge-Kutta de cuarto orden
% Datos del problema con los primeros valores iniciales
t0=0;
tf=500;
A1=0.4;
A2=0.3;
A3=0.4;
B1=0.37;
B2=0.05;
D=0.1;
C1=0.28;
C2=0.07;
p01=3.5;
d01=0.001;
e01=0.0002;
Y0=[p01 d01 e01]';
% Datos discretización
h=0.1;
N=(tf-t0)/h;
% Vectores de tiempo y soluciones app
t=t0:h:tf;
% Inicialización
Y(:,1)=Y0;
YY=Y0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*YY(1)-A2*YY(1)*YY(2)-A3*YY(1)*YY(3);-B1*YY(2)+B2*YY(1)*YY(2)-D*YY(2)*YY(3);-C1*YY(3)+C2*YY(1)*YY(3)-D*YY(2)*YY(3)];
% Calculamos k2=f(t_n+h/2,y_n+(1/2)*h*k1)
k2=[A1*(YY(1)+(h/2)*k1(1))-A2*(YY(1)+(h/2)*k1(1))*(YY(2)+(h/2)*k1(2))-A3*(YY(1)+(h/2)*k1(1))*(YY(3)+(h/2)*k1(3));-B1*(YY(2)+(h/2)*k1(2))+B2*(YY(1)+(h/2)*k1(1))*(YY(2)+(h/2)*k1(2))-D*(YY(2)+(h/2)*k1(2))*(YY(3)+(h/2)*k1(3));-C1*(YY(3)+(h/2)*k1(3))+C2*(YY(1)+(h/2)*k1(1))*(YY(3)+(h/2)*k1(3))-D*(YY(2)+(h/2)*k1(2))*(YY(3)+(h/2)*k1(3))];
% Calculamos k3=f(t_n+h/2,y_n+(1/2)*h*k2)
k3=[A1*(YY(1)+(h/2)*k2(1))-A2*(YY(1)+(h/2)*k2(1))*(YY(2)+(h/2)*k2(2))-A3*(YY(1)+(h/2)*k2(1))*(YY(3)+(h/2)*k2(3));-B1*(YY(2)+(h/2)*k2(2))+B2*(YY(1)+(h/2)*k2(1))*(YY(2)+(h/2)*k2(2))-D*(YY(2)+(h/2)*k2(2))*(YY(3)+(h/2)*k2(3));-C1*(YY(3)+(h/2)*k2(3))+C2*(YY(1)+(h/2)*k2(1))*(YY(3)+(h/2)*k2(3))-D*(YY(2)+(h/2)*k2(2))*(YY(3)+(h/2)*k2(3))];
% Calculamos k4=f(t_n+h,y_n+h*k3)
k4=[A1*(YY(1)+h*k3(1))-A2*(YY(1)+h*k3(1))*(YY(2)+h*k3(2))-A3*(YY(1)+h*k3(1))*(YY(3)+h*k3(3));-B1*(YY(2)+h*k3(2))+B2*(YY(1)+h*k3(1))*(YY(2)+h*k3(2))-D*(YY(2)+h*k3(2))*(YY(3)+h*k3(3));-C1*(YY(3)+h*k3(3))+C2*(YY(1)+h*k3(1))*(YY(3)+h*k3(3))-D*(YY(2)+h*k3(2))*(YY(3)+h*k3(3))];
YY=YY+(h/6)*(k1+2*k2+2*k3+k4);
Y(:,n+1)=YY;
end
figure (8)
hold on
plot(t,Y(1,:),'r','LineWidth',1.4)
plot(t,Y(2,:),'b','LineWidth',1.4)
plot(t,Y(3,:),'g','LineWidth',1.4)
hold off
% Datos del problema con los segundos valores iniciales
p02=3.5;
d02=0.00001;
e02=0.2;
Z0=[p02 d02 e02]';
% Inicialización
Z(:,1)=Z0;
ZZ=Z0;
for n=1:N
% Calculamos k1=f(t_n,y_n)
k1=[A1*ZZ(1)-A2*ZZ(1)*ZZ(2)-A3*ZZ(1)*ZZ(3);-B1*ZZ(2)+B2*ZZ(1)*ZZ(2)-D*ZZ(2)*ZZ(3);-C1*ZZ(3)+C2*ZZ(1)*ZZ(3)-D*ZZ(2)*ZZ(3)];
% Calculamos k2=f(t_n+h/2,y_n+(1/2)*h*k1)
k2=[A1*(ZZ(1)+(h/2)*k1(1))-A2*(ZZ(1)+(h/2)*k1(1))*(ZZ(2)+(h/2)*k1(2))-A3*(ZZ(1)+(h/2)*k1(1))*(ZZ(3)+(h/2)*k1(3));-B1*(ZZ(2)+(h/2)*k1(2))+B2*(ZZ(1)+(h/2)*k1(1))*(ZZ(2)+(h/2)*k1(2))-D*(ZZ(2)+(h/2)*k1(2))*(ZZ(3)+(h/2)*k1(3));-C1*(ZZ(3)+(h/2)*k1(3))+C2*(ZZ(1)+(h/2)*k1(1))*(ZZ(3)+(h/2)*k1(3))-D*(ZZ(2)+(h/2)*k1(2))*(ZZ(3)+(h/2)*k1(3))];
% Calculamos k3=f(t_n+h/2,y_n+(1/2)*h*k2)
k3=[A1*(ZZ(1)+(h/2)*k2(1))-A2*(ZZ(1)+(h/2)*k2(1))*(ZZ(2)+(h/2)*k2(2))-A3*(ZZ(1)+(h/2)*k2(1))*(ZZ(3)+(h/2)*k2(3));-B1*(ZZ(2)+(h/2)*k2(2))+B2*(ZZ(1)+(h/2)*k2(1))*(ZZ(2)+(h/2)*k2(2))-D*(ZZ(2)+(h/2)*k2(2))*(ZZ(3)+(h/2)*k2(3));-C1*(ZZ(3)+(h/2)*k2(3))+C2*(ZZ(1)+(h/2)*k2(1))*(ZZ(3)+(h/2)*k2(3))-D*(ZZ(2)+(h/2)*k2(2))*(ZZ(3)+(h/2)*k2(3))];
% Calculamos k4=f(t_n+h,y_n+h*k3)
k4=[A1*(ZZ(1)+h*k3(1))-A2*(ZZ(1)+h*k3(1))*(ZZ(2)+h*k3(2))-A3*(ZZ(1)+h*k3(1))*(ZZ(3)+h*k3(3));-B1*(ZZ(2)+h*k3(2))+B2*(ZZ(1)+h*k3(1))*(ZZ(2)+h*k3(2))-D*(ZZ(2)+h*k3(2))*(ZZ(3)+h*k3(3));-C1*(ZZ(3)+h*k3(3))+C2*(ZZ(1)+h*k3(1))*(ZZ(3)+h*k3(3))-D*(ZZ(2)+h*k3(2))*(ZZ(3)+h*k3(3))];
ZZ=ZZ+(h/6)*(k1+2*k2+2*k3+k4);
Z(:,n+1)=ZZ;
end
figure (9)
hold on
plot(t,Z(1,:),'r','LineWidth',1.4)
plot(t,Z(2,:),'b','LineWidth',1.4)
plot(t,Z(3,:),'g','LineWidth',1.4)
hold off
