[ 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
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
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.