Ejercicio: 1Eva_2021PAOI_T3 Interpolar, modelo de contagios 2020
literal a
Los datos de los pacientes casos graves entre las semanas 11 a la 20, que son el intervalo donde será válido el polinomio de interpolación son:
| semana | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
| casos graves | 1503 | 3728 | 7154 | 6344 | 4417 | 3439 | 2791 | 2576 | 2290 | 2123 |
de los cuales solo se usarán los indicados en el literal a : 11,13,16,18,20.
xi0 = [ 9, 10, 11, 12, 13, 14,
15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26 ])
fi0 = [ 1435, 1645, 1503, 3728, 7154, 6344,
4417, 3439, 2791, 2576, 2290, 2123,
2023, 2067, 2163, 2120, 2125, 2224 ])
xi = [ 11, 13, 16, 18, 20]
fi = [1503, 7154, 3439, 2576, 2123]
Se observa que los datos estan ordenados en forma ascendente respecto a la variable independiente, tambien se determina que no se encuentran equidistantes entre si (13-11=2, 16-13=3). Por lo que se descarta usar el método de diferencias finitas avanzadas.
Los métodos que se podrían usar con puntos no equidistantes en el eje semanas serían el método de diferencias divididas de Newton o el método de Lagrange.
Seleccionando por ejemplo, Diferencias divididas de Newton, donde primero se realiza la tabla:
| xi | fi | f[x1,x0] | f[x2,x1,x0] | f[x3,x2,x1,x0] | f[x4,x3,x2,x1,x0] |
| 11 | 1503 | =(7154-1503) /(13-11) = 2825.5 | =(-1238.33-2835.5) /(16-11) = -812.76 | =(161.36-(-812.76)) /(18-11) = 139.16 | =(-15.73-139.16) /(20-11) = -17.21 |
| 13 | 7154 | (3439-7154) /(16-13) = -1238.33 | (-431.5-(-1238.33)) /(18-13) = 161.36 | (51.25-161.36) /(20-13)= -15.73 | —- |
| 16 | 3439 | (2576-3439) /(18-16) = -431.5 | (-226.5-(-431.5)) /(20-16) = 51.25 | —- | |
| 18 | 2576 | (2123-2576) /(20-18) = -226.5 | —- | ||
| 20 | 2123 | —- |
con lo que se puede contruir el polinomio usando las diferencias divididas para el intervalo dado:
[2825.5 -812.76 139.16 -17.21]p_4(x) = 1503 + 2825.5(x-11) - 812.76(x - 13)(x - 11) + 139.16(x - 16)(x - 13)(x - 11) - 17.21(x - 18)(x - 16) (x - 13) (x - 11)
Simplificando el algoritmo se tiene:
p_4(x) = - 1172995.28 + 298304.50 x - 27840.50x^2 + 1137.36x^3 - 17.21 x^4
literal b
El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.
| xi | fi | p4(x) | |error| |
|---|---|---|---|
| 11 | 1503 | 1503 | 0 |
| 12 | 3728 | 6110.96 | 2382.96 |
| 13 | 7154 | 7154 | 0 |
| 14 | 6344 | 6293.18 | 50.81 |
| 15 | 4417 | 4776.52 | 359.52 |
| 16 | 3439 | 3439 | 0 |
| 17 | 2791 | 2702.53 | 88.46 |
| 18 | 2576 | 2576 | 0 |
| 19 | 2290 | 2655.22 | 365.22 |
| 20 | 2123 | 2123 | 0 |
literal c
Podría aplicarse uno de varios criterios, lo importante por lo limitado del tiempo en la evaluación son las conclusiones y recomendaciones expresadas en el literal e, basadas en lo realizado en los literales c y d. Teniendo como opciones:
– cambiar uno de los puntos selecionados, mateniendo así el grado del polinomio
– aumentar el número de puntos usados para armar el polinomio con grado mayor
– dividir el intervalo en uno o mas segmentos, con el correspondiente número de polinomios.
Se desarrolla la opción de cambiar uno de los puntos seleccionados, usando para esta ocasión como repaso la interpolación de Lagrange. Para los puntos se usa el punto con mayor error de la tabla del literal anterior y se elimina el punto penúltimo, es decir se usa la semana 12 en lugar de la semana 18 de la siguiente forma:
xi = [ 11, 12, 13, 16, 20] fi = [1503, 3728, 7154, 3439, 2123]p_4(x) = 1503 \frac{(x-12)(x-13)(x-16)(x-20)}{(11-12)(11-13)(11-16)(11-20)} + 3728\frac{(x-11)(x-13)(x-16)(x-20)}{(12-11)(12-13)(12-16)(12-20)} + 7154\frac{(x-11)(x-12)(x-16)(x-20)}{(13-11)(13-12)(13-16)(13-20)} + 3439\frac{(x-11)(x-12)(x-13)(x-20)}{(16-11)(16-12)(16-13)(16-20)} + 2123\frac{(x-11)(x-12)(x-13)(x-16)}{(20-11)(20-12)(20-13)(20-16)}
Simplificando el polinomio:
p_4(x) = \frac{46927445}{21} - \frac{1655552687}{2520} x + \frac{715457663}{10080}x^2 - \frac{8393347}{2520} x^3 + \frac{577153}{10080} x^4literal d
El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.
| xi | fi | p4(x) | error |
| 11 | 1503 | 1503 | 0 |
| 12 | 3728 | 3728 | 0 |
| 13 | 7154 | 7154 | 0 |
| 14 | 6344 | 8974.01 | 2630.01 |
| 15 | 4417 | 7755.22 | 3338.22 |
| 16 | 3439 | 3439 | 0 |
| 17 | 2791 | -2659.13 | -5450.13 |
| 18 | 2576 | -7849.45 | -10425.45 |
| 19 | 2290 | -8068.1 | -10358.1 |
| 20 | 2123 | 2123 | 0 |
literal e
El cambio aplicado a los puntos usados en el modelo del polinomio disminuyó el error entre las semanas 11 a 13. Sin embargo la magnitud del error aumentó para las semanas posteriores a la 13, es decir aumentó la distorsión de la estimación y se recomienda realizar otras pruebas para mejorar el modelo aplicando los otros criterios para determinar el que tenga mejor desempeño respecto a la medida de error.
Intrucciones en Python
Literal a y b. Desarrollado a partir del algoritmo desarrollado en clases:
# Polinomio interpolación # Diferencias Divididas de Newton # Tarea: Verificar tamaño de vectores, # verificar puntos equidistantes en x import numpy as np import sympy as sym import matplotlib.pyplot as plt # INGRESO , Datos de prueba xi0 = np.array([ 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26 ]) fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344, 4417, 3439, 2791, 2576, 2290, 2123, 2023, 2067, 2163, 2120, 2125, 2224 ]) xi1 = np.array([ 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 ]) fi1 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439, 2791, 2576, 2290, 2123 ]) xi = np.array([ 11, 13, 16, 18, 20]) fi = np.array([1503, 7154, 3439, 2576, 2123]) # PROCEDIMIENTO # Tabla de Diferencias Divididas Avanzadas titulo = ['i ','xi ','fi '] n = len(xi) ki = np.arange(0,n,1) tabla = np.concatenate(([ki],[xi],[fi]),axis=0) tabla = np.transpose(tabla) # diferencias divididas vacia dfinita = np.zeros(shape=(n,n),dtype=float) tabla = np.concatenate((tabla,dfinita), axis=1) # Calcula tabla, inicia en columna 3 [n,m] = np.shape(tabla) diagonal = n-1 j = 3 while (j < m): # Añade título para cada columna titulo.append('F['+str(j-2)+']') # cada fila de columna i = 0 paso = j-2 # inicia en 1 while (i < diagonal): denominador = (xi[i+paso]-xi[i]) numerador = tabla[i+1,j-1]-tabla[i,j-1] tabla[i,j] = numerador/denominador i = i+1 diagonal = diagonal - 1 j = j+1 # POLINOMIO con diferencias Divididas # caso: puntos equidistantes en eje x dDividida = tabla[0,3:] n = len(dfinita) # expresión del polinomio con Sympy x = sym.Symbol('x') polinomio = fi[0] for j in range(1,n,1): factor = dDividida[j-1] termino = 1 for k in range(0,j,1): termino = termino*(x-xi[k]) polinomio = polinomio + termino*factor # simplifica multiplicando entre (x-xi) polisimple = polinomio.expand() # polinomio para evaluacion numérica px = sym.lambdify(x,polisimple) # calcula errores en intervalo usado pfi1 = px(xi1) errado1 = np.abs(fi1-pfi1) # Puntos para la gráfica muestras = 101 a = np.min(xi) b = np.max(xi) pxi = np.linspace(a,b,muestras) pfi = px(pxi) # SALIDA np.set_printoptions(precision = 4) print('Tabla Diferencia Dividida') print([titulo]) print(tabla) print('dDividida: ') print(dDividida) print('polinomio: ') print(polinomio) print('polinomio simplificado: ' ) print(polisimple) print('errores en intervalo:') print(xi1) print(errado1) # Gráfica plt.plot(xi0,fi0,'o', label = 'Puntos') plt.plot(xi,fi,'ro', label = 'Puntos') for i in range(0,n,1): etiqueta = '('+str(xi[i])+','+str(fi[i])+')' plt.annotate(etiqueta,(xi[i],fi[i])) plt.plot(pxi,pfi, label = 'Polinomio') plt.legend() plt.xlabel('xi') plt.ylabel('fi') plt.title('Diferencias Divididas - Newton') plt.grid() plt.show()
Literal c y d. Se puede continuar con el algoritmo anterior. Como repaso se adjunta un método diferente al anterior.
# Interpolacion de Lagrange # divisores L solo para mostrar valores import numpy as np import sympy as sym import matplotlib.pyplot as plt # INGRESO , Datos de prueba xi0 = np.array([ 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26 ]) fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344, 4417, 3439, 2791, 2576, 2290, 2123, 2023, 2067, 2163, 2120, 2125, 2224 ]) xi2 = np.array([ 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 ]) fi2 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439, 2791, 2576, 2290, 2123 ]) xi = np.array([ 11, 12, 13, 16, 20]) fi = np.array([1503, 3728, 7154, 3439, 2123]) # PROCEDIMIENTO # Polinomio de Lagrange n = len(xi) x = sym.Symbol('x') polinomio = 0 divisorL = np.zeros(n, dtype = float) for i in range(0,n,1): # Termino de Lagrange numerador = 1 denominador = 1 for j in range(0,n,1): if (j!=i): numerador = numerador*(x-xi[j]) denominador = denominador*(xi[i]-xi[j]) terminoLi = numerador/denominador polinomio = polinomio + terminoLi*fi[i] divisorL[i] = denominador # simplifica el polinomio polisimple = polinomio.expand() # para evaluación numérica px = sym.lambdify(x,polisimple) # calcula errores en intervalo usado pfi2 = px(xi2) errado2 = np.abs(fi2-pfi2) # Puntos para la gráfica muestras = 101 a = np.min(xi) b = np.max(xi) pxi = np.linspace(a,b,muestras) pfi = px(pxi) # SALIDA print(' valores de fi: ',fi) print('divisores en L(i): ',divisorL) print() print('Polinomio de Lagrange, expresiones') print(polinomio) print() print('Polinomio de Lagrange: ') print(polisimple) print('errores en intervalo:') print(xi2) print(errado2) # Gráfica plt.plot(xi0,fi0,'o', label = 'Puntos') plt.plot(pxi,pfi, label = 'Polinomio') plt.legend() plt.xlabel('xi') plt.ylabel('fi') plt.title('Interpolación Lagrange') plt.show()