El ejercicio representa un sistema de ecuaciones diferenciales ordinarias, que serán resueltas usando Runge-Kutta de 2do Orden.
De compararse con la curva de contagios de Covid-19 se tienen diferencias en la población recuperada, pues el modelo se considera no letal por lo que no se contabiliza el número de fallecidos.
El módelo es el más básico y permite cambiar por ejemplo la tasa de infección, y se ve los cambios en la curva de infectados. Se puede observar lo que se indicaba como objetivo de «aplanar la curva» al disminuir la población expuesta mediante cambiar la tasa de infección al exponer más o menos población al contagio por iteacción entre «suceptibles» e «infectados.
Desarrollo analítico
Las fórmulas para el algoritmo se identifican como:
binfecta = 1.4
grecupera = 1/4
# 3Eva_2020PAOI_T2 Modelo epidemiológico no letal# Sistemas EDO con Runge-Kutta de 2do Ordenimport numpy as np
defrungekutta2_fgw(f,g,w,t0,x0,y0,z0,h,muestras):
tamano = muestras +1
tabla = np.zeros(shape=(tamano,4),dtype=float)
tabla[0] = [t0,x0,y0,z0]
ti = t0
xi = x0
yi = y0
zi = z0
for i inrange(1,tamano,1):
K1x = h * f(ti,xi,yi,zi)
K1y = h * g(ti,xi,yi,zi)
K1z = h * w(ti,xi,yi,zi)
K2x = h * f(ti+h, xi + K1x, yi+K1y, zi +K1z)
K2y = h * g(ti+h, xi + K1x, yi+K1y, zi +K1z)
K2z = h * w(ti+h, xi + K1x, yi+K1y, zi +K1z)
xi = xi + (1/2)*(K1x+K2x)
yi = yi + (1/2)*(K1y+K2y)
zi = zi + (1/2)*(K1z+K2z)
ti = ti + h
tabla[i] = [ti,xi,yi,zi]
tabla = np.array(tabla)
return(tabla)
# Programa# Parámetros de las ecuaciones
binfecta = 1.4
grecupera = 1/4
# Ecuaciones
f = lambda t,S,I,R : -binfecta*S*I
g = lambda t,S,I,R : binfecta*S*I - grecupera*I
w = lambda t,S,I,R : grecupera*I
# Condiciones iniciales
t0 = 0
S0 = 1.0
I0 = 0.001
R0 = 0.0
# parámetros del algoritmo
h = 1.0
muestras = 51
# PROCEDIMIENTO
tabla = rungekutta2_fgw(f,g,w,t0,S0,I0,R0,h,muestras)
ti = tabla[:,0]
Si = tabla[:,1]
Ii = tabla[:,2]
Ri = tabla[:,3]
# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, Si, Ii, Ri]')
print(tabla)
# Grafica tiempos vs poblaciónimport matplotlib.pyplot as plt
plt.plot(ti,Si, label='S')
plt.plot(ti,Ii, label='I')
plt.plot(ti,Ri, label='R')
plt.title('Modelo SIR')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()
Como en trazadores cúbicos, lo que se usa es un polinomio por cada tramo muestreado para una curva contínua, etc. se tiene que los polinomios deben tener valores iguales en los puntos el eje x = 0.4 y 0.6
Por lo que se evalua con los polinomios completos:
Valores que se usan en los extremos del polinomio S1(x) para crear un sistema de dos ecuaciones y determinar los valores de c y d, completando el polinomio.
la otra ecuación se podría obtener usando la propiedad que las primeras derivadas de los polinomios deben ser iguales en los puntos x=0,4 y x= 0.6
S’0(0.2) =S’1(0.2)
Tarea: Desarrollar la siguiente ecuación y resolver
Opción 2
Si no recuerda la propiedad anterior, puede optar por usar otros conceptos para aproximar el resultado.
Si para el tramo en que se busca el polinomio se puede retroceder un tamaño de paso x = 0.2 y evualuar usando S0(0.2), se obtiene otropunto de referencia para crear un polinomio que pase por los mismos puntos.
S0(0.2)=1+1.1186∗(0.2)+0.6938(0.2)3
Se aplica lo mismo para un tamaño de paso más adelante de x = 0.6 es x = 0.8m se evalua S2(0.8) y se tienen suficientes puntos para usar cualquier método de interpolación y determinar el polinomio para el tramo faltante.
que permite hacer una tabla de puntos, y usando por ejemplo el método de interpolación de Lagrange con x entre [0.2, 0.8] se obtiene otra forma del polinomio buscado:
A partir de la tabla del enunciado se realiza la tabla de diferencias finitas.
i
ti
fi
Δfi
Δ2fi
Δ3fi
Δ4fi
Δ5fi
1
0
0
32
-6
0
0
0
2
25
32
26
-6
0
0
3
50
58
20
-6
0
4
75
78
14
-6
5
100
92
8
6
125
100
Observando que a partir de la tercera diferencia finita los valores son cero, por lo que se usa la fórmula general de diferencias finitas divididas hasta el polinomio de grado 2.
# 3Eva_IIT2019_T1 Lanzamiento de Cohete# Tarea: Verificar tamaño de vectores# considerar puntos no equidistantes en eje timport numpy as np
import matplotlib.pyplot as plt
import sympy as sym
# INGRESO , Datos de prueba
ti = np.array([0.0, 25, 50, 75, 100, 125])
fi = np.array([0.0, 32, 58, 78, 92, 100])
# PROCEDIMIENTO# Tabla de diferencias finitas
titulo = ['i','ti','fi']
n = len(ti)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
ti = np.transpose([ti])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,ti,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
# Aumenta el título para cada columna
titulo.append('df'+str(c-2))
# calcula cada diferencia por fila
f = 0
while (f < diagonal):
tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
f = f+1
diagonal = diagonal - 1
c = c+1
# POLINOMIO con diferencias finitas# caso: puntos en eje t equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
t = sym.Symbol('t')
h = ti[1,0]-ti[0,0]
polinomio = fi[0,0]
for c inrange(1,n,1):
denominador = np.math.factorial(c)*(h**c)
factor = dfinita[0,c-1]/denominador
termino=1
for f inrange(0,c,1):
termino = termino*(t-ti[f])
polinomio = polinomio + termino*factor
# simplifica polinomio, multiplica los (t-ti)
polinomio = polinomio.expand()
# para evaluacion numérica
pt = sym.lambdify(t,polinomio)
# Puntos para la gráfica
a = np.min(ti)
b = np.max(ti)
muestras = 101
ti_p = np.linspace(a,b,muestras)
fi_p = pt(ti_p)
# SALIDAprint([titulo])
print(tabla)
print('polinomio:')
print(polinomio)
# Gráfica
plt.title('Interpolación polinómica')
plt.plot(ti,fi,'o', label = 'Puntos')
plt.plot(ti_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()
El resultado se interpreta mejor con una animación: (tarea)
Tarea: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación
Algorirmo en Python
# 3Eva_IT2019_T3 Difusión en sólidos# EDP parabólicas. método explícito,usando diferencias finitasimport numpy as np
import matplotlib.pyplot as plt
# INGRESO# Valores de frontera
Ta = 5
Tb = 0
T0 = 0
# longitud en x
a = 0
b = 0.1
# Constante K
K = 1/(1.6e-1)
# Tamaño de paso
dx = 0.02
dt = dx/100
# iteraciones en tiempo
n = 50
# PROCEDIMIENTO# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
m = len(xi)
ultimox = m-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)
# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0
u[ultimox,j] = Tb
# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb
# Calcula U para cada tiempo + dt
j = 0
whilenot(j>=ultimot):
u[0,j+1] = Ta
for i inrange(1,ultimox,1):
u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
u[m-1,j+1] = Tb
j=j+1
# SALIDAprint('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
# Gráfica
salto = int(n/10)
if (salto == 0):
salto = 1
for j inrange(0,n,salto):
vector = u[:,j]
plt.plot(xi,vector)
plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('phi[i,j]')
plt.title('Solución EDP parabólica')
plt.show()
La animación se complementa con lo mostrado en la sección de Unidades.
El ejercicio considera dos partes: interpolación e integración
a. Interpolación
Se requiere aproximar la función usando tres puntos. Para comprender la razón del método solicitado, se compara la función con dos interpolaciones:
a.1 Lagrange
a.2 Trazador cúbico sujeto
Observando la gráfica se aclara que en éste caso, una mejor aproximación se obtiene con el método de trazador cúbico sujeto. Motivo por lo que el tema tiene un peso de 40/100 puntos
