Referencia: Rodriguez 6.6.4 Pdf.221, 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 k en el intervalo [x0,xn]
\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]
El polinomio de interpolación se puede contruir 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.
Ejercicio
Se toman los datos del ejercicio de «diferencias finitas» y se grafican los resultados. Se resalta en la gráfica la distancia equidistante h entre los puntos consecutivos xi como los mostrados en la gráfica anterior.
xi = np.array([0.10, 0.2, 0.3, 0.4]) fi = np.array([1.45, 1.6, 1.7, 2.0])
y se usan los resulados previos de diferencias finitas:
Tabla Diferencia Finita: [['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. ]]
valores con los que se inicia el cálculo de la distancia constante entre puntos xi:
h = xi[1] – xi[0] = 0.2-0.1 = 0.1
Los valores de diferencias finitas a usar en la fórmula son de la primera fila desde la columna 3 en adelante.
dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]
La construción de la expresión del polinomio empieza con el valor del primer punto de la función f(0.1)
p_3(x) = f_0 = 1.45
para añadirle el primer término j=1, hay que calcular el factor:
factor = \frac{\Delta^1 f_0}{1!h^1} = \frac{0.15}{0.1} = 1.5con lo que el término se completa 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-02) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-02)Finalmente, añadiendo el tercer 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-02) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)Se puede seguir simpificando 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-02)+ + 41.667(x - 0.3)(x - 0.2)(x - 0.1)Seguir simplificando la expresión en papel, se deja como tarea.
La gráfica del polinomio obtenido comprueba que pasa por cada uno de los puntos dados para el ejercicio.
Algoritmo en Python
El polinomio se contruye usando el ejercicio de diferencias finitas.
Para contruir 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.
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 = 101 que son suficientes para que la gráfica se observe sin distorsión.
# Polinomio interpolación # Diferencias finitas avanzadas # Tarea: Verificar tamaño de vectores, # verificar puntos equidistantes en x import numpy as np import sympy as sym import matplotlib.pyplot as plt # INGRESO , Datos de prueba xi = np.array([0.10, 0.2, 0.3, 0.4]) fi = np.array([1.45, 1.6, 1.7, 2.0]) # PROCEDIMIENTO # 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 # POLINOMIO con diferencias Finitas avanzadas # caso: puntos equidistantes en eje x h = xi[1] - xi[0] dfinita = tabla[0,3:] n = len(dfinita) # expresión del polinomio con Sympy x = sym.Symbol('x') polinomio = fi[0] for j in range(1,n,1): denominador = np.math.factorial(j)*(h**j) factor = dfinita[j-1]/denominador termino = 1 for k in range(0,j,1): termino = termino*(x-xi[k]) polinomio = polinomio + termino*factor # simplifica multiplicando entre (x-xi) polisimple = polinomio.expand() # polinomio para evaluacion numérica px = sym.lambdify(x,polisimple) # Puntos para la gráfica muestras = 101 a = np.min(xi) b = np.max(xi) pxi = np.linspace(a,b,muestras) pfi = px(pxi) # SALIDA print('Tabla Diferencia Finita') print([titulo]) print(tabla) print('dfinita: ') print(dfinita) print('polinomio: ') print(polinomio) print('polinomio simplificado: ' ) print(polisimple) # Gráfica plt.plot(xi,fi,'o', label = 'Puntos') ##for i in range(0,n,1): ## plt.axvline(xi[i],ls='--', color='yellow') plt.plot(pxi,pfi, label = 'Polinomio') plt.legend() plt.xlabel('xi') plt.ylabel('fi') plt.title('Interpolación polinómica') plt.show()
teniendo como resultado:
Tabla Diferencia Finita [['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. ] 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
el polinomio de puede evaluar como px(valor) una vez que se convirtió a lambda:
>>> px(0.1) 1.4500000000000004 >>> px(0.2) 1.6000000000000025 >>> px0.3) 1.7000000000000042 >>>
Como tarea se recomienda realizar las gráficas comparativas entre métodos: