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