MateWiki - Contribuciones del usuario [es]

Modelos Epidemológicos Grupo 4-C

2015-03-11T15:25:39Z

Rojasrivero: /* Método de Heun */

{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
Suponemos que :
# La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
# La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con '''S=0''', es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.</math>: <math>b=0.3 \\ c=0.01</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50;
% consideramos el intervalo de tiempo [0,50] ya que el numero
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora '''S=100''' constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Repetimos el proceso para '''S=200''' constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)
%a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}
== Modelo Completo==
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
e I(0)=40 S(0)=10000.
En el primer caso utilizamos el siguiente código MatLab:
[[Archivo:graficafer1.png|400px|thumb|right|Modelo Completo]]

{{matlab|codigo=

% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
1.I=517.45
2.I=506.37
3.I=505.3
4.I=505.19
En el segundo caso utilizamos el código MatLab:
[[Archivo:graficafer2.png|400px|thumb|right|Modelo Completo]]
{{matlab|codigo=
% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=40;
S0=10000;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); % también N=length(t)-1;
N2=round((tN-t0)/h2); % también N=length(t)-1;
N3=round((tN-t0)/h3); % también N=length(t)-1;
N4=round((tN-t0)/h4); % también N=length(t)-1;
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
1.I=1.26+e004
2.I=9521.1
3.I=9469.6
4.I=9464.7
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta.
Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
{{matlab|codigo=
t0=0;
tN=25;
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
% RK4
K1=a*S.*r(i)-b.*r(i)-c.*r(i);
K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4)
end

%sacamos tabla de resultados
[t',I',r']
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');
}}

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. La interpretación más lógica para que la proporción de gente susceptible que se relaciona con los infectados vaya disminiyendo con el tiempo es que, debido al conocimiento de la epidemia, la gente empiece a tomar medidas para no estar en contacto con los infectados.Con lo que el sistema nos quedaría asi:

<math> \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.</math>
En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectado.

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
{{matlab|codigo=
%sistemas no lineales, no puedo usar matrices
clear ; clc;
t0=0;
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
a(i+1)=0.003/(1+t(i));
K1S=-a(i).*S(i).*I(i);
K2S=-a(i)*S(i).*I(i)+h*K1S;
S(i+1)=S(i)+(h/2).*(K1S+K2S);
K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off
}}

Es destacable ver que con un valor para la a constante y con estos valores iniciales parecidos, el modelo no es válido, pero con estas consideraciones ya si nos queda un resultado que se ajusta más a la realidad.

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]

{{matlab|codigo=

clear all;
clc;
h=0.0001;
t0=0;
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
for i=1:N
K1S=a*S(i)*I(i);
K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
I(i+1)=I(i)+(h/2)*(K1I+K2I);
S(i+1)=S(i)+(h/2)*(K1S+K2S);
[v,w]=max(I);
end
c=[c;[a,t(w)]];
end

}}

[[Categoría:Ecuaciones Diferenciales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:24:30Z

Rojasrivero:

{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
Suponemos que :
# La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
# La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con '''S=0''', es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.</math>: <math>b=0.3 \\ c=0.01</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50;
% consideramos el intervalo de tiempo [0,50] ya que el numero
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora '''S=100''' constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Repetimos el proceso para '''S=200''' constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)
%a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}
== Modelo Completo==
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
e I(0)=40 S(0)=10000.
En el primer caso utilizamos el siguiente código MatLab:
[[Archivo:graficafer1.png|400px|thumb|right|Modelo Completo]]

{{matlab|codigo=

% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
1.I=517.45
2.I=506.37
3.I=505.3
4.I=505.19
En el segundo caso utilizamos el código MatLab:
[[Archivo:graficafer2.png|400px|thumb|right|Modelo Completo]]
{{matlab|codigo=
% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=40;
S0=10000;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); % también N=length(t)-1;
N2=round((tN-t0)/h2); % también N=length(t)-1;
N3=round((tN-t0)/h3); % también N=length(t)-1;
N4=round((tN-t0)/h4); % también N=length(t)-1;
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
1.I=1.26+e004
2.I=9521.1
3.I=9469.6
4.I=9464.7
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta.
Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
{{matlab|codigo=
t0=0;
tN=25;
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
% RK4
K1=a*S.*r(i)-b.*r(i)-c.*r(i);
K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4)
end

%sacamos tabla de resultados
[t',I',r']
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');
}}

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math> \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.</math>
En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectados

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
{{matlab|codigo=
%sistemas no lineales, no puedo usar matrices
clear ; clc;
t0=0;
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
a(i+1)=0.003/(1+t(i));
K1S=-a(i).*S(i).*I(i);
K2S=-a(i)*S(i).*I(i)+h*K1S;
S(i+1)=S(i)+(h/2).*(K1S+K2S);
K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off
}}

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]

{{matlab|codigo=

% nos dicen que Imax se alcanza en t=5 en una experiencia
% Resolver el sistema para a(0.0005,0.002) y ver cual se ajusta a la exp
%usar un paso de 10^-4 y usar bucle para el valor max
%nos quedamos con el valor de a que da un maximo mas cerca de t=5
clear; clc;
h=0.0001;
t0=0;
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
for i=1:N
K1S=a*S(i)*I(i);
K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
I(i+1)=I(i)+(h/2)*(K1I+K2I);
S(i+1)=S(i)+(h/2)*(K1S+K2S);
[v,w]=max(I);
end
c=[c;[a,t(w)]];
end

}}

[[Categoría:Ecuaciones Diferenciales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Archivo:Graficafer2.png

2015-03-05T10:22:44Z

Rojasrivero:

Archivo:Graficafer1.png

2015-03-05T10:22:30Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:20:15Z

Rojasrivero:

{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
Suponemos que :
# La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
# La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con '''S=0''', es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.</math>: <math>b=0.3 \\ c=0.01</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50;
% consideramos el intervalo de tiempo [0,50] ya que el numero
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora '''S=100''' constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Repetimos el proceso para '''S=200''' constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)
%a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}
== Modelo Completo==
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
e I(0)=40 S(0)=10000.
En el primer caso utilizamos el siguiente código MatLab:
{{matlab|codigo=

% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
1.I=517.45
2.I=506.37
3.I=505.3
4.I=505.19
En el segundo caso utilizamos el código MatLab:
{{matlab|codigo=
% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=40;
S0=10000;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); % también N=length(t)-1;
N2=round((tN-t0)/h2); % también N=length(t)-1;
N3=round((tN-t0)/h3); % también N=length(t)-1;
N4=round((tN-t0)/h4); % también N=length(t)-1;
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
1.I=1.26+e004
2.I=9521.1
3.I=9469.6
4.I=9464.7
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta.
Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
{{matlab|codigo=
t0=0;
tN=25;
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
% RK4
K1=a*S.*r(i)-b.*r(i)-c.*r(i);
K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4)
end

%sacamos tabla de resultados
[t',I',r']
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');
}}

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math> \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.</math>
En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectados

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
{{matlab|codigo=
%sistemas no lineales, no puedo usar matrices
clear ; clc;
t0=0;
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
a(i+1)=0.003/(1+t(i));
K1S=-a(i).*S(i).*I(i);
K2S=-a(i)*S(i).*I(i)+h*K1S;
S(i+1)=S(i)+(h/2).*(K1S+K2S);
K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off
}}

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]

{{matlab|codigo=

% nos dicen que Imax se alcanza en t=5 en una experiencia
% Resolver el sistema para a(0.0005,0.002) y ver cual se ajusta a la exp
%usar un paso de 10^-4 y usar bucle para el valor max
%nos quedamos con el valor de a que da un maximo mas cerca de t=5
clear; clc;
h=0.0001;
t0=0;
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
for i=1:N
K1S=a*S(i)*I(i);
K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
I(i+1)=I(i)+(h/2)*(K1I+K2I);
S(i+1)=S(i)+(h/2)*(K1S+K2S);
[v,w]=max(I);
end
c=[c;[a,t(w)]];
end

}}

[[Categoría:Ecuaciones Diferenciales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:18:24Z

Rojasrivero:

{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
Suponemos que :
# La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
# La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con '''S=0''', es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.</math>: <math>b=0.3 \\ c=0.01</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50;
% consideramos el intervalo de tiempo [0,50] ya que el numero
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora '''S=100''' constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Repetimos el proceso para '''S=200''' constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)
%a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}
== Modelo Completo==
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
e I(0)=40 S(0)=10000.
En el primer caso utilizamos el siguiente código MatLab:
{{matlab|codigo=

% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
1.I=517.45
2.I=506.37
3.I=505.3
4.I=505.19
En el segundo caso utilizamos el código MatLab:
(código de TRABAJO2.m y su gráfica)
Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
1.I=1.26+e004
2.I=9521.1
3.I=9469.6
4.I=9464.7
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta.
Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
{{matlab|codigo=
t0=0;
tN=25;
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
% RK4
K1=a*S.*r(i)-b.*r(i)-c.*r(i);
K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4)
end

%sacamos tabla de resultados
[t',I',r']
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');
}}

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math> \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.</math>
En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectados

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
{{matlab|codigo=
%sistemas no lineales, no puedo usar matrices
clear ; clc;
t0=0;
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
a(i+1)=0.003/(1+t(i));
K1S=-a(i).*S(i).*I(i);
K2S=-a(i)*S(i).*I(i)+h*K1S;
S(i+1)=S(i)+(h/2).*(K1S+K2S);
K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off
}}

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]

{{matlab|codigo=

% nos dicen que Imax se alcanza en t=5 en una experiencia
% Resolver el sistema para a(0.0005,0.002) y ver cual se ajusta a la exp
%usar un paso de 10^-4 y usar bucle para el valor max
%nos quedamos con el valor de a que da un maximo mas cerca de t=5
clear; clc;
h=0.0001;
t0=0;
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
for i=1:N
K1S=a*S(i)*I(i);
K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
I(i+1)=I(i)+(h/2)*(K1I+K2I);
S(i+1)=S(i)+(h/2)*(K1S+K2S);
[v,w]=max(I);
end
c=[c;[a,t(w)]];
end

}}

[[Categoría:Ecuaciones Diferenciales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:11:57Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:11:10Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:03:20Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:02:41Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T10:01:41Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:57:55Z

Rojasrivero: /* Método Runge-Kutta */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:54:31Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:54:05Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:53:07Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:47:51Z

Rojasrivero: /* Método Runge-Kutta */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:44:26Z

Rojasrivero: /* Método de Euler y Trapecio */

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:40:38Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:39:43Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:39:18Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:27:35Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:26:43Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-05T09:25:39Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T22:59:22Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T22:49:28Z

Rojasrivero: /* Método de Heun */

Archivo:Ejercicio69.png

2015-03-04T22:48:58Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T22:46:32Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T22:34:41Z

Rojasrivero:

Archivo:Bonethijoputa.png

2015-03-04T22:32:22Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T18:26:25Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-04T12:11:10Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T12:09:13Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T12:07:27Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:59:39Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:59:19Z

Rojasrivero: /* Método de Heun */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:58:40Z

Rojasrivero: /* Método de Euler y Trapecio */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:57:40Z

Rojasrivero: /* Método Runge-Kutta */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:57:03Z

Rojasrivero: /* Método de Euler y Trapecio */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:55:46Z

Rojasrivero: /* Método de Euler y Trapecio */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:55:10Z

Rojasrivero: /* Método de Euler y Trapecio */

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:53:33Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:48:50Z

Rojasrivero:

Archivo:Graficazorra3.png

2015-03-04T11:45:46Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:45:16Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:43:50Z

Rojasrivero: /* Método Runge-Kutta */

Archivo:Graficazorra2.png

2015-03-04T11:41:44Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:41:25Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:39:08Z

Rojasrivero:

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:38:27Z

Rojasrivero: /* Método de Heun */

<gallery>
Archivo:Ejemplo.jpg|Descripción1
Archivo:Ejemplo.jpg|Descripción2
</gallery>
{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math>\frac{dS}{dt}=-aSI</math>:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con S=0, es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\frac{dI}{dt}=-bI-cI;</math>: <math>b=0.3 \\ c=0.01 \\ I_{0}=2000</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero.
Según se interpreta en la gráfica(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50; % consideramos el intervalo de tiempo [0,50] ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); % Vector para cáculo del tiempo que tarda en reducirse I(i) a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora S=100 constante a lo largo del tiempo y a=0,003, hacemos uso de un codigo de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor sera el tiempo que tarde en erradicarse.

CODIGO Y GRAFICA

Repetimos el proceso para S=200 constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente. Por tanto se requiere el uso de un código Matlab.

CODIGO Y GRAFICA

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math>\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I
\frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I</math>

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente codigo:

Basandonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a.

Vemos que el que mas se ajusta más a la experiencia es el...
[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Modelos Epidemológicos Grupo 4-C

2015-03-04T11:36:08Z

Rojasrivero:

<gallery>
Archivo:Ejemplo.jpg|Descripción1
Archivo:Ejemplo.jpg|Descripción2
</gallery>
{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math>\frac{dS}{dt}=-aSI</math>:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con S=0, es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\frac{dI}{dt}=-bI-cI;</math>: <math>b=0.3 \\ c=0.01 \\ I_{0}=2000</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero.
Según se interpreta en la gráfica(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50; % consideramos el intervalo de tiempo [0,50] ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); % Vector para cáculo del tiempo que tarda en reducirse I(i) a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora S=100 constante a lo largo del tiempo y a=0,003, hacemos uso de un codigo de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor sera el tiempo que tarde en erradicarse.

CODIGO Y GRAFICA

Repetimos el proceso para S=200 constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente. Por tanto se requiere el uso de un código Matlab.

CODIGO Y GRAFICA

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math>\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I
\frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I</math>

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente codigo:
[[Archivo:Graficazorra.png|600px|miniaturadeimagen|centro|Método de Heun]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1); % Vector para cáculo del tiempo que tarda en reducirse I(i) a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Basandonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a.

Vemos que el que mas se ajusta más a la experiencia es el...
[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]