1. Interpolación por Diferencias finitas avanzadas
Referencia: Rodríguez 6.6.4 p221, Burden 9Ed p129
Se usa en interpolación cuando los puntos en el "eje x" se encuentran igualmente espaciados, la diferencia entre puntos consecutivos xi es una constante denominada h.
h = xi+1 - xi
Una relación entre derivadas y diferencias finitas se establece mediante: f^{(n)}(z) = \frac{\Delta ^{n} f_{0}}{h^{n}} para algún z en el intervalo [x0,xn].
\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n)
El polinomio de interpolación se puede construir por medio de diferencias finitas avanzadas con las siguiente fórmula:
p_n (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...} + \frac{\Delta^n f_0}{n!h^n} (x - x_0)(x - x_1) \text{...} (x - x_{n-1})Observe que en la medida que se toman más puntos de muestra, el grado del polinomio puede ser mayor.
2. Ejercicio
Se toman los datos del ejercicio de diferencias finitas , observando que se requiere que el tamaño de paso h sea constante entre los puntos consecutivos xi.
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]| xi | 0.1 | 0.2 | 0.3 | 0.4 |
| fi | 1.45 | 1.6 | 1.7 | 2.0 |
Para éste ejercicio se comprueba que h = 0.1
3. Desarrollo Analítico
Se inicia el cálculo de la distancia entre puntos xi, comprobando que debe ser constante h, además que los valores xi deben encontrarse ordenados de forma ascendente:
h = xi[1] - xi[0] = 0.2-0.1 = 0.1
h = xi[2] - xi[1] = 0.3-0.2 = 0.1
h = xi[3] - xi[2] = 0.4-0.3 = 0.1
Se usan los resultados previos del ejercicio con diferencias finitas:
Tabla Diferencia Finita:
[['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. ]]Como los tamaños de paso en xi son constantes, h=0.1, es posible usar el método. Para la construcción del polinomio, se usan los valores de diferencias finitas de la primera fila desde la columna 3 en adelante.
dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]Se empieza con el valor del primer punto de la función f(0.1), j=0.
p3(x) = f0 = 1.45
para añadirle el término j=1, más el factor por calcular:
completando el término como:
término = factor(x-x_0) = 1.5(x-0.1) p_3(x) = 1.45 + 1.5(x-0.1)
Para el segundo término j=2, se repite el proceso:
factor = \frac{\Delta^2 f_0}{2!h^2} = \frac{-0.05}{2!(0.1)^2} = -2.5 término = factor(x-x_0) = -2.5(x-0.1)(x-0.2) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-0.2)
Finalmente, añadiendo el término j=3 cuyo cálculo es semejante a los anteriores, se deja como tarea.
El resultado del método es:
p_3(x)= 1.45 + 1.5(x-0.1) -2.5(x-0.1)(x-0.2) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)Se puede seguir simplificando la respuesta, por ejemplo usando solo el término de grado con 1.5(x-0.1) se tiene que:
p_3(x) = 1.3 + 1.5x -2.5(x-0.1)(x-0.2)+ + 41.667(x - 0.3)(x - 0.2)(x - 0.1)Seguir simplificando la expresión en papel, los detalles se deja como tarea, se obtiene:
p_3(x) = 1+ 6.833x 27.5x^2+41.666x^3La gráfica del polinomio obtenido comprueba que pasa por cada uno de los puntos del ejercicio.
4. Algoritmo en Python
El polinomio se construye usando el ejercicio de diferencias finitas.
Para construir la expresión del polinomio añadiendo los términos de la fórmula, se define la variable simbólica x con Sympy.
Para simplificar el polinomio resultante con las expresiones de multiplicación, se utiliza la instrucción sym.expand().
En caso de requerir evaluar la fórmula con un vector de datos se la convierte a la forma lambda para evaluación numérica.
teniendo como resultado:
dtramos: [0.1 0.1 0.1]
Tabla 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. ]
Diferencias Finitas Avanzadas
polinomio:
1.5*x + 41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1) - 2.50000000000001*(x - 0.2)*(x - 0.1) + 1.3
polinomio simplificado:
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333335*x + 0.999999999999999
p(x):
3 2
41.6666666666667⋅x - 27.5⋅x + 6.83333333333335⋅x + 0.999999999999999Instrucciones en Python
# Polinomio 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
# INGRESO
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
titulo = 'Diferencias Finitas Avanzadas'
# PROCEDIMIENTO
casicero = 1e-15 # redondear a cero
# Matrices como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Tabla de Diferencias Finitas
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)
# Calular tabla de Diferencias Finitas
diagonal = n-1 # derecha a izquierda
j = 3 # inicia en columna 3
while (j < m): # columna
tabla_etiq.append('d'+str(j-2)+'f')
i = 0 # fila
while (i < diagonal): # antes de diagonal
tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
if abs(tabla[i,j])<casicero: # casicero revisa
tabla[i,j]=0
i = i + 1
diagonal = diagonal - 1
j = j + 1
# POLINOMIO con diferencias Finitas Avanzadas
# caso: puntos equidistantes en eje x
dtramos = np.diff(xi,1) # tramos xi
dfinita = tabla[0,3:] # diferencias finitas
h = xi[1] - xi[0] # suponer tramos equidistantes
# polinomio con diferencias Finitas Avanzadas
x = sym.Symbol('x')
polinomio = fi[0] +0*x # sym.S.Zero en Sympy
for i in range(1,n,1):
denominador = math.factorial(i)*(h**i)
factor = dfinita[i-1]/denominador
termino = 1
for j in range(0,i,1):
termino = termino*(x-xi[j])
polinomio = polinomio + termino*factor
polisimple = polinomio.expand() # simplifica los (x-xi)
px = sym.lambdify(x,polisimple) # evaluacion numerica
# SALIDA
print('dtramos:',dtramos)
print('Tabla Diferencias Finitas')
print([tabla_etiq])
print(tabla)
print('dfinita: ')
print(dfinita)
print(titulo)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)
print('p(x):')
sym.pprint(polisimple)
el polinomio de puede evaluar como px(valor) una vez que se convierte a la forma lambda para usar con Numpy:
>>> px(0.1)
1.4500000000000004
>>> px(0.2)
1.6000000000000025
>>> px0.3)
1.7000000000000042
>>> 5. Gráfica de polinomio de interpolación
Las instrucciones son semejantes a las presentadas en polinomio de interpolación.
Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores xi, definiendo el número de muestras = 21 que son suficientes para que la gráfica se observe sin distorsión.
# GRAFICA --------------
import matplotlib.pyplot as plt
muestras = 21 # resolución gráfica
a = np.min(xi) # intervalo [a,b]
b = np.max(xi)
xk = np.linspace(a,b,muestras)
yk = px(xk)
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xk,yk, label='p(x)')
# entorno de grafica
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(titulo)
plt.tight_layout()
plt.show()
Tarea: se recomienda realizar las gráficas comparativas entre métodos, debe mostrar la diferencia con los métodos que requieren el tamaño de paso equidistante h, y los que no lo requieren. Permite detectar errores de selección de método para interpolación.
6. Algoritmo como función
Se reordena el algoritmo para usar funciones. El bloque de salida es semejante al resultado anterior. Las funciones pueden ser reutilizadas de ser necesarias en el siguiente método de interpolación.
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
dfinitas:
[ 0.15 -0.05 0.25 0. ]
polinomio:
1.30000000000000 +
1.5*x +
-2.50000000000001*(x - 0.2)*(x - 0.1) +
41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1)
polinomio simplificado:
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333335*x + 0.999999999999999
p(x):
3 2
41.6666666666667⋅x - 27.5⋅x + 6.83333333333335⋅x + 0.999999999999999Se añade un bloque en forma de funciones para la revisión que los puntos xi se presenten en orden ascendente y sean equidistantes.
El método para las tablas de diferencias finitas es muy semejante al de diferencias finitas divididas, se integra en una sola función con el parámetro de selección de tipo_tabla, parámetro usado también para generar el título de la gráfica.
En caso de existir inconsistencias o errores se almacenan en una lista de mensajes "msj".
Instrucciones en Python
# Diferencias finitas Avanzadas/Divididas de Newton
# funciones diferencias_tabla, revisa_orden, tramos_equidistantes,
# Tarea: Verificar tamaño de vectores
import numpy as np
import math
import sympy as sym
# INGRESO
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
tipo_tabla = 'finitas' # finitas o divididas
# Algoritmos como Funciones
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])
# PROGRAMA ---------------------
# PROCEDIMIENTO
casicero = 1e-12
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
xi_ascendente,msj = revisa_orden(xi)
h_iguales,msj_h = tramos_equidistantes(xi, casicero = casicero)
if len(msj_h)>0:
msj.extend(msj_h) # mensajes de error unificados
tablaDif = diferencias_tabla(xi,fi,tipo=tipo_tabla,
vertabla=True, casicero = casicero)
tabla = tablaDif[0]
tabla_etiq = tablaDif[1]
dfinita = tabla[0,3:] # diferencias finitas
h = xi[1] - xi[0] # suponer tramos equidistantes
# polinomio con diferencias Finitas Avanzadas/Divididas
x = sym.Symbol('x')
polinomio = fi[0] +0*x # sym.S.Zero en Sympy
for i in range(1,n,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-xi[j])
polinomio = polinomio + termino*factor
polisimple = polinomio.expand() # simplifica los (x-xi)
px = sym.lambdify(x,polisimple) # evaluacion numerica
# SALIDA
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('d'+tipo_tabla+': ')
print(dfinita)
print('polinomio: ')
#print(polinomio)
terminos = sym.Add.make_args(polinomio)
n_term = len(terminos)
for i in range(0,n_term,1):
if i<(n_term-1):
print(terminos[i],'+')
else:
print(terminos[i])
print()
print('polinomio simplificado: ' )
print(polisimple)
print('p(x):')
sym.pprint(polisimple)
Como el bloque de gráfica se usa en otros métodos de la unidad, se incorpora la verificación de lista de errores 'msj' usando el bloque 'try-except'.
La gráfica se genera considerando el tipo de interpolación que se ha seleccionado. También se ajusta el título acorde al tipo de tabla seleccionada en el bloque de ingreso.
# GRAFICA --------------
import matplotlib.pyplot as plt
muestras = 21 # resolución gráfica
titulo = 'Interpolación: Diferencias '
if tipo_tabla == 'finitas':
titulo = titulo + 'Finitas Avanzadas'
if tipo_tabla == 'divididas':
titulo = titulo + 'Divididas de Newton'
a = np.min(xi) # intervalo [a,b]
b = np.max(xi)
xk = np.linspace(a,b,muestras)
yk = px(xk)
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xk,yk, label='p(x)')
try: # existen mensajes de error
msj_existe = len(msj)
except NameError:
msj = []
if len(msj)>0 and tipo_tabla == 'finitas': # tramos con error
untipo = msj[0][0]
donde = msj[0][1] # indice error
plt.plot(xi[donde:donde+2],fi[donde:donde+2],'ro',label=untipo)
plt.plot(xi[donde:donde+2],fi[donde:donde+2],'r--')
# entorno de grafica
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(titulo)
plt.tight_layout()
plt.show()