4.2.1 Interpolación por Diferencias finitas avanzadas (h) con Python

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

1. Interpolación por Diferencias finitas avanzadas

Referencia: Rodríguez 6.6.4 p221, Burden 9Ed p129

Se usa en interpolación cuando los puntos en el "eje x" se encuentran igualmente espaciados, la diferencia entre puntos consecutivos x_i es una constante denominada h.

h = x_i+1 - x_i

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 [x₀,x_n].

$\frac{\Delta ^{n} f_{0}}{h^{n}}$ es una aproximación para la n-ésima derivada f⁽ⁿ⁾

El polinomio de interpolación se puede construir 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.

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

2. Ejercicio

Se toman los datos del ejercicio de diferencias finitas , observando que se requiere que el tamaño de paso h sea constante entre los puntos consecutivos x_i.

xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

x_i	0.1	0.2	0.3	0.4
f_i	1.45	1.6	1.7	2.0

Para éste ejercicio se comprueba que h = 0.1

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

3. Desarrollo Analítico

Se inicia el cálculo de la distancia entre puntos xi, comprobando que debe ser constante h, además que los valores xi deben encontrarse ordenados de forma ascendente:

h = xi[1] - xi[0] = 0.2-0.1 = 0.1
h = xi[2] - xi[1] = 0.3-0.2 = 0.1
h = xi[3] - xi[2] = 0.4-0.3 = 0.1

Se usan los resultados previos del ejercicio con diferencias finitas:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'd1f', 'd2f', 'd3f', 'd4f']]
[[ 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.  ]]

Como los tamaños de paso en xi son constantes, h=0.1, es posible usar el método. Para la construcción del polinomio, se usan los valores de diferencias finitas de la primera fila desde la columna 3 en adelante.

dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]

Se empieza con el valor del primer punto de la función f(0.1), j=0.

p₃(x) = f₀ = 1.45

para añadirle el término j=1, más el factor por calcular:

factor = \frac{\Delta^1 f_0}{1!h^1} = \frac{0.15}{1!(0.1)} = 1.5

completando el término 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-0.2)

p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-0.2)

Finalmente, añadiendo el 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-0.2) +

+ 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Se puede seguir simplificando 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-0.2)+

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

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

4. Algoritmo en Python

El polinomio se construye usando el ejercicio de diferencias finitas.

Para construir 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 para evaluación numérica.

Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores x_i, definiendo el número de muestras = 21 que son suficientes para que la gráfica se observe sin distorsión.

Se añade un bloque para la revisión que los puntos xi sean equidistantes y se encuentren en orden.

# Diferencias finitas Avanzadas - Interpolación
# Tarea: Verificar tamaño de vectores
import numpy as np
import math
import sympy as sym

# INGRESO
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

# PROCEDIMIENTO
casicero = 1e-15
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)

msj = [] # mensajes de error

# xi revisar orden ascendente
tramos  = np.diff(xi,1)  # diferencias en xi
donde = -1  # suponer ordenados
i = 0
while i<(n-1) and donde<0:
    if tramos[i]<0:
        donde = i+1
        msj.append(['sin orden en',donde])
    i = i+1

# xi tramos desiguales o no equidistantes
if (n-1)>=2: # len(tramos)>=2
    dtramos = np.diff(tramos,1) # cero o menores a casicero
    errado  = np.max(np.abs(dtramos)) # |error| mayor
    donde = -1 # suponer equidistantes
    if errado>=casicero:  # no equidistantes
        donde = np.argmax(np.abs(dtramos))+1
        msj.append(['no equidistantes',donde])


# Tabla de Diferencias Finitas
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 Finitas
diagonal = n-1 # derecha a izquierda
j = 3  # inicia en columna 3
while (j < m): # columna
    tabla_etiq.append('d'+str(j-2)+'f')
    i = 0 # fila
    while (i < diagonal): # antes de diagonal
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]

        if abs(tabla[i,j])<casicero: # casicero revisa
                tabla[i,j]=0
        i = i + 1
    diagonal = diagonal - 1
    j = j + 1

dfinita = tabla[0,3:] # diferencias finitas
h = xi[1] - xi[0]     # suponer tramos equidistantes

# polinomio con diferencias Finitas Avanzadas
x = sym.Symbol('x')
polinomio = fi[0] +0*x # en Sympy
for i in range(1,n,1):
    denominador = math.factorial(i)*(h**i)
    factor = dfinita[i-1]/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) # evaluación numérica

# SALIDA
if len(msj)>0: # mensajes de error
    print('Revisar tramos, d1x:',tramos)
    for unmsj in msj:
        print('Tramos',unmsj[0],
              'desde i:',unmsj[1])
else:
    print('Tramo h:',h)
print('Tabla Diferencia Finita')
print([tabla_etiq])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

teniendo como resultado:

Tramo h: 0.1
Tabla Diferencia Finita
[['i', 'xi', 'fi', 'd1f', 'd2f', 'd3f', 'd4f']]
[[ 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 convierte a la forma lambda para usar con Numpy:

>>> px(0.1)
1.4500000000000004
>>> px(0.2)
1.6000000000000025
>>> px0.3)
1.7000000000000042
>>>

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

5. Gráfica de polinomio de interpolación

Las instrucciones son semejantes a las presentadas en polinomio de interpolación. Se añade instrucciones para el caso en que los puntos no sean equidistantes en el eje de las x.

Como el bloque de gráfica se usa en otros métodos de la unidad, se incorpora la verificación de lista de errores 'msj' usando el bloque 'try-except'.

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

titulo = 'Interpolación: Diferencias Finitas Avanzadas'
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)')

try: # existen mensajes de error
    msj_existe = len(msj)
except NameError:
    msj = []
    
if len(msj)>0: # 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()

Tarea: se recomienda realizar las gráficas comparativas entre métodos, debe mostrar la diferencia con los métodos que requieren el tamaño de paso equidistante h, y los que no lo requieren. Permite detectar errores de selección de método para interpolación.

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]
..

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. Las funciones pueden ser reutilizadas de ser necesarias en el siguiente método de interpolación.

Tabla de Diferencias finitas
[['i', 'xi', 'fi', 'd1f', 'd2f', 'd3f', 'd4f']]
[[ 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.  ]]
xi_ascendente: True
tramos xi Equidistantes: True
Tramo h: 0.1
polinomio: 
1.30000000000000 +
1.5*x +
-2.50000000000001*(x - 0.2)*(x - 0.1) +
41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1)

polinomio simplificado: 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333335*x + 0.999999999999999

Instrucciones en Python

# Diferencias finitas Avanzadas - Interpolación
# funciones dif_finitas
# Tarea: Verificar tamaño de vectores
import numpy as np
import math
import sympy as sym

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)
    # Calular 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 ---------------------

# INGRESO
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

tipo_tabla = 'finitas'

# 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)   
tabla,titulo = diferencias_tabla(xi,fi,tipo=tipo_tabla,
                                 vertabla=True, casicero = casicero)

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 con el mismo bloque de instrucciones de la sección gráfica.

[ Dif_Finitas avanzadas ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ gráfica ]
[ función ]