5.1.1 Diferencias finitas avanzadas polinomio

Referencia: Rodriguez 6.6.4 Pdf.221

En la tabla de diferencias finitas,  si los valores en el eje x se encuentran igualmente espaciados, la diferencia es una constante denominada 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 f(n) en el intervalo [x0,xn]

Polinomio de interpolación de diferencias finitas avanzadas:

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})

El polinomio se puede realizar usando el ejercicio de diferencias finitas. Se crea la variable simbólica con sympy, se genera el polinomio añadiendo los términos de la operación.

Para simplificar se expanden las operaciones de multiplicación, y se obtiene el polinomio.

En caso de requerir evaluar la fórmula con un vector de datos, de preferencia se la convierte a la forma lambda.

# Diferencias finitas avanzadas para polinomio interpolación
# Referencia Rodriguez 6.6.4 Pdf.221
# Tarea: Verificar tamaño de vectores
#        considerar puntos no equidistantes en eje x
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# 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)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
xi = np.transpose([xi])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
    # Aumenta el título para cada columna
    titulo.append('df'+str(c-2))
    # calcula cada diferencia por fila
    f = 0
    while (f < diagonal):
        tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
        f = f+1
    
    diagonal = diagonal - 1
    c = c+1

# POLINOMIO con diferencias finitas
# caso: puntos en eje x equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
x = sym.Symbol('x')
h = xi[1,0]-xi[0,0]
polinomio = fi[0,0]
for c in range(1,n,1):
    denominador = np.math.factorial(c)*(h**c)
    factor = dfinita[0,c-1]/denominador
    termino=1
    for f  in range(0,c,1):
        termino = termino*(x-xi[f])
    polinomio = polinomio + termino*factor
# simplifica polinomio, multiplica los (x-xi)
polinomio = polinomio.expand()
# para evaluacion numérica
px = sym.lambdify(x,polinomio)

# Puntos para la gráfica
a = np.min(xi)
b = np.max(xi)
muestras = 101
xi_p = np.linspace(a,b,muestras)
fi_p = px(xi_p)

# SALIDA
print([titulo])
print(tabla)
print('polinomio:')
print(polinomio)

# Gráfica
plt.title('Interpolación polinómica')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(xi_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()

teniendo como resultado:

[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 1.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 2.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 3.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 4.    0.4   2.    0.    0.    0.    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 expandido
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333335*x + 0.999999999999999
>>> 

el polinomio de puede evaluar como p(valor) una vez que se convirtió a lambda:

>>> p(0.1)
1.4500000000000004
>>> p(0.2)
1.6000000000000025
>>> p(0.3)
1.7000000000000042
>>> 

Como tarea se recomienda realizar las gráficas comparativas entre métodos: