[ 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 la librería 'pillow'.
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
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', 'd1f', 'd2f']
[[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', 'd1f', 'd2f', 'd3f']
[[ 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', 'd1f', 'd2f', 'd3f', 'd4f']
[[ 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 math 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('d'+str(j-2)+'f') # 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 = 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='pillow') # ani.save(narchivo+'_video.mp4') plt.draw() plt.show()
[ Dif_Finitas ] [ Dif_Finitas avanzadas ] [ Dif_Divididas_Newton ]
Interpolación por Diferencias Divididas de Newton
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.