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.