4.2.2 Interpolación por Diferencias divididas de Newton (Δx) con Python

Dif_Divididas_Newton

Ejercicio

Analítico

Algoritmo

gráfica

función

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.

2. Ejercicio

Referencia: 1Eva2008TII_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.

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.

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.

# Diferencias Divididas de Newton
# Tarea: verificar tamaño de vectores
import numpy as np
import sympy as sym
 
# INGRESO, Datos de prueba
xi = [3.2 , 3.8 , 4.2 , 4.5 ]
fi = [5.12, 6.42, 7.25, 6.85]
 
# 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 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)
# Calular tabla de Diferencias Divididas
diagonal = n-1 # derecha a izquierda
j = 3  # inicia en columna 3
while (j < m): # columna
    tabla_etiq.append('F['+str(j-2)+']') 
    i = 0 # fila de columna
    paso = j-2 # inicia en 1
    while (i < diagonal): # antes de diagonal
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]

        denominador = (xi[i+paso]-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

# Polinomio de Diferencias Finitas Avanzadas
n = len(xi)
x = sym.Symbol('x')
dDividida = tabla[0,3:] # diferencias Divididas
polinomio = fi[0] +0*x # en Sympy
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
    
polisimple = sym.expand(polinomio)

px = sym.lambdify(x,polisimple)

# SALIDA
print('Tabla Diferencia Dividida-Newton')
print([tabla_etiq])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polisimple: ')
print(polisimple)

El resultado del algoritmo se muestra a continuación:

Tabla Diferencia Dividida-Newton
[['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
>>>

5. Gráfica de polinomio de interpolación

# GRAFICA  --------------
import matplotlib.pyplot as plt

titulo = 'Interpolación Diferencias Divididas - Newton'
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)')

plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(titulo)
plt.show()

6. Algoritmo como función

Se reordena el algoritmo para usar funciones, la gráfica es la misma que para el algoritmo sin funciones anterior

El bloque de salida es semejante al resultado anterior.

# Diferencias finitas Avanzadas/Dividisas 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 = [3.2 , 3.8 , 4.2 , 4.5 ]
fi = [5.12, 6.42, 7.25, 6.85]
  
tipo_tabla = 'divididas' # 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)

La gráfica se genera considerando el tipo de interpolación que se ha seleccionado:

# 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--')
     
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(titulo)
plt.show()

4.2.2 Interpolación por Diferencias divididas de Newton (Δx) con Python

1. Diferencias divididas de Newton

2. Ejercicio

3. Desarrollo analítico

4. Algoritmo en Python

5. Gráfica de polinomio de interpolación

6. Algoritmo como función

Métodos Numéricos

Análisis Numérico