4.1.4 Interpolación por Diferencias divididas de Newton con Python

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

1. Diferencias divididas de Newton

Referencia: Chapra 18.1.3 p508, Rodríguez 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	x_i	f[x_i]	Primero	Segundo	Tercero
0	x₀	f[x₀]	f[x₁,x₀]	f[x₂,x₁,x₀]	f[x₃,x₂,x₁,x₀]
1	x₁	f[x₁]	f[x₂,x₁]	f[x₃,x₂,x₁]
2	x₂	f[x₂]	f[x₃,x₂]
3	x₂	f[x₃]

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.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

2. Ejercicio

Referencia: 1Eva_IIT2008_T3_MN Ganancia en inversión

Se dispone de los datos (x, f(x)), en donde x es un valor de inversión y f(x) es un valor de ganancia, ambos en miles de dólares:

inversión	ganancia
3.2	5.12
3.8	6.42
4.2	7.25
4.5	6.85

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 equidistantes entre si.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

3. Desarrollo analítico

Se toman los datos de la tabla como arreglos para el algoritmo

xi = [3.2 , 3.8 , 4.2 , 4.5 ]
fi = [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 columna o cuarta diferencia dividida es cero por no disponer de datos para hacer el cálculo.

i	x_i	f[x_i]	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$	$\frac{-4.869-(-0.0917)}{4.5-3.2} =-3.6749$
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 interpolación:

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

se 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^3

El polinomio en el intervalo de xi, en la gráfica muestra que pasa por todos los puntos.

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

4. 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 polinomio.

# 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 = [3.2, 3.8, 4.2, 4.5]
fi = [5.12, 6.42, 7.25, 6.85]

# PROCEDIMIENTO
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
        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
>>>

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] ||
[ Dif_Divididas_Newton ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]