4.9 Interpolar – Mascota por 50 años

Ejercicio de interpolación

  • Desarrollar por grupos
  • integrar soluciones en un mismo gráfico
  • mostrar polinomios y rangos de existencia

# caparazon superior
xi = [113., 117, 134, 153, 169, 184, 194, 203]
fi = [127., 141, 161, 166, 160, 155, 140, 132]
puntos = [[xi,fi]]

# caparazon inferior 1
xi = [113., 123, 130, 149, 182, 197, 208, 211]
fi = [127., 116, 112, 107, 110, 114, 112, 108]
puntos.append([xi,fi])

# caparazon inferior 2
xi = [107., 114, 120, 143, 170, 192, 210]
fi = [124., 114, 108, 99, 99, 106, 107]
puntos.append([xi,fi])

# Patas 01
xi = [110., 116, 120, 121, 134, 143, 153, 168, 173, 177, 181, 188, 194, 201, 205, 210, 214, 218]
fi = [92., 86, 89, 86, 87, 95, 91, 93, 85, 85, 83, 86, 88, 86, 91, 87, 92, 92]
puntos.append([xi,fi])

# Patas 02
xi = [109., 115, 117, 120, 125, 130, 134, 138, 142, 145, 150, 152, 154]
fi = [91., 87, 82, 82, 78, 78, 81, 78, 76, 79, 77, 80, 86]
puntos.append([xi,fi])

# Patas 03
xi = [172., 175, 182, 186, 191, 195, 201, 205, 210, 211, 217]
fi = [86., 79, 79, 77, 82, 78, 81, 79, 82, 85, 90]
puntos.append([xi,fi])

# Patas 04
xi = [113., 118, 122, 127, 132, 136, 140, 144, 149, 152]
fi = [88., 90, 88, 86, 85, 83, 87, 84, 86, 85]
puntos.append([xi,fi])

# Rabito 01
xi = [97., 102, 108, 111, 117, 121]
fi = [120., 113, 108, 105, 102, 102]
puntos.append([xi,fi])

# cabeza 01
xi = [194., 196, 203, 210, 211, 209]
fi = [177., 167, 157, 149, 138, 135]
puntos.append([xi,fi])

# cabeza 02
xi = [195., 199, 207, 214, 219, 224, 229, 234, 239, 242, 244]
fi = [177., 185, 190, 190, 188, 193, 192, 185, 175, 172, 168]
puntos.append([xi,fi])

# cabeza 03
xi = [220., 226, 234, 239, 241]
fi = [164., 162, 163, 163, 163]
puntos.append([xi,fi])

# cabeza 04
xi = [203., 211, 214, 219, 223, 224, 225, 230, 236, 241]
fi = [115., 119, 125, 124, 127, 137, 146, 149, 154, 162]
puntos.append([xi,fi])

# cabeza 05
xi = [208., 212, 215, 219, 221, 225, 228]
fi = [174., 177, 178, 179, 182, 183, 178]
puntos.append([xi,fi])

# cabeza 06
xi = [206., 210, 214, 218, 220, 224, 229, 233, 232]
fi = [179., 182, 182, 181, 174, 171, 170, 175, 180]
puntos.append([xi,fi])

4.8 Interpolar – Mascota descansando

Referencia: Burden 9th Ed. Ejercicio 32 p164

La parte superior de una mascota descansando se describe con los puntos presentados en la tabla.  Los puntos se usan para ejercicios de interpolación.

xi = [1,2,5,6,7,8,10,13,17,20,23,24,25,27,27.7,28,29,30]
fi = [3.0,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5,7.0,6.1,5.6,5.8,5.2,4.1,4.3,4.1,3.0]

4.7 Interpolar – Pato en pleno vuelo

Referencia: Burden Cap 3. Ejemplo 1.

La figura muestra un joven pato en pleno vuelo. Para aproximar el perfil de la parte superior del pato, se presentan 21 puntos a lo largo de la curva de aproximación relativos a un sistema de coordenadas sobrepuestas.

xi = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fi = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]


a) Realice la interpolación entre puntos usando la interpolación de Lagrange para la parte superior del pato en vuelo.

b) ¿Es posible crear un solo polinomio para los 21 puntos presentados? Explique su respuesta en pocas palabras.

c) ¿En concepto, cuando considera necesario crear un nuevo polinomio?

d) Adjunte en los archivos para: la gráfica que interpola los puntos, y los polinomios que la generan.

e) (Puntos extra: 10, para promediar con lección o taller)

Se adjuntan los demás puntos del perfil del pato en pleno vuelo, presente el o los polinomios que interpolan la figura (gráfica y polinomios), las despuestas a literales b,c,d deben también ser válidas para el ejercicio completo.

# Perfil Superior del pato
xiA = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fiA = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]

# Perfil inferior cabeza
xiB = [0.817, 0.897, 1.022, 1.191, 1.510, 1.834, 2.264, 2.962, 3.624, 4.202, 4.499, 4.779, 5.109, 5.527]
fiB = [1.180, 1.065, 1.023, 1.010, 1.032, 1.085, 1.192, 1.115, 1.087, 1.100, 0.830, 0.608, 0.350, 0.106]

# Perfil Ala superior
xiC = [4.659, 4.865, 5.085, 5.261, 5.387, 5.478, 5.527]
fiC = [-5.161, -4.741, -3.933, -2.951, -1.970, -0.981, 0.106]

# Perfil Ala inferior
xiD = [4.659, 4.750, 4.990, 5.289, 5.560, 5.839, 6.113, 6.606, 6.916, 7.305, 7.563, 7.802, 7.983]
fiD = [-5.161, -5.259, -5.284, -5.268, -5.161, -4.982, -4.769, -4.286, -3.911, -3.213, -2.670, -2.176, -1.655]

# Perfil inferior posterior
xiE = [8.141, 8.473, 8.832, 9.337, 9.887, 10.572, 10.995, 11.501, 11.923, 12.364, 12.763, 13.300]
fiE = [-1.138, -0.434, -0.514, -0.494, -0.382, -0.005, -0.090, -0.085, -0.030, 0.093, 0.120, 0.250]

# Perfil ojo superior
xiF = [2.663, 2.700, 2.805, 2.886]
fiF = [2.202, 2.279, 2.293, 2.222]

# Perfil ojo inferior
xiG = [2.663, 2.720, 2.826, 2.886]
fiG = [2.202, 2.130, 2.143, 2.222]

4.6 Métodos de interpolación con gráficos animados en Python

[ Dif_Finitas ] [ Dif_Finitas avanzadas ] [ Dif_Divididas_Newton ]

Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para Interpolación polinómica, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

En los algoritmos, se convierten las partes en funciones, que se usan para generar los polinomios para cada grado y se incorporan en un gráfico animado con los resultados presentados.

Para la gráfica animada se añaden las instrucciones siguientes:

Se determina el intervalo para el eje x usando los valores mínimos y máximos del vector xi. Se toma como muestras para el polinomio las suficientes para una línea suave, que pueden ser mayores a 21 y se crean los valores para el polinomio en p_xi.

Para cada polinomio se guarda en una tabla los valores px(p_xi)del polinomio evaluado el vector p_xi

La gráfica (graf_ani) se crea en una ventana (fig_ani), inicializando con los puntos [xi,fi] en color rojo, y configurando los parámetros base para el gráfico.

Se usan procedimientos para crear unatrama() para cada polinomio y en cada cambio se limpia la trama manteniendo la base con limpiatrama().

En caso de requerir un archivo gif animado se proporciona un nombre de archivo. Para crear el archivo se requiere de un complemento ‘imagemagick‘ a ser instalado.

otros ejemplos de animación en el curso de Fundamentos de Programación:

Movimiento circular – Una partícula, animación con matplotlib-Python


Interpolación por Diferencias Finitas Avanzadas

Diferencias Finitas Avanzadas GIFanimado

En la función para interpolación se añade la verificación de tamaños de paso equidistantes o iguales. En el caso de no tener tamaños de paso equidistantes

tabla de diferencias finitas
['i', 'xi', 'fi', 'df1', 'df2']
[[0.   0.1  1.45 0.15 0.  ]
 [1.   0.2  1.6  0.   0.  ]]
dfinita:  [0.15 0.  ]
1.45 +
+( 0.15 / 0.1 )* ( x - 0.1 )
polinomio simplificado
1.5*x + 1.3

tabla de diferencias finitas
['i', 'xi', 'fi', 'df1', 'df2', 'df3']
[[ 0.    0.1   1.45  0.15 -0.05  0.  ]
 [ 1.    0.2   1.6   0.1   0.    0.  ]
 [ 2.    0.3   1.7   0.    0.    0.  ]]
dfinita:  [ 0.15 -0.05  0.  ]
1.45 +
+( 0.15 / 0.1 )* ( x - 0.1 )
+( -0.05 / 0.02 )*  (x - 0.2)*(x - 0.1) 
polinomio simplificado
-2.5*x**2 + 2.25*x + 1.25

tabla de diferencias finitas
['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]
dfinita:  [ 0.15 -0.05  0.25  0.  ]
1.45 +
+( 0.15 / 0.1 )* ( x - 0.1 )
+( -0.05 / 0.02 )*  (x - 0.2)*(x - 0.1) 
+( 0.25 / 0.006 )*  (x - 0.3)*(x - 0.2)*(x - 0.1) 
polinomio simplificado
41.66667*x**3 - 27.500002*x**2 + 6.8333337*x + 0.99999998

Polinomios con Diferencias Finitas Avanzadas
px_0 =
1.45

px_1 =
1.5*x + 1.3

px_2 =
       2                
- 2.5*x  + 2.25*x + 1.25

px_3 =
          3              2                           
41.66667*x  - 27.500002*x  + 6.8333337*x + 0.99999998

El resumen de las instrucciones se presentan a continuación.

# Interpolación -Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym

def dif_finitas(xi,fi, precision=5, vertabla=False):
    '''Genera la tabla de diferencias finitas
    resultado en: tabla,titulo
    Tarea: verificar tamaño de vectores
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    # Tabla de Diferencias Finitas
    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 finitas 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('df'+str(j-2))
        # cada fila de columna
        i = 0
        while (i < diagonal):
            tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
            i = i+1
        diagonal = diagonal - 1
        j = j+1
    if vertabla==True:
        np.set_printoptions(precision)
        print('tabla de diferencias finitas')
        print(titulo)
        print(tabla)
    return(tabla, titulo)

def pasosEquidistantes(xi, casicero = 1e-15):
    ''' Revisa tamaños de paso h en vector xi.
    True:  h son equidistantes,
    False: h tiene tamaño de paso diferentes y dónde.
    '''
    xi = np.array(xi,dtype=float)
    n = len(xi)
    # revisa tamaños de paso equidistantes
    h_iguales = True
    if n>3: 
        dx = np.zeros(n,dtype=float)
        for i in range(0,n-1,1): # calcula hi como dx
            dx[i] = xi[i+1]-xi[i]
        for i in range(0,n-2,1): # revisa diferencias
            dx[i] = dx[i+1]-dx[i]
            if dx[i]<=casicero: # redondea cero
                dx[i]=0
            if abs(dx[i])>0:
                h_iguales=False
                print('tamaños de paso diferentes en i:',i+1,',',i+2)
        dx[n-2]=0
    return(h_iguales)

def interpola_dfinitasAvz(xi,fi, vertabla=False,
                          precision=5, casicero = 1e-15):
    '''Interpolación de diferencias finitas
    resultado: polinomio en forma simbólica,
    redondear a cero si es menor que casicero 
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    n = len(xi)
    # revisa tamaños de paso equidistantes
    h_iguales = pasosEquidistantes(xi, casicero)
    if vertabla==True:
        np.set_printoptions(precision)
    # POLINOMIO con diferencias Finitas avanzadas
    x = sym.Symbol('x')
    polisimple = sym.S.NaN # expresión del polinomio con Sympy
    if h_iguales==True:
        tabla,titulo = dif_finitas(xi,fi,vertabla)
        h = xi[1] - xi[0]
        dfinita = tabla[0,3:]
        if vertabla==True:
            print('dfinita: ',dfinita)
            print(fi[0],'+')
        n = len(dfinita)
        polinomio = fi[0]
        for j in range(1,n,1):
            denominador = np.math.factorial(j)*(h**j)
            factor = np.around(dfinita[j-1]/denominador,precision)
            termino = 1
            for k in range(0,j,1):
                termino = termino*(x-xi[k])
            if vertabla==True:
                txt1='';txt2=''
                if n<=2 or j<=1:
                    txt1 = '('; txt2 = ')'
                print('+(',np.around(dfinita[j-1],precision),
                      '/',np.around(denominador,precision),
                      ')*',txt1,termino,txt2)
            polinomio = polinomio + termino*factor
        # simplifica multiplicando entre (x-xi)
        polisimple = polinomio.expand() 
    if vertabla==True:
        print('polinomio simplificado')
        print(polisimple)
    return(polisimple)

# INGRESO , Datos de prueba
xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
# tabla polinomios
n = len(xi)
px_tabla = [fi[0]]
for grado in range(1,n,1):
    polinomio = interpola_dfinitasAvz(xi[:grado+1],fi[:grado+1],
                                      vertabla=True, precision=4)
    print('',)
    px_tabla.append(polinomio)

# SALIDA
print('Polinomios con Diferencias Finitas Avanzadas')
for grado in range(0,n,1):
    print('px_'+str(grado)+' =') #, px_tabla[grado])
    sym.pprint(px_tabla[grado])
    print()


Parte adicional para la gráfica con GIF animado

# GRAFICA CON ANIMACION polinomios --------
import matplotlib.pyplot as plt
import matplotlib.animation as animation

unmetodo = 'Diferencias Finitas Avanzadas'
narchivo = 'DifFinitasAvz' # nombre archivo
muestras = 51 # de cada p(X)

# Puntos para la gráfica
a = np.min(xi)
b = np.max(xi)
p_xi = np.linspace(a,b,muestras)

# lineas por grado de polinomio
x = sym.Symbol('x')
px_lineas = np.zeros(shape=(n,muestras), dtype=float)
for grado in range(0,n,1):
    polinomio = px_tabla[grado]
    px = sym.utilities.lambdify(x,polinomio,'numpy')
    px_lineas[grado] = px(p_xi)

# Parametros de trama/foto
retardo = 800   # milisegundos entre tramas
tramas = len(px_lineas)
ymax = np.max(fi)
ymin = np.min(fi)
deltax = 0.1*np.abs(b-a)
deltay = 0.1*np.abs(ymax-ymin)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
# Función Base
fx_linea, = graf_ani.plot(xi,fi,'o',color='red')
# Polinomios de px_tabla grado = 0
px_unalinea, = graf_ani.plot(p_xi, px_lineas[0],
                         label='grado: 0')
# Configura gráfica
graf_ani.set_xlim([a-deltax,b+deltax])
graf_ani.set_ylim([ymin-deltay,ymax+deltay])
graf_ani.axhline(0, color='k')  # Linea horizontal en cero
graf_ani.set_title('Polinomio - '+unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('p(x)')
graf_ani.grid()

# Cuadros de texto en gráfico
txt_x = (b+a)/2
txt_y = ymax
txt_poli = graf_ani.text(txt_x, txt_y,'p(x):',
                      horizontalalignment='center')
txt_grado = graf_ani.text(txt_x, txt_y-deltay,'grado:',
                        horizontalalignment='center')
# Nueva Trama
def unatrama(i,p_xi,pxi): 
    # actualiza cada linea
    px_unalinea.set_xdata(p_xi)
    px_unalinea.set_ydata(pxi[i])
    unpolinomio = px_tabla[i]
    if unpolinomio == sym.S.NaN:
        unpolinomio = 'h pasos no equidistantes' 
    etiquetap = 'p'+str(i)+'(x) = '+str(unpolinomio)
    px_unalinea.set_label(etiquetap)
    # actualiza texto
    txt_poli.set_text(etiquetap)
    txt_grado.set_text('Grado: '+str(i))
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numerocolor = i%10
        lineacolor = 'C'+str(numerocolor)
    px_unalinea.set_color(lineacolor)
    
    return (px_unalinea, txt_poli, txt_grado)
# Limpia Trama anterior
def limpiatrama(): 
    px_unalinea.set_ydata(np.ma.array(p_xi, mask=True))
    px_unalinea.set_label('')
    txt_poli.set_text('')
    txt_grado.set_text('')
    return (px_unalinea,txt_poli, txt_grado)

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (p_xi,px_lineas),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

[ Dif_Finitas ] [ Dif_Finitas avanzadas ] [ Dif_Divididas_Newton ]


Interpolación por Diferencias Divididas de Newton

Diferencias Divididas Newton GIF animado

tabla de diferencias divididas
['i', 'xi', 'fi', 'F[1]', 'F[2]']
[[0.   0.1  1.45 1.5  0.  ]
 [1.   0.2  1.6  0.   0.  ]]
difDividida:  [1.5 0. ]
1.45 +
+( 1.5 )*( x - 0.1 )
polinomio simplificado
1.5*x + 1.3

tabla de diferencias divididas
['i', 'xi', 'fi', 'F[1]', 'F[2]', 'F[3]']
[[ 0.    0.1   1.45  1.5  -2.5   0.  ]
 [ 1.    0.2   1.6   1.    0.    0.  ]
 [ 2.    0.3   1.7   0.    0.    0.  ]]
difDividida:  [ 1.5 -2.5  0. ]
1.45 +
+( 1.5 )*( x - 0.1 )
+( -2.5 )* (x - 0.2)*(x - 0.1) 
polinomio simplificado
-2.5*x**2 + 2.25*x + 1.25

tabla de diferencias divididas
['i', 'xi', 'fi', 'F[1]', 'F[2]', 'F[3]', 'F[4]']
[[ 0.00000e+00  1.00000e-01  1.45000e+00  1.50000e+00 -2.50000e+00
   5.00000e+00  0.00000e+00]
 [ 1.00000e+00  2.00000e-01  1.60000e+00  1.00000e+00  3.33067e-15
   0.00000e+00  0.00000e+00]
 [ 2.00000e+00  3.00000e-01  1.70000e+00  1.00000e+00  0.00000e+00
   0.00000e+00  0.00000e+00]
 [ 3.00000e+00  6.00000e-01  2.00000e+00  0.00000e+00  0.00000e+00
   0.00000e+00  0.00000e+00]]
difDividida:  [ 1.5 -2.5  5.   0. ]
1.45 +
+( 1.5 )*( x - 0.1 )
+( -2.5 )* (x - 0.2)*(x - 0.1) 
+( 5.0 )* (x - 0.3)*(x - 0.2)*(x - 0.1) 
polinomio simplificado
5.0*x**3 - 5.5*x**2 + 2.8*x + 1.22

Polinomios con Diferencias Divididas Newton
px_0 =
1.45

px_1 =
1.5*x + 1.3

px_2 =
       2                
- 2.5*x  + 2.25*x + 1.25

px_3 =
     3        2               
5.0*x  - 5.5*x  + 2.8*x + 1.22
>>> 

Instrucciones en Python

# Interpolación -Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
import numpy as np
import sympy as sym

def dif_divididas(xi,fi, vertabla=False,
                  precision=5, casicero = 1e-15):
    '''Genera la tabla de diferencias divididas
    resultado en: tabla, titulo
    Tarea: verificar tamaño de vectores
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    # Tabla de Diferencias Divididas
    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
            if np.abs(tabla[i,j])<= casicero:
                tabla[i,j] = 0
            i = i+1
        diagonal = diagonal - 1
        j = j+1

    if vertabla==True:
        np.set_printoptions(precision)
        print('tabla de diferencias divididas')
        print(titulo)
        print(tabla)
    return(tabla, titulo)

def interpola_difDividida(xi,fi, vertabla=False,
                          precision=5, casicero = 1e-15):
    '''Interpolación de diferencias finitas
    resultado: polinomio en forma simbólica,
    redondear a cero si es menor que casicero 
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)

    # Tabla de Diferencias Divididas
    tabla,titulo = dif_divididas(xi,fi,vertabla,
                                 precision,casicero)
    dDividida = tabla[0,3:]
    n = len(dDividida)
    x = sym.Symbol('x')
    polinomio = fi[0]
    if vertabla==True:
        print('difDividida: ',dDividida)
        print(fi[0],'+')
    for j in range(1,n,1):
        factor = np.around(dDividida[j-1],precision)
        termino = 1
        for k in range(0,j,1):
            termino = termino*(x-xi[k])
        if vertabla==True:
            txt1='';txt2=''
            if n<=2 or j<=1:
                txt1 = '('; txt2 = ')'
            print('+(',factor,')*'+txt1,termino,txt2)
        polinomio = polinomio + termino*factor
    
    # simplifica multiplicando entre (x-xi)
    polisimple = polinomio.expand() 
    if vertabla==True:
        print('polinomio simplificado')
        print(polisimple)
    return(polisimple)

# INGRESO , Datos de prueba
xi = [0.10, 0.2, 0.3, 0.6]
fi = [1.45, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
# tabla polinomios
n = len(xi)
px_tabla = [fi[0]]
for grado in range(1,n,1):
    polinomio = interpola_difDividida(xi[:grado+1],fi[:grado+1],
                                          vertabla=True)
    print('',)
    px_tabla.append(polinomio)

# SALIDA
print('Polinomios con Diferencias Divididas Newton')
for grado in range(0,n,1):
    print('px_'+str(grado)+' =') #, px_tabla[grado])
    sym.pprint(px_tabla[grado])
    print()

Para la gráfica animada se usa el mismo bloque de instrucciones del método de Diferencias Finitas avanzadas, solo requiere cambiar el nombre del método y el nombre para el archivo GIF animado.


Interpolación por el Método de Lagrange

Interpola con Lagrange
 1.45 * (x - 0.4)*(x - 0.3)*(x - 0.2) / (-0.2 + 0.1)*(-0.3 + 0.1)*(-0.4 + 0.1) 
+ 1.6 * (x - 0.4)*(x - 0.3)*(x - 0.1) / (-0.1 + 0.2)*(-0.3 + 0.2)*(-0.4 + 0.2) 
+ 1.7 * (x - 0.4)*(x - 0.2)*(x - 0.1) / (-0.1 + 0.3)*(-0.2 + 0.3)*(-0.4 + 0.3) 
+ 2.0 * (x - 0.3)*(x - 0.2)*(x - 0.1) / (-0.1 + 0.4)*(-0.2 + 0.4)*(-0.3 + 0.4) 
polinomio simplificado
41.666667*x**3 - 27.5*x**2 + 6.833333*x + 1.0
Interpolación con Lagrange
41.666667*x**3 - 27.5*x**2 + 6.833333*x + 1.0
>>> 

Instrucciones en Python

# Interpolación - Lagrange
# Tarea: Verificar tamaño de vectores,
import numpy as np
import sympy as sym

def interpola_Lagrange(xi,fi,vertabla=False,
                       precision=6, casicero = 1e-15):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio con Lagrange
    if vertabla==True:
        print('Interpola con Lagrange')
    polinomio = sym.S.Zero
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        numerador = 1
        denominador = 1
        for j  in range(0,n,1):
            if (j!=i):
                numerador = numerador*(x-xi[j])
                denominador = denominador*(sym.UnevaluatedExpr(xi[i])-xi[j])
        if vertabla==True:
            txt0='' ; txt1='('; txt2=')'
            if i>0:
                txt0='+'
            if n>2:
                txt1=''; txt2=''
            print(txt0,fi[i],'*'+txt1,numerador,txt2+'/'+ txt1,
                  denominador,txt2)
        #factor = np.around(fi[i]/float(denominador.doit()),precision) 
        polinomio = polinomio + (fi[i]/denominador.doit())*numerador
    # Expande el polinomio
    polisimple = polinomio.expand()
    polisimple = redondea_coef(polisimple, precision)
    if vertabla==True:
        print('polinomio simplificado')
        print(polisimple)
    return(polisimple)

def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
    ''' redondea coeficientes de términos suma de una ecuacion
    '''
    tipo = type(ecuacion)
    tipo_eq = False
    if tipo == sym.core.relational.Equality:
        RHS = ecuacion.rhs
        ecuacion = ecuacion.lhs
        tipo = type(ecuacion)
        tipo_eq = True

    if tipo == sym.core.add.Add: # términos suma de ecuacion
        term_sum = sym.Add.make_args(ecuacion)
        ecuacion = sym.S.Zero
        for term_k in term_sum:
            # factor multiplicativo de termino suma
            term_mul = sym.Mul.make_args(term_k)
            producto = sym.S.One
            for factor in term_mul:
                if not(factor.has(sym.Symbol)):
                    factor = np.around(float(factor),precision)
                    if (abs(factor)%1)<casicero: # si es entero
                        factor = int(factor)
                producto = producto*factor
            ecuacion = ecuacion + producto
    if tipo == sym.core.mul.Mul: # termino único, busca factores
        term_mul = sym.Mul.make_args(ecuacion)
        producto = sym.S.One
        for factor in term_mul:
            if not(factor.has(sym.Symbol)):
                factor = np.around(float(factor),precision)
                print(factor)
                if (abs(factor)%1)<casicero: # si es entero
                    factor = int(factor)
            producto = producto*factor
        ecuacion = producto
    if tipo == float: # si es entero
        if (abs(ecuacion)%1)<casicero: 
            ecuacion = int(ecuacion)
    if tipo_eq:
        ecuacion = sym.Eq(ecuacion,RHS)
    return(ecuacion)

# INGRESO , Datos de prueba
xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
# tabla polinomios
polisimple = interpola_Lagrange(xi,fi,
                                   vertabla=True)

# SALIDA
print('Interpolación con Lagrange')
print(polisimple)

Para la gráfica animada se usa el mismo bloque de instrucciones del método de Diferencias Finitas avanzadas, solo requiere cambiar el nombre del método y el nombre para el archivo GIF animado.


 

4.5 Interpolación paramétrica con Python

Referencia: Rodriguez 6.9.2 p236, Burden 9Ed 3.6 p164

En algunos casos, los datos (x,y) no tienen una relación de tipo funcional y(x), entonces no se pueden aplicar directamente los métodos de interpolación revisados.

Por ejemplo, en la trayectoria del balón en el «gol imposible», la gráfica de la trayectoria en el espacio o sus proyecciones en los planos dependen del parámetro tiempo «t»en lugar de una relación de x,y,z

Referencia: 1Eva_IT2018_T4 El gol imposible

Tabla de datos:

ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]

Sin embargo si las coordenadas (x,y) se expresan como funciones de otra variable t denominada parámetro, entonces los puntos x(t), y(t) tienen relación funcional, y se pueden construir polinomios de interpolación.

Solución propuesta: s1Eva_IT2018_T4 El gol imposible


Ejemplo

Las coordinadas x(t) y y(t) del recorrido de un cohete registradas en los instantes t fueron:

ti  = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

Usaremos un algoritmo en Python para mostrar la trayectoria x,y para el problema planteado.

Al realizar la interpolación de los puntos para obtener polinomios que dependen de «t» se obtiene:

px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

polinomios con los que se puede realizar la gráfica px(t), py(t) en forma separada. Pero para comprender mejor la trayectoria del cohete, se utiliza la gráfica px(t) vs py(t) en el intervalo t entre[0,3]

Las intrucciones para mostrar el resultado son:

# interpolación paramétrica
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
ti = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

# PROCEDIMIENTO
# interpolando con lagrange
px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

t = np.arange(0,3,0.01)
puntosx = px(t)
puntosy = py(t)

# Salida
plt.plot(puntosx,puntosy)
plt.show()

4.4 Trazadores Cúbicos Natural o libre con Python

Referencia:  Burden 3.5 p105, Chapra 18.6.3 p532, Rodríguez 6.11.1 p244

Tiene como objetivo incorporar condiciones adicionales para la función en los extremos del intervalo donde se encuentran los puntos o «nodos».

Por ejemplo, si los datos son los de una trayectoria en un experimento de física, se podría disponer de la aceleración el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.

El método busca obtener un polinomio de tercer grado para cada sub-intervalo o tramo entre «nodos» consecutivos (j,  j+1), de la forma:

S_j(x_j) = a_j + b_j (x-x_j) + c_j(x-xj)^2 + d_j(x-x_j)^3

para cada j = 0, 1, …, n-1

Para n datos existen n-1 tramos, cuatro incógnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1) incógnitas para todo el intervalo .

Para los términos (xj+1xj) de los tramos que se usan varias veces en el desarrollo, se simplifican como hj:

h_j = x_{j+1} - x_{j} S_j(x_j) = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3

Para generar el sistema de ecuaciones, se siguen los siguientes criteros:

1.  Los valores de la función deben ser iguales en los nodos interiores

S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1})

2. Las primeras derivadas en los nodos interiores deben ser iguales

S'_j(x_{j+1}) = S'_{j+1}(x_{j+1})

3. Las segundas derivadas en los nodos interiores deben ser iguales

S''_j(x_{j+1}) = S''_{j+1}(x_{j+1})

4. El primer y último polinomio deben pasar a través de los puntos extremos

f(x_0) = S_0(x_0) = a_0 f(x_n) = S_n(x_n) = a_n

5. Una de las condiciones de frontera se satisface:

5a. frontera libre o natural: Las segundas derivadas en los nodos extremos son cero

S''(x_0) = S''(x_n) = 0

5b. frontera sujeta: las primeras derivadas en los nodos extremos son conocidas

S'(x_0) = f'(x_0)

S'(x_n) = f'(x_n)

Tarea: Revisar en los textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.

Ubicar en el texto las ecuaciones resultantes, que son las que se aplicarán en el algoritmo.


Algoritmo en Python

El algoritmo parte de lo desarrollado para «trazadores lineales», donde se presentaron los bloques para:

  • construir el trazador en una tabla de polinomios por tramos
  • graficar los trazadores en cada tramo,
  • evaluar cada uno de ellos en cada tramo con muestras suficientes para presentar una buena precisión en la gráfica

Del algoritmo básico se modifica entonces el bloque del cálculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores cúbicos naturales (enviado a revisar como tarea).

# Trazador cúbico natural
# Condición: S''(x_0) = S''(x_n) = 0
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def traza3natural(xi,yi):
    n = len(xi)
    
    # Valores h
    h = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        h[j] = xi[j+1] - xi[j]
    
    # Sistema de ecuaciones
    A = np.zeros(shape=(n-2,n-2), dtype = float)
    B = np.zeros(n-2, dtype = float)
    S = np.zeros(n, dtype = float)

    A[0,0] = 2*(h[0]+h[1])
    A[0,1] = h[1]
    B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])

    for i in range(1,n-3,1):
        A[i,i-1] = h[i]
        A[i,i] = 2*(h[i]+h[i+1])
        A[i,i+1] = h[i+1]
        factor21 = (yi[i+2]-yi[i+1])/h[i+1]
        factor10 = (yi[i+1]-yi[i])/h[i]
        B[i] = 6*(factor21 - factor10)
        
    A[n-3,n-4] = h[n-3]
    A[n-3,n-3] = 2*(h[n-3]+h[n-2])
    factor12 = (yi[n-1]-yi[n-2])/h[n-2]
    factor23 = (yi[n-2]-yi[n-3])/h[n-3]
    B[n-3] = 6*(factor12 - factor23)
    
    # Resolver sistema de ecuaciones S
    r = np.linalg.solve(A,B)
    for j in range(1,n-1,1):
        S[j] = r[j-1]
    S[0] = 0
    S[n-1] = 0
    
    # Coeficientes
    a = np.zeros(n-1, dtype = float)
    b = np.zeros(n-1, dtype = float)
    c = np.zeros(n-1, dtype = float)
    d = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        a[j] = (S[j+1]-S[j])/(6*h[j])
        b[j] = S[j]/2
        factor10 = (yi[j+1]-yi[j])/h[j]
        c[j] = factor10 - (2*h[j]*S[j]+h[j]*S[j+1])/6
        d[j] = yi[j]
    
    # Polinomio trazador
    x = sym.Symbol('x')
    px_tabla = []
    for j in range(0,n-1,1):

        pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
        pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
        
        pxtramo = pxtramo.expand()
        px_tabla.append(pxtramo)
    
    return(px_tabla)

# PROGRAMA -----------------------
# INGRESO , Datos de prueba
xi = np.array([0.1 , 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.8, 1.7, 2.0])
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = traza3natural(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print(' x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

con lo que se obtiene el resultado por cada tramo

Polinomios por tramos: 
 x = [0.1,0.2]
-146.666666666667*x**3 + 44.0*x**2 + 0.566666666666666*x + 1.1
 x = [0.2,0.3]
283.333333333333*x**3 - 214.0*x**2 + 52.1666666666667*x - 2.34
 x = [0.3,0.4]
-136.666666666667*x**3 + 164.0*x**2 - 61.2333333333333*x + 9.0
>>>

Para el trazado de la gráfica se añade al final del algoritmo:

# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)
    
    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)

    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'ro', label='puntos')
plt.plot(xtraza,ytraza, label='trazador'
         , color='blue')
plt.title('Trazadores Cúbicos Naturales')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

Si los polinomios no se igualan entre los nodos, tampoco sus velocidades y aceleraciones; podría considerar la experiencia semejante a variaciones de velocidad y aceleración en una trayectoria como la mostrada en los siguientes videos.

Car sales woman scares customers. maxman.tv. 4 enero 2016

4.3 Trazadores lineales (Splines) grado1 con Python

[ Trazador Lineal ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


1. Trazadores lineales (Splines) grado1

Referencia: Chapra 18.6.1 p525

El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (spline) para dibujar curvas suaves a través de un conjunto de puntos.

La unión más simple entre dos puntos es una línea recta. El método crea un polinomio para cada par de puntos consecutivos en el intervalo, por lo que el resultado será una tabla de polinomios.

Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales.

f_0(x) = f(x_0) + m_0(x-x_0), x_0\leq x\leq x_1 f_1(x) = f(x_1) + m_1(x-x_1), x_1\leq x\leq x_2

f_n(x) = f(x_{n-1}) + m_{n-1}(x-x_{n-1}), x_{n-1}\leq x\leq x_n

donde

m_i = \frac{f(x_{i+1}) - f(x_i)}{(x_{i+1}-x_i)}

Observe que la expresión de f(x) para un tramo entre dos puntos es el polinomio de grado 1 realizado con diferencia finita avanzadas  o las diferencias divididas.

Las ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y xn. Al localizar primero el intervalo dentro del cual está el punto, puede seleccionar el polinomio que corresponde a ese tramo.

[ Trazador Lineal ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


2. Ejercicio

Datos de los puntos como ejemplo para el algoritmo

xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]

[ Trazador Lineal ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

..


3. Desarrollo analítico

El método con trazadores lineales, permite plantear los bloques necesarios para manejar varios polinomios, uno por cada tramo entre dos puntos dentro del intervalo del problema

0.1 ≤ x≤ 0.2

m_0= \frac{1.8 - 1.45}{0.2-0.1} = 3.5 f_0(x) = 1.45 + 3.5(x-0.1)

0.2 ≤ x ≤ 0.3

m_1= \frac{1.7 - 1.8}{0.3-0.2} = -1 f_1(x) = 1.8 + -1(x-0.2)

0.3 ≤ x ≤ 0.4

m_2= \frac{2 - 1.7}{0.4-0.3} = 3 f_3(x) = 1.7 + 3(x-0.3)

El método permitirá disponer de un punto de partida para trazadores de mayor grado, por ejemplo los cúbicos.

[ Trazador Lineal ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


4. Algoritmo en Python

Los polinomios de cada tramo se almacenan en una tabla, cada uno puede ser utilizado individualmente en su respectivo tramo, por ejemplo para realizar la gráfica de línea entre tramos.

# Trazador (spline) lineal, grado 1
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def trazalineal(xi,fi):
    n = len(xi)
    x = sym.Symbol('x')
    px_tabla = []
    
    tramo = 1
    while not(tramo>=n):
        # con 1ra diferencia finita avanzada 
        numerador = fi[tramo]-fi[tramo-1]
        denominador = xi[tramo]-xi[tramo-1]
        m = numerador/denominador
        pxtramo = fi[tramo-1] + m*(x-xi[tramo-1])
        px_tabla.append(pxtramo)
        tramo = tramo + 1
        
    return(px_tabla)

# PROGRAMA
# INGRESO , Datos de prueba
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = trazalineal(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print('  x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

Se obtiene como resultado:

Polinomios por tramos: 
  x = [0.1,0.2]
3.5*x + 1.1
  x = [0.2,0.3]
-1.0*x + 2.0
  x = [0.3,0.4]
3.0*x + 0.8
>>> 

Para añadir la gráfica se añaden las instrucciones para:

  • evaluar el polinomio en cada tramo
  • concatenar los resultados de todos los tramos en los vectores xtraza, ytraza.
  • poner en la gráfica los puntos del problema y las líneas que genera cada polinomio
# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)

    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)

    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'o', label='puntos')
plt.plot(xtraza,ytraza, label='trazador')
plt.title('Trazadores lineales (splines)')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

[ Trazador Lineal ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

4.2 Interpolación polinómica de Lagrange con Python

[ Interpola Lagrange ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


1. Interpolación polinómica de Lagrange

Referencia: Chapra 18.2 p516, Burden 3.1 p80, Rodríguez 6.2 p195

El polinomio de interpolación de Lagrange reformula el polinomio de interpolación de Newton evitando el cálculo de la tabla de diferencias divididas. El método de Lagrange tolera las diferencias entre distancias x de los puntos de muestra.

El polinomio de Lagrange se construye a partir de las fórmulas:

f_{n} (x) = \sum_{i=0}^{n} L_{i} (x) f(x_{i}) L_{i} (x) = \prod_{j=0, j \neq i}^{n} \frac{x-x_j}{x_i - x_j}

Donde una vez que se han seleccionado los puntos a usar que generan la misma cantidad de términos que puntos.

[ Interpola Lagrange ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


2. Ejercicio

Dados los 4 puntos en la tabla se requiere un polinomio de grado 3.

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

[ Interpola Lagrange ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

..


3. Desarrollo analítico

Se calculan cuatro términos de Lagrange según las fórmulas,

término 1

L_{0} (x) = \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)}

término 2

L_{1} (x) = \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)}

término 3

L_{2} (x) = \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)}

término 4

L_{3} (x) = \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

se construye el polinomio usando la fórmula para fn(x) para cada valor fi,

p_3(x) = 1 L_{0} (x) + 1.6 L_{1} (x) + 1.7 L_{2} (x) + 2 L_{3} (x) p_3(x) = 1 \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)} + + 1.6 \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)} + 1.7 \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)} + 2 \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

La simplificación de la expresión del polinomio se puede dejar como tarea o realizarla con Sympy con la instrucción sym.expand().

p_3(x) = 41.6666 x^3 - 27.5 x^2 + 6.8333 x + 1

Una forma de verificar que el polinomio es correcto es usar un punto original xi[i] y comprobar que la evaluación tiene como resultado fi[i].

[ Interpola Lagrange ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


3. Algoritmo en Python

En el algoritmo en Python, para construir las expresiones de cada término de Lagrange se usa la forma simbólica con Sympy, con los datos:

xi = [0, 0.2, 0.3, 0.4]
fi = [1, 1.6, 1.7, 2.0]

En cada término Li(x) se usan todos los elementos i , excepto el mismo elemento i, en el numerador y denominador de la expresión.

En polinomio se agrupan todos los términos multiplicados por su respectivo valor fi[i].

    valores de fi:  [1.  1.6 1.7 2. ]
divisores en L(i):  [-0.024  0.004 -0.003  0.008]

Polinomio de Lagrange, expresiones
400.0*x*(x - 0.4)*(x - 0.3) 
  - 566.666666666667*x*(x - 0.4)*(x - 0.2) 
  + 250.0*x*(x - 0.3)*(x - 0.2) 
  - 41.6666666666667*(x - 0.4)*(x - 0.3)*(x - 0.2)

Polinomio de Lagrange: 
41.666666666667*x**3 - 27.5*x**2 + 6.8333333333334*x + 1.0

Las expresiones del polinomio contiene los binomios en su forma básica, para resolver y simplificar las ecuaciones se usa polinomio.expand().

Para realizar la gráfica del polinomio es conveniente convertirlo a la forma lambda con Numpy, de esta forma se evalúa en una sola linea todos los puntos para el intervalo [a,b] en x.

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = [0, 0.2, 0.3, 0.4]
fi = [1, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
# 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)

# 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)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

[ Interpola Lagrange ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

4.1.4 Interpolación por Diferencias divididas de Newton con Python

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


1. Diferencias divididas de Newton

Referencia: Chapra 18.1.3 p508, Rodríguez 6.7 p223, Burden 9Ed 3.3 p124.

El método se usa en el caso que los puntos en el «eje x» se encuentran espaciados de forma arbitraria y provienen de una función desconocida pero supuestamente diferenciable.

La n-ésima diferencia dividida finita es:

f[x_{n}, x_{n-1}, \text{...}, x_{1}, x_{0}] = \frac{f[x_{n}, x_{n-1}, \text{...}, x_{1}]- f[x_{n-1}, x_{n-2}, \text{...}, x_{0}]}{x_{n}-x_{0}}

Para lo cual se debe interpretar la tabla de diferencias divididas, también como una aproximación a una derivada:

i xi f[xi] Primero Segundo Tercero
0 x0 f[x0] f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
1 x1 f[x1] f[x2,x1] f[x3,x2,x1]
2 x2 f[x2] f[x3,x2]
3 x2 f[x3]

En la fórmula del polinomio, las diferencias divididas sirven para evaluar los coeficientes de cada término adicional para aumentar el grado. Semejante al proceso realizado para  «diferencias finitas divididas»:

f_n(x) = f(x_0)+(x-x_0)f[x_1,x_0] + + (x-x_0)(x-x_1)f[x_2,x_1,x_0] + \text{...}+ + (x-x_0)(x-x_1)\text{...}(x-x_{n-1})f[x_n,x_{n-1},\text{...},x_0]

Este polinomio de interpolación se conoce como de Newton en diferencias divididas.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


2. Ejercicio

Referencia: 1Eva_IIT2008_T3_MN Ganancia en inversión inversionGanancia01

Se dispone de los datos (x, f(x)), en donde x es un valor de inversión y f(x) es un valor de ganancia, ambos en miles de dólares:

inversión ganancia
3.2 5.12
3.8 6.42
4.2 7.25
4.5 6.85

Considere que los valores invertidos en materia prima para producción, dependen de la demanda de del producto en el mercado, motivo por el que los valores de inversión no guardan distancias equidistantes entre si.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

..


3. Desarrollo analítico

Se toman los datos de la tabla como arreglos para el algoritmo

xi = [3.2 , 3.8 , 4.2 , 4.5 ]
fi = [5.12, 6.42, 7.25, 6.85]

Con los datos se llena la tabla de diferencias divididas, donde por simplicidad, se escriben las operaciones en cada casilla. La última columna o cuarta diferencia dividida es cero por no disponer de datos para hacer el cálculo.

i xi f[xi] Primero Segundo Tercero
0 3.2 5.12 \frac{6.42-5.12}{3.8-3.2} =2.1667 \frac{2.075-2.1667}{4.2-3.2} =-0.0917 \frac{-4.869-(-0.0917)}{4.5-3.2} =-3.6749
1 3.8 6.42 \frac{7.25-6.42}{4.2-3.8} =2.075 \frac{-1.3333-2.075}{4.5-3.8} =-4.869
2 4.2 7.25 \frac{6.85-7.25}{4.5-4.2} =-1.3333
3 4.5 6.85

Las diferencias divididas de la primera fila son los valores usados para la expresión del polinomio de interpolación:

ddividida = tabla[0,3:] 
          = [ 2.1667, -0.0917, -3.6749, 0. ]

La expresión del polinomio inicia con el primer valor de la función f(x0).

p_3(x) = f_0 = 5.12

se calcula el primer término usando el factor de la primera diferencia dividida y se añade a la expresión del polinomio.

término = (x-x_0)f[x1,x0] = (x-3.2)2.1667 p_3(x) = 5.12 + 2.1667(x-3.2)

con el siguiente valor de ddividida[] se procede de manera semejante:

término = (x-x_0)(x-x_1)f[x1,x0] = = (x-3.2)(x-3.8)(-0.0917) p_3(x) = 5.12 + 2.1667(x-3.2)+ + (-0.0917)(x-3.2)(x-3.8)

Realizando el cálculo del tercer y último termino con diferencias divididas,  el polinomio obtenido es:

p_3(x) = 5.12 + 2.1667(x-3.2) + + (-0.0917)(x-3.2)(x-3.8) + + (-3.6749)(x - 3.2)(x - 3.8)(x - 4.2)

que simplificando al multiplicar entre los términos (x-xi):

p_3(x) = 184.7569 - 149.9208 x + 41.0673 x^2 - 3.6749 x^3

El polinomio en el intervalo de xi, en la gráfica muestra que pasa por todos los puntos.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


4. Algoritmo en Python

El algoritmo para  interpolación de Diferencias Divididas o Newton, considera reutilizar el procedimiento de cálculo de diferencias finitas incorporando la parte del denominador.

La creación de la expresión del polinomio también es semejante a la usada para diferencias finitas avanzadas.

Se incorpora la parte gráfica para observar los resultados en el intervalo xi, con el número de muestras = 101 para tener una buena resolución de la línea del polinomio.

# 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
xi = [3.2, 3.8, 4.2, 4.5]
fi = [5.12, 6.42, 7.25, 6.85]

# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)

# Tabla de Diferencias Divididas
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)

# 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)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()

El resultado del algoritmo se muestra a continuación:

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]']]
[[ 0.      3.2     5.12    2.1667 -0.0917 -3.6749  0.    ]
 [ 1.      3.8     6.42    2.075  -4.869   0.      0.    ]
 [ 2.      4.2     7.25   -1.3333  0.      0.      0.    ]
 [ 3.      4.5     6.85    0.      0.      0.      0.    ]]
dDividida: 
[ 2.1667 -0.0917 -3.6749  0.    ]
polinomio: 
2.16666666666667*x - 3.67490842490842*(x-4.2)*(x-3.8)*(x-3.2)
 - 0.0916666666666694*(x - 3.8)*(x - 3.2) - 1.81333333333334
polinomio simplificado: 
-3.67490842490842*x**3 + 41.0673076923077*x**2 
- 149.920860805861*x + 184.756923076923
>>> 

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

4.1.2 Interpolación por Diferencias finitas avanzadas con Python

[ Dif_Finitas ] || [ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] || [ Dif_Divididas_Newton ]
..


1. Interpolación por Diferencias finitas avanzadas

Referencia: Rodríguez 6.6.4 p221, Burden 9Ed p129

Se usa en interpolación cuando los puntos en el «eje x» se encuentran igualmente espaciados, la diferencia entre puntos consecutivos xi es una constante denominada h.

h = xi+1 – xi


Una relación entre derivadas y diferencias finitas se establece mediante:

f^{(n)}(z) = \frac{\Delta ^{n} f_{0}}{h^{n}} para algún k en el intervalo [x0,xn]

\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]

El polinomio de interpolación se puede construir por medio de diferencias finitas avanzadas con las siguiente fórmula:

p_n (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...} + \frac{\Delta^n f_0}{n!h^n} (x - x_0)(x - x_1) \text{...} (x - x_{n-1})

Observe que en la medida que se toman más puntos de muestra, el grado del polinomio puede ser mayor.

[ Dif_Finitas ] || [ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] || [ Dif_Divididas_Newton ]
..


2. Ejercicio

Se toman los datos del ejercicio de diferencias finitas , observando que se requiere que el tamaño de paso h  sea constante entre los puntos consecutivos xi.

xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

Para éste ejercicio se comprueba que h = 0.1

[ Dif_Finitas ] || [ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] || [ Dif_Divididas_Newton ]

..


3. Desarrollo analítico

Se usan los resultados previos del ejercicio con diferencias finitas:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]

Valores con los que se inicia el cálculo de la distancia entre puntos xi, que debe ser constante:

h = xi[1] – xi[0] = 0.2-0.1 = 0.1
h = xi[2] – xi[1] = 0.3-0.2 = 0.1
h = xi[3] – xi[2] = 0.4-0.3 = 0.1

Para la construcción del polinomio, se usan los valores de diferencias finitas de la primera fila desde la columna 3 en adelante.

dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]

Se empieza con el valor del primer punto de la función f(0.1)

p3(x) = f0 = 1.45

para añadirle el primer término j=1, más el factor por calcular:

factor = \frac{\Delta^1 f_0}{1!h^1} = \frac{0.15}{1!(0.1)} = 1.5

completando el término como:

término = factor(x-x_0) = 1.5(x-0.1) p_3(x) = 1.45 + 1.5(x-0.1)

Para el segundo término j=2, se repite el proceso:

factor = \frac{\Delta^2 f_0}{2!h^2} = \frac{-0.05}{2!(0.1)^2} = -2.5 término = factor(x-x_0) = -2.5(x-0.1)(x-02) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-02)

Finalmente, añadiendo el tercer término j=3 cuyo cálculo es semejante a los anteriores, se deja como tarea.

El resultado del método es:

p_3(x)= 1.45 + 1.5(x-0.1) -2.5(x-0.1)(x-02) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Se puede seguir simplificando la respuesta, por ejemplo usando solo el término de grado con 1.5(x-0.1) se tiene que:

p_3(x) = 1.3 + 1.5x -2.5(x-0.1)(x-02)+ + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Seguir simplificando la expresión en papel, se deja como tarea.

La gráfica del polinomio obtenido comprueba que pasa por cada uno de los puntos dados para el ejercicio.

[ Dif_Finitas ] || [ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] || [ Dif_Divididas_Newton ]
..


4. Algoritmo en Python

El polinomio se construye usando el ejercicio de diferencias finitas.

Para construir la expresión del polinomio añadiendo los términos de la fórmula, se define la variable simbólica x con Sympy.

Para simplificar el polinomio resultante con las expresiones de multiplicación, se utiliza la instrucción sym.expand().

En caso de requerir evaluar la fórmula con un vector de datos se la convierte a la forma lambda.

Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores xi, definiendo el número de muestras = 101 que son suficientes para que la gráfica se observe sin distorsión.

# Polinomio interpolación
# Diferencias finitas avanzadas
# 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
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(xi,dtype=float)
# Tabla de Diferencias Finitas
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 finitas 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('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = 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):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    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)

# 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('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()

teniendo como resultado:

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]
dfinita: 
[ 0.15 -0.05  0.25  0.  ]
polinomio: 
1.5*x + 41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1) 
- 2.50000000000001*(x - 0.2)*(x - 0.1) + 1.3
polinomio simplificado: 
41.6666666666667*x**3 - 27.5*x**2 
+ 6.83333333333335*x + 0.999999999999999

el polinomio de puede evaluar como px(valor) una vez que se convierte a la forma lambda para usar con Numpy:

>>> px(0.1)
1.4500000000000004
>>> px(0.2)
1.6000000000000025
>>> px0.3)
1.7000000000000042
>>> 

Tarea: se recomienda realizar las gráficas comparativas entre métodos, debe mostrar la diferencia con los métodos que requieren el tamaño de paso equidistante h, y los que no lo requieren. Permite detectar errores de selección de método para interpolación.

[ Dif_Finitas ] || [ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] || [ Dif_Divididas_Newton ]