Categoría: Sol_1Eva 2021-

  • s1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

    Ejercicio: 1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

    literal a

    Se plantea el sistema de ecuaciones de acuerdo a la cantidad de medicamentos que se administra a cada tipo de paciente, siendo las variables X=[a,b,c,d] de acuerdo a la tabla:

    Niños Adolescentes Adultos Adultos Mayores Medicamentos /semana
    Medicamento_A 0.3 0.4 1.1 4.7 3500
    Medicamento_B 1 3.9 0.15 0.25 3450
    Medicamento_C 0 2.1 5.6 1.0 6500
    0.3 a + 0.4 b + 1.1 c + 4.7 d = 3500 a + 3.9 b + 0.15 c + 0.25 d = 3450 0 a +2.1 b + 5.6 c + 1.0 d = 6500

    Mostrando que el número de incógnitas no es igual al número de ecuaciones, por lo que sería necesario de disponer de otra ecuación, o  usar una variable libre.


    literal b

    Siguiendo las indicaciones sobre la variable libre que sea para el grupo de pacientes niños, te tiene que

    0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K 3.9 b + 0.15 c + 0.25 d = 3450 - K 2.1 b + 5.6 c + 1.0 d = 6500 - 0 K

    y haciendo K = 100, se convierte en :

    0.4 b + 1.1 c + 4.7 d = 3500 - 0.3(100) = 3470 3.9 b + 0.15 c + 0.25 d = 3450 - 100 = 3350 2.1 b + 5.6 c + 1.0 d = 6500 - 0 = 6500

    literal c

    Antes de desarrollar el ejercicio para el algoritmo se requiere convertir el sistema de ecuaciones a su forma matricial. Por lo que se plantea la matriz aumentada:

    \begin{pmatrix} 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

    Se realiza el pivoteo parcial por filas, que al aplicar en la primera columa, se intercambia la primera y segunda fila, de tal forma que el mayor valor quede en la primera casilla de la diagonal.

    \begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

    se continua el mismo proceso para la segunda casilla de la diagonal hacia abajo, se aplica también pivoteo parcial,

    \begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \end{pmatrix}

    Con lo que el sistema queda listo para resolver por cualquier método.

    Si seleccionamos Gauss-Seidel, el vector inicial de acuerdo al enunciado será X0=[100,100,100]

    y las ecuaciones a resolver son:

    b = \frac{3350 - 0.15 c - 0.25 d}{3.9} c = \frac{6500 - 2.1 b - 1.0 d}{5.6} d = \frac{3470- 0.4 b - 1.1 c}{4.7}

    iteración 1

    X0=[100,100,100]

    b = \frac{3350 - 0.15(100) - 0.25 (100)}{3.9} = 848.71 c = \frac{6500 - 2.1 (848.71 ) - 1.0 (100)}{5.6} = 824.58 d = \frac{3470- 0.4 (848.71 ) - 1.1 (824.58)}{4.7} = 473.07

    X1=[848.71, 824.58, 473.07]

    diferencia = [848.71, 824.58, 473.07] - [100,100,100]

    diferencia = [748.71, 724.58, 373.07]

    error = max|diferencia| = 748.71

    resultado del error es mayor que tolerancia de 0.01, se continua con la siguiente iteración.

    iteración 2

    X1=[848.71, 824.58, 473.07]

    b = \frac{3350 - 0.15 (824.58) - 0.25 (473.07)}{3.9} = 796.93 c = \frac{6500 - 2.1 (796.93) - 1.0 (473.07)}{5.6} = 777.38 d = \frac{3470- 0.4 (796.93) - 1.1 (777.38)}{4.7} = 488.53

    X2 = [796.93, 777.38, 488.53]

    diferencia = [796.93, 777.38, 488.53]  - [848.71, 824.58, 473.07]

    diferencia = [-51.78 ,  -47.2, 15.52]

    error = max|diferencia| = 51.78

    iteración 3

    X2 = [796.93, 777.38, 488.53]

    b = \frac{3350 - 0.15 (777.38) - 0.25 (488.53)}{3.9} = 797.75 c = \frac{6500 - 2.1 (797.75) - 1.0 (488.53)}{5.6} = 774.31 d = \frac{3470- 0.4 (797.75) - 1.1 (774.31)}{4.7} = 489.18

    x3 = [ 797.75, 774.31, 489.18]

    diferencias = [ 797.75, 774.31, 489.18] - [796.93, 777.38, 488.53]

    diferencias = [0.82, -3.07, 0.65]

    error = max|diferencia| = 3.07

    Observación, el error disminuye en cada iteración, por lo que el sistema converge.

    solución con el algoritmo, solo se toma la parte entera de la respuesta, pues los pacientes son números enteros.

    respuesta X: [797.83, 774.16, 489.20]
    

    con lo que el primer resultado es:

     X = [797, 774, 489]

    literal d

    Si la cantidad de pacientes presentados en una seman es la que se indica en el enunciado [350,1400,1500,1040], se determina que el sistema hospitalario estatal no podrá atender a todos los pacientes al compara con la capacidad encontrada en el literal c.

    Hay un exceso de 2129 pacientes encontrados por grupo:

    exceso = [350, 1400, 1500, 1040] - [100, 797, 774, 489] = [250, 603, 726, 551]

    con lo que se recomienda una medida de confinamiento para disminuir los contagios y poder atender a los pacientes que se presenten.

    literal e

    Siendo K = 0, al estar vacunados todos los niños, y suponiendo que la vacuna es una cura, se tiene:

    0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K = 3500 3.9 b + 0.15 c + 0.25 d = 3450 - K = 3450 2.1 b + 5.6 c + 1.0 d = 6500 - 0 k = 6500

    pivoteado parcialmente por filas se tiene:

    \begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3450 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3500 \end{pmatrix}

    Resolviendo por un método directo:

    se obtiene:

    respuesta X: [823.46, 763.35, 495.94]

    que al compararse con la capacidad anterior en números enteros se encuentra una diferencia de un incremento de 21 pacientes neto ante la condición de usar todos los medicamentos disponibles.

    diferencia = [823, 763, 495] - [797, 774, 489] = [26, -11, 6]


     Instrucciones en Python

    # Método de Gauss-Seidel
    # solución de sistemas de ecuaciones
    # por métodos iterativos
    
    import numpy as np
    
    # INGRESO
    A = np.array([[0.4, 1.10, 4.7 ],
                  [3.9, 0.15, 0.25],
                  [2.1, 5.60, 1.0  ]])
    k = 100
    B = np.array([3500.-0.3*k,3450-1*k,6500-0*k])
    
    X0  = np.array([100.0,100,100])
    
    tolera = 0.001
    iteramax = 100
    
    # PROCEDIMIENTO
    # Matriz aumentada
    Bcolumna = np.transpose([B])
    AB  = np.concatenate((A,Bcolumna),axis=1)
    AB0 = np.copy(AB)
    
    # Pivoteo parcial por filas
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no está en diagonal
        if (dondemax !=0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:]  = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal
    
    A = np.copy(AB[:,:n])
    B = np.copy(AB[:,n])
    
    # Gauss-Seidel
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    #  valores iniciales
    X = np.copy(X0)
    diferencia = np.ones(n, dtype=float)
    errado = 2*tolera
    
    itera = 0
    while not(errado<=tolera or itera>iteramax):
        # por fila
        for i in range(0,n,1):
            # por columna
            suma = 0 
            for j in range(0,m,1):
                # excepto diagonal de A
                if (i!=j): 
                    suma = suma-A[i,j]*X[j]
            
            nuevo = (B[i]+suma)/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        print(X)
        errado = np.max(diferencia)
        itera = itera + 1
    
    # Respuesta X en columna
    X = np.transpose([X])
    
    # revisa si NO converge
    if (itera>iteramax):
        X=0
    # revisa respuesta
    verifica = np.dot(A,X)
    
    # SALIDA
    print('respuesta X: ')
    print(X)
    print('verificar A.X=B: ')
    print(verifica)
    
  • s1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

    Ejercicio: 1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

    Literal a

    Evaluando las sucesión de la forma recursiva:

    xi
    [-0.45       -0.4383     -0.4458     -0.441      -0.4441 
     -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
     -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
    errores
    [ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
     -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
      1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]
    

    literal b

    Se puede afirmar que converge, observe la diferencia entre cada dos valores consecutivos de la sucesión ... (continuar de ser necesario)

    literal c

    Para el algoritmo se requiere la función f(x) y su derivada f'(x)

    f(x) = x +ln(x+2) f'(x) = 1 + \frac{1}{x+2}

    x0 = -0.45

    itera = 1

    x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} x_{1} = x_0 -\frac{f(x_0)}{f'(x_0)} = -0.45 -\frac{-0.45+ln(-0.45+2)}{1 + \frac{1}{-0.45+2}}

    x1 = -0.44286

    error = |x1-x0| = |-0.44286 -(-0.45)| = 0.007139

    itera = 2

    x_{2} = -0.4428 -\frac{-0.4428 + ln(-0.45+2)}{1 + \frac{1}{-0.4428+2}}

    x1 = -0.44286

    error = |x1-x0| = |-0.44285 -(-4.4286)| = 6.4394e-06

    con lo que se cumple el valor de tolerancia y no se requiere otra iteración

    la raiz se encuentra en x = -0.44286

    Solución con algoritmo

    xi
    [-0.45       -0.4383     -0.4458     -0.441      -0.4441 
     -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
     -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
    errores
    [ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
     -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
      1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]
    ['xi', 'xnuevo', 'tramo']
    [[-4.5000e-01 -4.4286e-01  7.1392e-03]
     [-4.4286e-01 -4.4285e-01  6.4394e-06]]
    raiz en:  -0.44285440100759543
    con error de:  6.439362322474551e-06
    >>> 
    

    Instrucciones en Python

    import numpy as np
    import matplotlib.pyplot as plt
    
    # literal a sucesión en forma recursiva
    def secuenciaL(n,x0):
        if n == 0:
            xn = x0
        if n>0:
            xn = np.log(1/(2+secuenciaL(n-1,x0)))
        return(xn)
    x0 = -0.45
    n = 15
    xi = np.zeros(n,dtype=float)
    xi[0] = x0
    errado = np.zeros(n,dtype=float)
    for i in range(1,n,1):
        xi[i] = secuenciaL(i,x0)
        errado[i-1] = xi[i] - xi[i-1]
       
    np.set_printoptions(precision=4)
    print('xi: ')
    print(xi)
    print('errado: ')
    print(errado)
    
    #Grafica literal a y b
    plt.plot(xi)
    plt.plot(xi,'o')
    plt.xlabel('n')
    plt.ylabel('xn')
    plt.show()
    
    # Método de Newton-Raphson
    # Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
    
    import numpy as np
    
    # INGRESO
    fx  = lambda x: x+np.log(x+2)
    dfx = lambda x: 1+ 1/(x+2)
    
    x0 = -0.45
    tolera = 0.0001
    
    a = -0.5
    b = -0.2
    muestras = 21
    # PROCEDIMIENTO
    tabla = []
    tramo = abs(2*tolera)
    xi = x0
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo  = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,tramo])
        xi = xnuevo
    
    # convierte la lista a un arreglo.
    tabla = np.array(tabla)
    n = len(tabla)
    
    # para la gráfica
    xj = np.linspace(a,b,muestras)
    fj = fx(xj)
    
    # SALIDA
    print(['xi', 'xnuevo', 'tramo'])
    np.set_printoptions(precision = 4)
    print(tabla)
    print('raiz en: ', xi)
    print('con error de: ',tramo)
    
    plt.plot(xj,fj)
    plt.axhline(0, color='grey')
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.show()
    

    litera d: tarea

  • s1Eva_2021PAOI_T3 Interpolar, modelo de contagios 2020

    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^4

    literal 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()