MateWiki - Contribuciones del usuario [es]

T4. Modelos Epidemiológicos

2015-03-06T12:04:11Z

Luisld: /* Dificultad Método del Trapecio */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
t=t0:h:tN;
y=y0;
S(1)=y(1);
I(1)=y(2);
%Parámetros
a=0.003;
b=0.3;
c=0.01;
%Euler
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(1)
hold on
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off
%R-k
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
hold on
figure(2)
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1)=0.003/(1+t(i));
k1=-a(i)*S(i)*I(i);
k2=-a(i)*S(i)*I(i)+k1*h;
S(i+1)=S(i)+h/2*(k1+k2);
l1=a(i)*S(i)*I(i)-I(i)*(b+c);
l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T12:03:59Z

Luisld: /* Dificultad Método del Trapecio */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
t=t0:h:tN;
y=y0;
S(1)=y(1);
I(1)=y(2);
%Parámetros
a=0.003;
b=0.3;
c=0.01;
%Euler
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(1)
hold on
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off
%R-k
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(2)
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1)=0.003/(1+t(i));
k1=-a(i)*S(i)*I(i);
k2=-a(i)*S(i)*I(i)+k1*h;
S(i+1)=S(i)+h/2*(k1+k2);
l1=a(i)*S(i)*I(i)-I(i)*(b+c);
l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:57:57Z

Luisld: /* Dificultad Método del Trapecio */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
%Parámetros
a=0.003;
b=0.3;
c=0.01;
%Euler
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(1)
hold on
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off
%R-k
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(2)
plot(t,S,'g')
plot(t,I,'r')
legend('Susceptible','Infectados','Location','best')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1)=0.003/(1+t(i));
k1=-a(i)*S(i)*I(i);
k2=-a(i)*S(i)*I(i)+k1*h;
S(i+1)=S(i)+h/2*(k1+k2);
l1=a(i)*S(i)*I(i)-I(i)*(b+c);
l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:50:14Z

Luisld: /* Dificultad Método del Trapecio */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
%Parámetros
a=0.003;
b=0.3;
c=0.01;
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1)=0.003/(1+t(i));
k1=-a(i)*S(i)*I(i);
k2=-a(i)*S(i)*I(i)+k1*h;
S(i+1)=S(i)+h/2*(k1+k2);
l1=a(i)*S(i)*I(i)-I(i)*(b+c);
l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:48:05Z

Luisld: /* Método Heun, Variable a dependiente de t, a(t) */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1)=0.003/(1+t(i));
k1=-a(i)*S(i)*I(i);
k2=-a(i)*S(i)*I(i)+k1*h;
S(i+1)=S(i)+h/2*(k1+k2);
l1=a(i)*S(i)*I(i)-I(i)*(b+c);
l2 =(a(i)*S(i)*I(i))-I(i)*(b+c)+l1*h;
I(i+1)=I(i)+h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:46:58Z

Luisld: /* Método Heun, Variable a dependiente de t, a(t) */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
k1 = -a(i)*S(i)*I(i);
k2 = -a(i)*S(i)*I(i) + k1*h;
S(i+1) = S(i) + h/2*(k1+k2);
l1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
l2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + l1*h;
I(i+1) = I(i) + h/2*(l1+l2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

Archivo:Apt5b.jpg

2015-03-06T11:34:32Z

Luisld: 800 20 0.1 rk

800 20 0.1 rk

Archivo:Apt5a.jpg

2015-03-06T11:34:05Z

Luisld: Euler con 800 y 20 0.1

Euler con 800 y 20 0.1

T4. Modelos Epidemiológicos

2015-03-06T11:26:46Z

Luisld: /* S0=10000 y I0=40 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>12680.8</math> y se va a producir transcurridas <math>12</math> horas.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>9521.1</math> y se va a producir transcurridas <math>8</math> horas y <math>24</math> minutos.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001. El máximo de enfermos esperado en este caso es de <math>9469.6</math> y se va a producir transcurridas <math>8</math> horas y <math>12</math> minutos.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.El máximo de enfermos esperado en este caso es de <math>9464.7</math> y se va a producir transcurridas <math>8</math> horas y <math>11</math> minutos.]]

El máximo de enfermos esperado será en esta ocasión de <math>9465</math> y se va a producir al cabo de <math>8</math> horas y <math>11</math> minutos.

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:22:51Z

Luisld: /* S0=800, I0=20 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.]]

El maximo de enfermos esperado y cuando se va a producir es................

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:21:52Z

Luisld: /* S0=800, I0=20 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir transcurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir transcurridos <math>2</math> días y <math>17</math> horas.]]

Por tanto en esta primera simulación de poblaciones observamos que cuando reducimos el paso <math>h</math> obtenemos mayor precisión en la población así pues podemos indicar que el máximo alcanzado será entonces de <math>505</math>, obtenido a los <math>2</math> días y <math>17</math> horas aproximadamente.

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.]]

El maximo de enfermos esperado y cuando se va a producir es................

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:18:41Z

Luisld: /* S0=800, I0=20 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de <math>517.44</math> y se va a producir trancurridos <math>2</math> días y <math>19</math> horas.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01. El máximo de enfermos esperado en este caso es de <math>506.3</math> y se va a producir trancurridos <math>2</math> días y <math>17</math> horas.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.El máximo de enfermos esperado en este caso es de <math>505.3</math> y se va a producir trancurridos <math>2</math> días y <math>17</math> horas.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001. El máximo de enfermos esperado en este caso es de <math>505.2</math> y se va a producir trancurridos <math>2</math> días y <math>17</math> horas.]]

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.]]

El maximo de enfermos esperado y cuando se va a producir es................

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:15:34Z

Luisld: /* S0=800, I0=20 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1. El máximo de enfermos esperado en este caso es de y se va a producir trancurridos]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001.]]

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.]]

El maximo de enfermos esperado y cuando se va a producir es................

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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

T4. Modelos Epidemiológicos

2015-03-06T11:06:21Z

Luisld: /* S=0 */

{{ TrabajoED | Modelos Epidemiológicos, desarrollo de una epidemia. Grupo 28 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Luis Doreste Henríquez 127 Luis López Díaz 239 Enrique Martínez Mur 271}}

== Interpretación de los diferentes parámetros ==
En el desarrollo de la epidemia encontramos dos tipos distintos de población que influyen a la evolución de la misma y que se rigen por las siguientes ecuaciones

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>

Identificando primero las variables y luego los parámetros diferenciamos:
* T: tiempo
* S: población que puede ser infectada por la enfermedad. Parámetro en función del tiempo.
* I: población infectada por la enfermedad. Parámetro en función del tiempo.
* a: interacción entre el número de individuos infectados y no infectados.
* b: fallecimientos a causa de la enfermedad.
* c: población que consigue recuperarse de la enfermedad.

Además, en el desarrollo de la epidemia consideraremos que la tasa de cambio depende tanto de los fallecimientos como de la población curada y que la tasa de población susceptible a ser infectada es proporcional a la interacción entre ambas poblaciones.

Si queremos entender las ecuaciones por las que se rige nuestra epidemia:

:<math>\frac{\mathrm{d} S}{\mathrm{d} T}=-aSI </math>
Representa la exposición a la que se encuentra la población frente a la epidemia, es decir, la variación de la población susceptible a ser infectada, <math>S</math>, a causa de la interacción directa de ésta con la población ya infectada, <math>I</math>, en función del paramtro <math>a</math> antes descrito. Así observamos que es una función no lineal (producto de las variables <math>S</math> y <math>I</math>) y decreciente en el tiempo.

:<math>\frac{\mathrm{d} I}{\mathrm{d} T}=aSI-bI-cI </math>
En este caso se representa la variación de la población infectada con un primer término similar a la primera ecuación pero de signo contrario lo que significa que este término hace crecer el número de infectados por la epidemia y otro posterior ya únicamente dependiente de la población infectada, <math>I</math>, que esta acompañado por <math>-(b+c)</math> que representas las muertes y curas a causa de la enfermedad.

==Método de Euler y Trapecio ==
Consideramos que los parámetros del problema en el desarrollo de los metodos de Euler y del Trapecio son
<math>I0 = 2000</math>, <math>b = 0.3</math> y <math>c = 0.01</math> y lo resolveremos en ambos casos para una discretización <math>h = 0.1</math>
===S=0===

{{matlab|codigo=
%Apartado2
clear all
t0=0; tN=40;
I0=2000;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
I=zeros(1,N+1);
I(1)=I0;
for i=1:N
I(i+1)=I(i)-h*(0.3+0.001)*I(i);
end
%Gráficas
[t',I'];
[minimo,indice]=min(abs(I-500))
f=I(indice)
p=t(find(I==(500-minimo)))
hold on
plot(p,I,'r')
plot(t,f,'b')
plot(t,I)
hold off
sprintf('Los 500 enfermos se alcanza al cabo de %d días', p)
}}
[[Archivo:Apt2.jpg|300px|thumb|centro|Evolución de la población de los Infectados con ausencia de población susceptible a ser infectada. Se observa claramente la desaparición de la población al completo con el paso del tiempo.]]
La población de infectados se reduce hasta un cuarto del total inicial, siendo ese valor igual a <math>500</math>, en un tiempo en torno a <math>5</math> segundos, información sacada directamente de la interpretación de la gráfica aportada.
<math>5</math>

===S=100===
En este segundo caso necesitamos también el valor del parametro <math>a=0.003</math>.

{{matlab|codigo=
%Apartado 3
clear all
t0=0; tN=40; %Intervalo t
S0=100;
I0=2000;
%Valores iniciales I y S
h=0.1;
t=t0:h:tN;
a=0.003; %Dato problema
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S0;I(n)]+h*[-a*S0*I(n);a*S0*I(n)-(b+c)*I(n)];
S(n+1)=S0;
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
}}
 
[[Archivo:Apt3.jpg|300px|thumb|left|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=100). Se observa que la población susceptible se mantiene constante en el tiempo mientras que la población de infectados desciende linealmente en el tiempo.]]
[[Archivo:Apt3b.jpg|300px|thumb|centro|Evolución de la población de los Infectados en presencia de población susceptible a ser infectada (S=200). En este caso se observa que la población susceptible acaba por desaparecer mientras que la poblacion infectada crece de forma exponencial.]]
 

== Euler Completo ==
En el modelo de Euler completo consideramos ahora los datos iniciales <math>S0=800</math>, <math>I0=20</math> y <math>S0=10000</math> y <math>I0=40</math> para tiempos menores de <math>40</math> días y discretización temporal <math>h=0.1, h=0.01,h=0.001,h=0.0001</math>.

{{matlab|codigo=
%Apartado 4
clear all
t0=0; tN=40; %Intervalo t
S0=input('N población susceptible: ');
I0=input('N infectados: ');
%Valores iniciales I y S
h=input('Valor del paso: ');
t=t0:h:tN;
a=0.003;
b=0.3;
c=0.01;
N=(tN-t0)/h;
S(1)=S0;
I(1)=I0;
%Resolucion forma matricial
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
hold on
figure(1)
plot(t,S,'r')
plot(t,I,'b')
hold off
maxi=max(I);
pos=find(I==maxi);
tiempo=t(pos);
t1=sprintf('El número máximo de infectados será %d',maxi);
t2=sprintf(' y se dará transcurridos %d días',tiempo);
t3=[t1 t2];
disp(t3)

}}
=== S0=800, I0=20 ===

[[Archivo:Apt4 1 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 1 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 1 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 1 4.jpg|300px|thumb|centro|h=0.0001.]]

=== S0=10000 y I0=40===
[[Archivo:Apt4 2 1.jpg|300px|thumb|left|h=0.1.]]
[[Archivo:Apt4 2 2.jpg|300px|thumb|centro|h=0.01.]]

[[Archivo:Apt4 2 3.jpg|300px|thumb|left|h=0.001.]]
[[Archivo:Apt4 2 4.jpg|300px|thumb|centro|h=0.0001.]]

El maximo de enfermos esperado y cuando se va a producir es................

== Runge-Kutta 4º Orden ==

=== Comparativa con Euler ===
Si realizamos una comparativa entre ambos métodos, encontramos la similitud de que ambos métodos están basados en iteraciones. El método de Euler busca la pendiente de la recta tangente empleando un paso muy pequeño. Por su parte, el método de Runge Kutta emplea cuatro parámetros, que son las pendientes de las cuatro tangentes a la curva integral, y empleando un paso mayor. Podemos concluir, por tanto, que el método de Runge Kutta supone una mejora con respecto al de Euler.

=== Dificultad Método del Trapecio ===

Cuando se emplea un método implícito, como podría ser el método trapezoidal, para resolver un problema, la mayor dificultad que podemos encontrar es el hecho de que deberemos despejar de una forma manual aquellos parámetros necesarios para poder posteriormente programar la solución. Puede que, en algún caso, el proceso sea bastante complicado, llegando incluso a ser imposible despejar.

{{matlab|codigo=

%Apartado5
clear all
t0=0; tN=30; %Intervalo t
y1=input('N población susceptible: ');
y2=input('N infectados de inicio: ');
y0=[y1;y2]; %Valores iniciales S e I
h=input('valor del paso: ');
N=(tN-t0)/h;
y=y0;
S(1)=y(1);
I(1)=y(2);
a=0.003; b=0.3; c=0.01; % Valores de los parámetros
%Resolución Runge Kutta
for n=1:N
k1=[-a*y(1)*y(2);a*y(1)*y(2)-(b+c)*y(2)];
k2=[-a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));a*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2)-(b+c)*(y(2)+1/2*h*k1(2)))];
k3=[-a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));a*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2)-(b+c)*(y(2)+1/2*h*k2(2)))];
k4=[-a*(y(1)+h*k3(1))*(y(2)+h*k3(2));a*(y(1)+h*k3(1))*(y(2)+h*k3(2)-(b+c)*(y(2)+h*k3(2)))];
y=y+h/6*(k1+2*k2+2*k3+k4);
S(n+1)=y(1);
I(n+1)=y(2);
end
t=t0:h:tN;
hold on
figure(1)
plot(t,S,'b')
plot(t,I,'r')
hold off
for n=1:N
A=[S(n);I(n)]+h*[-a*S(n)*I(n);a*S(n)*I(n)-(b+c)*I(n)];
S(n+1)=A(1);
I(n+1)=A(2);
end
%Gráfica
figure(2)
hold on
plot(t,S,'b')
plot(t,I,'r')
hold off

}}

GRAFICAS FALTA SACAR

== Método Heun, Variable a dependiente de t, a(t) ==

{{matlab|codigo=
%Apartado6
t0=0; tN=40;
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
S=zeros(1,N+1);
I=zeros(1,N+1);
S(1)=1600;
I(1)=40;
a=zeros(1,N+1);
a(1)=0.003/(1+t0); %a(t)
b=0.3; c=0.01;
for i=1:N
a(i+1) = 0.003/(1+t(i));
K1 = -a(i)*S(i)*I(i);
K2 = -a(i)*S(i)*I(i) + K1*h;
S(i+1) = S(i) + h/2*(K1+K2);
L1 = (a(i)*S(i)*I(i)) -I(i)*(b+c);
L2 = (a(i)*S(i)*I(i)) -I(i)*(b+c) + L1*h;
I(i+1) = I(i) + h/2*(L1+L2);
end
%Gráficas
hold on
plot (t,S,'b')
plot (t,I,'r')
hold off
}}

[[Archivo:Apt6.jpg|300px|thumb|centro|]]

El coeficiente de variable a(t) lo interpretamos.....

[[Categoría:Ecuaciones Diferentiales]]

[[Categoría:ED14/15]]

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