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. ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
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. ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
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. ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
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.6667⋅x - 27.5⋅x + 6.8333⋅x + 1El resumen de las instrucciones se presentan a continuación.
# Interpolación grafica animada.gif
# - Diferencias finitas avanzadas
# - Diferencias finitas divididas
# Tarea: Verificar tamaño de vectores,
import numpy as np
import math
import sympy as sym
def revisa_orden(xi,orden='up'):
''' Revisa orden en vector xi.
orden='up': orden ascendente
orden='down' : orden descendente
resultado:
True: xi ordenado,
False: xi desordenado y dónde.
'''
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
n = len(xi)
msj = [] # mensajes de error
# xi revisar orden ascendente o descendente
tramos = np.diff(xi,1) # diferencias en xi
ordenado = True # suponer ordenados
k = 1 # 1: ascendente
if orden=='down':
k = -1 # -1: descendente
donde = -1
i = 0
while i<(n-1) and donde<0:
if k*tramos[i]<0:
ordenado = False
donde = i+1
msj.append(['sin orden',donde])
i = i+1
return([ordenado,msj])
def tramos_equidistantes(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.
'''
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
n = len(xi)
msj = [] # mensajes de error
# revisa tamaños de paso desiguales o no equidistantes
h_iguales = True
if (n-1)>=2: # al menos dos tramos
dtramos = np.diff(xi,2) # cero o menores a casicero
errado = np.max(np.abs(dtramos)) # |error| mayor
if errado>=casicero: # no equidistantes
h_iguales=False
donde = np.argmax(np.abs(dtramos))+1
msj.append(['no equidistantes',donde])
return([h_iguales,msj])
def diferencias_tabla(xi,fi,tipo='finitas', vertabla=False,
casicero = 1e-15, precision=4):
'''Genera la tabla de diferencias finitas o divididas
tipo = 'finitas, tipo = 'divididas'
resultado en: tabla, titulo
Tarea: verificar tamaño de vectores
'''
# revisa tipo de tabla
tipolist = ['finitas','divididas']
if not(tipo in tipolist):
print('error de tipo, seleccione:',tipolist)
return()
prefijo = ['d','f']
if tipo=='divididas':
prefijo = ['F[',']']
# Matrices como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Tabla de Diferencias Finitas/Divididas
tabla_etiq = ['i','xi','fi']
ki = np.arange(0,n,1) # filas
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)
n,m = np.shape(tabla)
# Calcular tabla de Diferencias Finitas/Divididas
diagonal = n-1 # derecha a izquierda
j = 3 # inicia en columna 3
while (j < m): # columna
tabla_etiq.append(prefijo[0]+str(j-2)+prefijo[1])
# paso = j-2 # inicia en 1
i = 0 # fila
while (i < diagonal): # antes de diagonal
tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
if tipo=='divididas':
# denominador = xi[i+paso]-xi[i]
denominador = xi[i+j-2]-xi[i]
tabla[i,j] = tabla[i,j]/denominador
if abs(tabla[i,j])<casicero: # casicero revisa
tabla[i,j]=0
i = i + 1
diagonal = diagonal - 1
j = j + 1
if vertabla==True:
np.set_printoptions(precision)
print('Tabla de Diferencias',tipo)
print([tabla_etiq])
print(tabla)
return([tabla, tabla_etiq])
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)
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)
# PROGRAMA ----------------------------
# INGRESO , Datos de prueba
xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
tipo_tabla = 'finitas' # 'divididas'
vertabla = True
precision = 4
# PROCEDIMIENTO
casicero = 1e-12
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
# tabla polinomios
n = len(xi)
px_tabla = [fi[0]]
for grado in range(1,n,1):
xj = xi[:grado+1]
fj = fi[:grado+1]
m = len(xj)
xi_ascendente,msj = revisa_orden(xj)
h_iguales,msj_h = tramos_equidistantes(xj, casicero = casicero)
if len(msj_h)>0:
msj.extend(msj_h)
tabla,titulo = diferencias_tabla(xj,fj,tipo=tipo_tabla,
vertabla=vertabla, casicero = casicero)
dfinita = tabla[0,3:] # diferencias finitas
h = xj[1] - xj[0] # suponer tramos equidistantes
# polinomio con diferencias Finitas Avanzadas/ divididas
x = sym.Symbol('x')
polinomio = fj[0] +0*x # sym.S.Zero en Sympy
for i in range(1,m,1):
factor = dfinita[i-1] # diferencias divididas
if tipo_tabla=='finitas': # diferencias Finitas Avanzadas
denominador = math.factorial(i)*(h**i)
factor = factor/denominador
termino = 1
for j in range(0,i,1):
termino = termino*(x-xj[j])
polinomio = polinomio + termino*factor
polisimple = polinomio.expand() # simplifica los (x-xi)
px = sym.lambdify(x,polisimple) # evaluacion numerica
try: # intenta redondear coeficientes a precisiion
polisimple = redondea_coef(polisimple, precision,casicero)
except NameError:
print("redondea_coef() no se encuentra definida.")
px_tabla.append(polisimple)
if vertabla==True:
print('xi_ascendente:',xi_ascendente)
print('tramos xi Equidistantes:',h_iguales)
if len(msj)>0: # mensajes de error
print('Revisar tramos, d1x:',np.diff(xi,1))
for unmsj in msj:
print('Tramos',unmsj[0],
'desde i:',unmsj[1])
else:
print('Tramo h:',h)
print('',)
# SALIDA
if tipo_tabla == 'finitas':
print('Polinomios con Diferencias Finitas Avanzadas')
if tipo_tabla == 'divididas':
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()
Parte adicional para la gráfica con GIF animado
# GRAFICA CON ANIMACION polinomios --------
import matplotlib.pyplot as plt
import matplotlib.animation as animation
if tipo_tabla == 'finitas':
unmetodo = 'Diferencias Finitas Avanzadas'
narchivo = 'DifFinitasAvz' # nombre archivo
if tipo_tabla == 'divididas':
unmetodo = 'Diferencias Divididas Newton'
narchivo = 'DifDividNewton' # 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()
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. ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
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. ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
Tabla de Diferencias divididas
[['i', 'xi', 'fi', 'F[1]', 'F[2]', 'F[3]', 'F[4]']]
[[ 0. 0.1 1.45 1.5 -2.5 5. 0. ]
[ 1. 0.2 1.6 1. 0. 0. 0. ]
[ 2. 0.3 1.7 1. 0. 0. 0. ]
[ 3. 0.6 2. 0. 0. 0. 0. ]]
xi_ascendente: True
tramos xi Equidistantes: False
Revisar tramos, d1x: [0.1 0.1 0.3]
Tramos no equidistantes desde i: 2
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⋅x - 5.5⋅x + 2.8⋅x + 1.22
>>> Instrucciones en Python
Se actualiza la sección de ingreso del programa, para obtener el resultado mostrado.
# PROGRAMA ----------------------------
# INGRESO , Datos de prueba
xi = [0.10, 0.2, 0.3, 0.6]
fi = [1.45, 1.6, 1.7, 2.0]
tipo_tabla = 'divididas' # 'divididas'
vertabla = True
precision = 4
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.