Referencia: Chapra 18.1.3 p508 pdf532, Rodriguez 6.7 p223, Burden 9Ed 3.3 p124.
El método se usa en el caso que los puntos en el «eje x» se encuentran espaciados de forma arbitraria y provienen de una función desconocida pero supuestamente diferenciable.
La n-ésima diferencia dividida finita es:
f[x_{n}, x_{n-1}, \text{...}, x_{1}, x_{0}] = \frac{f[x_{n}, x_{n-1}, \text{...}, x_{1}]- f[x_{n-1}, x_{n-2}, \text{...}, x_{0}]}{x_{n}-x_{0}}Para lo cual se debe interpretar la tabla de diferencias divididas, también como una aproximación a una derivada:
| i | xi | f[xi] | Primero | Segundo | Tercero |
| 0 | x0 | f[x0] | f[x1,x0] | f[x2,x1,x0] | f[x3,x2,x1,x0] |
| 1 | x1 | f[x1] | f[x2,x1] | f[x3,x2,x1] | |
| 2 | x2 | f[x2] | f[x3,x2] | ||
| 3 | x2 | f[x3] |
En la fórmula del polinomio, las diferencias divididas sirven para evaluar los coeficientes de cada término adicional para aumentar el grado. Semejante al proceso realizado para «diferencias finitas divididas»:
f_n(x) = f(x_0)+(x-x_0)f[x_1,x_0] + + (x-x_0)(x-x_1)f[x_2,x_1,x_0] + \text{...}+ + (x-x_0)(x-x_1)\text{...}(x-x_{n-1})f[x_n,x_{n-1},\text{...},x_0]Este polinomio de interpolación se conoce como de Newton en diferencias divididas.
Ejercicio
El ejercicio perminte mostrar con mayor detalle las operaciones del procedimiento. Para el ejercicio, se toma como referencia el descrito en el enlace, observando que las distancias entre puntos para el «eje x» NO son equidistantes,
Referencia: 1Eva_IIT2008_T3_MN Ganancia en inversión
Los datos del problema muestran varios valores de inversión xi que generan cada uno ganancias fi.
Considere que los valores invertidos en materia prima para producción, dependen de la demanda de del producto en el mercado, motivo por el que los valores de inversión no guardan distancias equidistante entre si.
Siguiento del formato usado en los algoritmos, los datos son:
xi = np.array([3.2 , 3.8 , 4.2 , 4.5 ]) fi = np.array([5.12, 6.42, 7.25, 6.85])
Con los datos se llena la tabla de diferencias divididas, donde por simplicidad, se escriben las operaciones en cada casilla. La última de deja como tarea.:
| i | xi | f[xi] | Primero | Segundo | Tercero |
| 0 | 3.2 | 5.12 | \frac{6.42-5.12}{3.8-3.2} =2.1667 | \frac{2.075-2.1667}{4.2-3.2} =-0.0917 | f[x3,x2,x1,x0] |
| 1 | 3.8 | 6.42 | \frac{7.25-6.42}{4.2-3.8} =2.075 | \frac{-1.3333-2.075}{4.5-3.8} =-4.869 | |
| 2 | 4.2 | 7.25 | \frac{6.85-7.25}{4.5-4.2} =-1.3333 | ||
| 3 | 4.5 | 6.85 |
Las diferencias divididas de la primera fila son los valores usados para la expresión del polinomio de interpolacion:
ddividida = tabla[0,3:] = [ 2.1667, -0.0917, -3.6749, 0. ]
La expresión del polinomio inicia con el primer valor de la función f(x0).
p_3(x) = f_0 = 5.12se calcula el primer término usando el factor de la primera diferencia dividida y se añade a la expresión del polinomio.
término = (x-x_0)f[x1,x0] = (x-3.2)2.1667 p_3(x) = 5.12 + 2.1667(x-3.2)con el siguiente valor de ddividida[] se procede de manera semejante:
término = (x-x_0)(x-x_1)f[x1,x0] = = (x-3.2)(x-3.8)(-0.0917) p_3(x) = 5.12 + 2.1667*(x-3.2)+ + (-0.0917)(x-3.2)(x-3.8)Realizando el cálculo del tercer y último termino con diferencias divididas, el polinomio obtenido es:
p_3(x) = 5.12 + 2.1667(x-3.2) + + (-0.0917)(x-3.2)(x-3.8) + + (-3.6749)(x - 3.2)(x - 3.8)(x - 4.2)que simplificando al multiplicar entre los términos (x-xi):
p_3(x) = 184.7569 - 149.9208 x + 41.0673 x^2 - 3.6749 x^3El polinomio graficado en el intervalo de xi muestra que pasa por todos los puntos.
Algoritmo en Python
El algoritmo para interpolación de diferencias divididas o Newton, considera reutilizar el procedimiento de cálculo de diferencias finitas incorporando la parte del denominador.
La creación de la expresión del polinomio también es semejante a la usada para diferencias finitas avanzadas.
Se incorpora la parte gráfica para observar los resultados en el intervalo xi, con el número de muestras = 101 para tener una buena resolución de la línea del polinomo.
# Polinomio interpolación # Diferencias Divididas de Newton # 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([3.2, 3.8, 4.2, 4.5]) fi = np.array([5.12, 6.42, 7.25, 6.85]) # PROCEDIMIENTO # Tabla de Diferencias Divididas Avanzadas 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 i = i+1 diagonal = diagonal - 1 j = j+1 # POLINOMIO con diferencias Divididas # caso: puntos equidistantes en eje x dDividida = 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): factor = dDividida[j-1] 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 np.set_printoptions(precision = 4) print('Tabla Diferencia Dividida') print([titulo]) print(tabla) print('dDividida: ') print(dDividida) 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('Diferencias Divididas - Newton') plt.show()
El resultado del algoritmo se muestra a continuación:
Tabla Diferencia Dividida [['i ', 'xi ', 'fi ', 'F[1]', 'F[2]', 'F[3]', 'F[4]']] [[ 0. 3.2 5.12 2.1667 -0.0917 -3.6749 0. ] [ 1. 3.8 6.42 2.075 -4.869 0. 0. ] [ 2. 4.2 7.25 -1.3333 0. 0. 0. ] [ 3. 4.5 6.85 0. 0. 0. 0. ]] dDividida: [ 2.1667 -0.0917 -3.6749 0. ] polinomio: 2.16666666666667*x - 3.67490842490842*(x-4.2)*(x-3.8)*(x-3.2) - 0.0916666666666694*(x - 3.8)*(x - 3.2) - 1.81333333333334 polinomio simplificado: -3.67490842490842*x**3 + 41.0673076923077*x**2 - 149.920860805861*x + 184.756923076923 >>>