Interpolar – Mascota 50 años

Ejercicio de interpolación

  • Desarrollar por grupos
  • integrar soluciones en un mismo gráfico
  • mostrar polinomios y rangos de existencia

# caparazon superior
xi = [113., 117, 134, 153, 169, 184, 194, 203]
fi = [127., 141, 161, 166, 160, 155, 140, 132]
puntos = [[xi,fi]]

# caparazon inferior 1
xi = [113., 123, 130, 149, 182, 197, 208, 211]
fi = [127., 116, 112, 107, 110, 114, 112, 108]
puntos.append([xi,fi])

# caparazon inferior 2
xi = [107., 114, 120, 143, 170, 192, 210]
fi = [124., 114, 108, 99, 99, 106, 107]
puntos.append([xi,fi])

# Patas 01
xi = [110., 116, 120, 121, 134, 143, 153, 168, 173, 177, 181, 188, 194, 201, 205, 210, 214, 218]
fi = [92., 86, 89, 86, 87, 95, 91, 93, 85, 85, 83, 86, 88, 86, 91, 87, 92, 92]
puntos.append([xi,fi])

# Patas 02
xi = [109., 115, 117, 120, 125, 130, 134, 138, 142, 145, 150, 152, 154]
fi = [91., 87, 82, 82, 78, 78, 81, 78, 76, 79, 77, 80, 86]
puntos.append([xi,fi])

# Patas 03
xi = [172., 175, 182, 186, 191, 195, 201, 205, 210, 211, 217]
fi = [86., 79, 79, 77, 82, 78, 81, 79, 82, 85, 90]
puntos.append([xi,fi])

# Patas 04
xi = [113., 118, 122, 127, 132, 136, 140, 144, 149, 152]
fi = [88., 90, 88, 86, 85, 83, 87, 84, 86, 85]
puntos.append([xi,fi])

# Rabito 01
xi = [97., 102, 108, 111, 117, 121]
fi = [120., 113, 108, 105, 102, 102]
puntos.append([xi,fi])

# cabeza 01
xi = [194., 196, 203, 210, 211, 209]
fi = [177., 167, 157, 149, 138, 135]
puntos.append([xi,fi])

# cabeza 02
xi = [195., 199, 207, 214, 219, 224, 229, 234, 239, 242, 244]
fi = [177., 185, 190, 190, 188, 193, 192, 185, 175, 172, 168]
puntos.append([xi,fi])

# cabeza 03
xi = [220., 226, 234, 239, 241]
fi = [164., 162, 163, 163, 163]
puntos.append([xi,fi])

# cabeza 04
xi = [203., 211, 214, 219, 223, 224, 225, 230, 236, 241]
fi = [115., 119, 125, 124, 127, 137, 146, 149, 154, 162]
puntos.append([xi,fi])

# cabeza 05
xi = [208., 212, 215, 219, 221, 225, 228]
fi = [174., 177, 178, 179, 182, 183, 178]
puntos.append([xi,fi])

# cabeza 06
xi = [206., 210, 214, 218, 220, 224, 229, 233, 232]
fi = [179., 182, 182, 181, 174, 171, 170, 175, 180]
puntos.append([xi,fi])

Interpolar – Mascota descansando

Referencia: Burden 9th Ed. Ejercicio 32 p164

La parte superior de una mascota descansando se describe con los puntos presentados en la tabla.  Los puntos se usan para ejercicios de interpolación.

xi = [1,2,5,6,7,8,10,13,17,20,23,24,25,27,27.7,28,29,30]
fi = [3.0,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5,7.0,6.1,5.6,5.8,5.2,4.1,4.3,4.1,3.0]

4.4.1 Trazadores Cúbicos Natural o libre

Referencia:  Burden 9Ed 3.5 p147 pdf.151, Chapra 18.6.3 p532 pdf556, Rodriguez 6.11.1 p244

Tiene como objetivo incorporar condiciones adicionales para la función en los extremos del intervalo donde se encuentran los puntos o «nodos».

Por ejemplo, si los datos son los de una trayectoria en un experimento de física, se podría disponer de la aceleración el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.

El método busca obtener un polinómio de tercer grado para cada sub-intervalo o tramo entre «nodos» consecutivos (j,  j+1), de la forma:

S_j(x_j) = a_j + b_j (x-x_j) + c_j(x-xj)^2 + d_j(x-x_j)^3

para cada j = 0, 1, …, n-1

Para n datos existen n-1 tramos, cuatro incógnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1) incógnitas para todo el intervalo .

Para los términos (xj+1xj) de los tramos que se usan varias veces en el desarrollo, se simplifican como hj:

h_j = x_{j+1} - x_{j} S_j(x_j) = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3

Para generar el sistema de ecuaciones, se siguen los siguientes criteros:

1.  Los valores de la función deben ser iguales en los nodos interiores

S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1})

2. Las primeras derivadas en los nodos interiores deben ser iguales

S'_j(x_{j+1}) = S'_{j+1}(x_{j+1})

3. Las segundas derivadas en los nodos interiores deben ser iguales

S''_j(x_{j+1}) = S''_{j+1}(x_{j+1})

4. El primer y último polinomio deben pasar a través de los puntos extremos

f(x_0) = S_0(x_0) = a_0 f(x_n) = S_n(x_n) = a_n

5. Una de las condiciones de frontera se satisface:

5a. frontera libre o natural: Las segundas derivadas en los nodos extremos son cero

S''(x_0) = S''(x_n) = 0

5b. frontera sujeta: las primeras derivadas en los nodos extremos son conocidas

S'(x_0) = f'(x_0)

S'(x_n) = f'(x_n)

Tarea: Revisar en los textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.

Ubicar en el texto las ecuaciones resultantes, que son las que se aplicarán en el algoritmo.


Algoritmo en Python

El algoritmo parte de lo desarrollado para «trazadores lineales», donde se presentaron los bloques para:

  • construir el trazador en una tabla de polinomios por tramos,
  • graficar los trazadores en cada tramo, requientdo,
  • evaluar cada uno de ellos en cada tramo con muestras suficientes para presentar una buena precisión en la gráfica.

Del algoritmo básico se modifica entonces el bloque del cálculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores cúbicos naturales (enviado a revisar como tarea).

# Trazador cúbico natural
# Condición: S''(x_0) = S''(x_n) = 0
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def traza3natural(xi,yi):
    n = len(xi)
    # Valores h
    h = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        h[j] = xi[j+1] - xi[j]
    
    # Sistema de ecuaciones
    A = np.zeros(shape=(n-2,n-2), dtype = float)
    B = np.zeros(n-2, dtype = float)
    S = np.zeros(n, dtype = float)

    A[0,0] = 2*(h[0]+h[1])
    A[0,1] = h[1]
    B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])

    for i in range(1,n-3,1):
        A[i,i-1] = h[i]
        A[i,i] = 2*(h[i]+h[i+1])
        A[i,i+1] = h[i+1]
        factor21 = (yi[i+2]-yi[i+1])/h[i+1]
        factor10 = (yi[i+1]-yi[i])/h[i]
        B[i] = 6*(factor21 - factor10)
        
    A[n-3,n-4] = h[n-3]
    A[n-3,n-3] = 2*(h[n-3]+h[n-2])
    factor12 = (yi[n-1]-yi[n-2])/h[n-2]
    factor23 = (yi[n-2]-yi[n-3])/h[n-3]
    B[n-3] = 6*(factor12 - factor23)
    
    # Resolver sistema de ecuaciones S
    r = np.linalg.solve(A,B)
    for j in range(1,n-1,1):
        S[j] = r[j-1]
    S[0] = 0
    S[n-1] = 0
    
    # Coeficientes
    a = np.zeros(n-1, dtype = float)
    b = np.zeros(n-1, dtype = float)
    c = np.zeros(n-1, dtype = float)
    d = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        a[j] = (S[j+1]-S[j])/(6*h[j])
        b[j] = S[j]/2
        factor10 = (yi[j+1]-yi[j])/h[j]
        c[j] = factor10 - (2*h[j]*S[j]+h[j]*S[j+1])/6
        d[j] = yi[j]
    
    # Polinomio trazador
    x = sym.Symbol('x')
    px_tabla = []
    for j in range(0,n-1,1):

        pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
        pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
        
        pxtramo = pxtramo.expand()
        px_tabla.append(pxtramo)
    
    return(px_tabla)

# PROGRAMA de prueba
# INGRESO , Datos de prueba
xi = np.array([0.1 , 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.8, 1.7, 2.0])
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = traza3natural(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print(' x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

con lo que se obtiene el resultado por cada tramo

Polinomios por tramos: 
 x = [0.1,0.2]
-146.666666666667*x**3 + 44.0*x**2 + 0.566666666666666*x + 1.1
 x = [0.2,0.3]
283.333333333333*x**3 - 214.0*x**2 + 52.1666666666667*x - 2.34
 x = [0.3,0.4]
-136.666666666667*x**3 + 164.0*x**2 - 61.2333333333333*x + 9.0
>>>

Para el trazado de la gráfica se añade al final del algoritmo:

# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)
    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)
    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'ro', label='puntos')
plt.plot(xtraza,ytraza, label='trazador'
         , color='blue')
plt.title('Trazadores Cúbicos Naturales')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

Si los polinomios no se igualan entre los nodos, tampoco sus velocidades y aceleraciones, podría considerar la experiencian semejante a variaciones de velocidad y aceleración en una trayectoria:

4.4 Trazadores lineales (Splines) grado1

Referencia: Chapra 18.6.1 p525 pdf549

El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (spline) para dibujar curvas suaves a través de un conjunto de puntos.

La unión más simple entre dos puntos es una línea recta. El método crea un polinomio para cada par de puntos consecutivos en el intervalo, por lo que el resultado será una tabla de polinomios.

Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales.

f(x) = f(x_0) + m_0(x-x_0), x_0\leq x\leq x_1 f(x) = f(x_1) + m_1(x-x_1), x_1\leq x\leq x_2

f(x) = f(x_{n-1}) + m_{n-1}(x-x_{n-1}), x_{n-1}\leq x\leq x_n

donde

m_i = \frac{f(x_{i+1}) - f(x_i)}{(x_{i+1}-x_i)}

Observe que la expresión de f(x) para un tramo entre dos puntos es el polinomio de grado 1 realizado con diferencia finita avanzadas  o las diferencias divididas.

Las ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y xn. Al localizar primero el intervalo dentro del cual está el punto, puede seleccionar el polinomio que corresponde a ese tramo.

Polinomios por tramos: 
 x = [0.1,0.2]    p(x) = 3.5*x + 1.1
 x = [0.2,0.3]    p(x) = -1.0*x + 2.0
 x = [0.3,0.4]    p(x) = 3.0*x + 0.8

Algoritmo en Python

Datos de los puntos como ejemplo para el algoritmo

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

El método con trazadores lineales, permite plantear los bloques necesarios para manejar varios polinomios, uno por cada tramo entre dos puntos dentro del intervalo del problema

El método permitirá disponer de un punto de partida para trazadores de mayor grado, por ejemplo los cúbicos.

Los polinomios de cada tramo se almacenan en una tabla, cada uno puede ser utilizado individualmente en su respectivo tramo, por ejemplo para realizar la gráfica de línea entre tramos.

# Trazador (spline) lineal, grado 1
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def trazalineal(xi,fi):
    n = len(xi)
    x = sym.Symbol('x')
    px_tabla = []
    tramo = 1
    while not(tramo>=n):
        # con 1ra diferencia finita avanzada 
        numerador = fi[tramo]-fi[tramo-1]
        denominador = xi[tramo]-xi[tramo-1]
        m = numerador/denominador
        pxtramo = fi[tramo-1] + m*(x-xi[tramo-1])
        px_tabla.append(pxtramo)
        tramo = tramo + 1
    return(px_tabla)

# PROGRAMA
# INGRESO , Datos de prueba
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = trazalineal(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print('  x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

Se obtiene como resultado:

Polinomios por tramos: 
  x = [0.1,0.2]
3.5*x + 1.1
  x = [0.2,0.3]
-1.0*x + 2.0
  x = [0.3,0.4]
3.0*x + 0.8
>>> 

Para añadir la gráfica se añaden las instrucciones para:

  • evaluar el polinomio en cada tramo
  • concatenar los resultados de todos los tramos en los vectores xtraza, ytraza.
  • poner en la gráfica los puntos del problema y las líneas que genera cada polinomio
# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)
    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)
    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'o', label='puntos')
plt.plot(xtraza,ytraza, label='trazador')
plt.title('Trazadores lineales (splines)')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

4.5 Interpolación paramétrica

Referencia: Rodriguez 6.9.2 p236, Burden 9Ed 3.6 p164

En algunos casos, los datos (x,y) no tienen una relación de tipo funcional y(x), entonces no se pueden aplicar directamente los métodos de interpolación revisados.

Por ejemplo, en la trayectoria del balón en el «gol imposible», la gráfica de la trayectoria en el espacio o sus proyecciones en los planos dependen del parámetro tiempo «t»en lugar de una relación de x,y,z

Referencia: 1Eva_IT2018_T4 El gol imposible

Sin embargo si las coordenadas (x,y) se expresan como funciones de otra variable t denomiada parámetro, entonces los puntos x(t), y(t) tienen relación funcional, y se pueden construir polinomios de interpolación.


Ejemplo

Las coordinadas x(t) y y(t) del recorrido de un cohete registradas en los instantes t fueron:

ti  = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

Usaremos un algoritmo en Python para mostrar la trayectoria x,y para el problema planteado.

Al realizar la interpolación de los puntos para obtener polinomios que dependen de «t» se obtiene:

px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

polinomios con los que se puede realizar la gráfica px(t), py(t) en forma separada. Pero para comprender mejor la trayectoria del cohete, se utiliza la gráfica px(t) vs py(t) en el intervalo t entre[0,3]

Las intrucciones para mostrar el resultado son:

# interpolación paramétrica
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
ti = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

# PROCEDIMIENTO
# interpolando con lagrange
px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

t = np.arange(0,3,0.01)
puntosx = px(t)
puntosy = py(t)

# Salida
plt.plot(puntosx,puntosy)
plt.show()

4.2.2 Diferencias divididas de Newton

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

4.2.1 Diferencias finitas avanzadas polinomio

Referencia: Rodriguez 6.6.4 Pdf.221, Burden 9Ed p129

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

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 k en el intervalo [x0,xn]

\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]

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


Ejercicio

Se toman los datos del ejercicio de «diferencias finitas» y se grafican los resultados. Se resalta en la gráfica la distancia equidistante h entre los puntos consecutivos xi como los mostrados en la gráfica anterior.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

y se usan los resulados previos de diferencias finitas:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 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.  ]]

valores con los que se inicia el cálculo de la distancia constante entre puntos xi:

h = xi[1] – xi[0] = 0.2-0.1 = 0.1

Los valores de diferencias finitas a usar en la fórmula son de la primera fila desde la columna 3 en adelante.

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

La construción de la expresión del polinomio empieza con el valor del primer punto de la función f(0.1)

p_3(x) = f_0 = 1.45

para añadirle el primer término j=1, hay que calcular el factor:

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

con lo que el término se completa 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-02) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-02)

Finalmente, añadiendo el tercer 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-02) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Se puede seguir simpificando 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-02)+ + 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.


Algoritmo en Python

El polinomio se contruye usando el ejercicio de diferencias finitas.

Para contruir 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.

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

# Polinomio interpolación
# Diferencias finitas avanzadas
# 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([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)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
# diferencias finitas 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('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = 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):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    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
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
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('Interpolación polinómica')
plt.show()

teniendo como resultado:

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 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 convirtió a lambda:

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

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

4.2 Diferencias finitas

Referencia: Rodriguez 6.5 Pdf.211, Chapra 18.1.3 pdf.509

Establece relaciones simples entre los puntos dados que describen la funcion f(x). Las tabla de diferencias finitas es un elemento base para varios métodos de interpolación, por lo que se trata como un tema inicial.

Para un ejemplo se toman los siguientes puntos:

xi 0.10 0.2 0.3 0.4
fi 1.45 1.6 1.7 2.0

La tabla de diferencias finitas se construye tomando los datos y sus índices como parte de las primeras columnas.

i xi fi Δ1fi Δ2fi Δ3fi Δ4fi
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. 0.

Cada casila de diferencia finita Δjfi se obtiene restando los dos valores consecutivos de la columna anterior. Por ejemplo, para la primera columna:

\Delta ^1 f_0 = 1.6-1.45 = 0.15 \Delta ^1 f_1 = 1.7-1.6 = 0.1 \Delta ^1 f_2 = 2.0-1.7 = 0.3

y la siguiente Δ1f3 no es posible calcular.

Siguiendo el mismo procedimiento se calculan los valores de la siguiente columna Δ2fi como las diferencias de la columna anterior, así sucesivamente.

Si la función f(x) de donde provienen los datos es un polinonio de grado n, entonces la n-ésima diferencia finita será una constante, y las siguientes diferencias se anularán.


Algoritmo en Python

Para crear la tabla de diferencias finitas, las primeras columnas requieren concatenar los valores de los índice ixi y fi.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

Los índices i se crean en un vector ki, pues la variable i es usada como fila en matrices, por lo que evitamos confundirlas al usar la variable.

A la matriz con las tres columnas descritas, se le añade a la derecha una matriz de nxn para calcular las diferencias.

Se calculan las diferencias para cada columna, realizando la operación entre los valores de cada fila. Considere que no se realizan cálculos desde la diagonal hacia abajo en la tabla, los valores quedan como cero.

Al final se muestra el título y el resultado de la tabla.

# Tabla de Diferencias finitas
# resultado en: [título,tabla]
# Tarea: verificar tamaño de vectores
import numpy as np

# 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)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
# diferencias finitas 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('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# SALIDA
print('Tabla Diferencia Finita: ')
print([titulo])
print(tabla)

el resultado de aplicar el algoritmo es:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 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.  ]]

>>> 

Gráfica de puntos con líneas horizontales

Para tener una referencia visual sobre las primeras diferencias finitas, en una gráfica se trazan las líneas horizontales que pasan por cada punto. Para las segundas diferencias se debe graficar las primeras diferencias finitas vs xi repitiendo el proceso. Las líneas de distancia se añadieron con un editor de imágenes.

las intrucciones adicionales al algoritmo anterior para añadir la gráfica son:

# Gráfica
import matplotlib.pyplot as plt
for i in range(0,n,1):
    plt.axhline(fi[i],ls='--', color='yellow')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencia Finita')
plt.show()

Interpolar – Pato en pleno vuelo

Referencia: Burden Cap 3. Ejemplo 1.

La figura muestra un joven pato en pleno vuelo. Para aproximar el perfil de la parte superior del pato, se presentan 21 puntos a lo largo de la curva de aproximación relativos a un sistema de coordenadas sobrepuestas.

xi = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fi = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]


a) Realice la interpolación entre puntos usando la interpolación de Lagrange para la parte superior del pato en vuelo.

b) ¿Es posible crear un solo polinomio para los 21 puntos presentados? Explique su respuesta en pocas palabras.

c) ¿En concepto, cuando considera necesario crear un nuevo polinomio?

d) Adjunte en los archivos para: la gráfica que interpola los puntos, y los polinomios que la generan.

e) (Puntos extra: 10, para promediar con lección o taller)

Se adjuntan los demás puntos del perfil del pato en pleno vuelo, presente el o los polinomios que interpolan la figura (gráfica y polinomios), las despuestas a literales b,c,d deben también ser válidas para el ejercicio completo.

# Perfil Superior del pato
xiA = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fiA = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]

# Perfil inferior cabeza
xiB = [0.817, 0.897, 1.022, 1.191, 1.510, 1.834, 2.264, 2.962, 3.624, 4.202, 4.499, 4.779, 5.109, 5.527]
fiB = [1.180, 1.065, 1.023, 1.010, 1.032, 1.085, 1.192, 1.115, 1.087, 1.100, 0.830, 0.608, 0.350, 0.106]

# Perfil Ala superior
xiC = [4.659, 4.865, 5.085, 5.261, 5.387, 5.478, 5.527]
fiC = [-5.161, -4.741, -3.933, -2.951, -1.970, -0.981, 0.106]

# Perfil Ala inferior
xiD = [4.659, 4.750, 4.990, 5.289, 5.560, 5.839, 6.113, 6.606, 6.916, 7.305, 7.563, 7.802, 7.983]
fiD = [-5.161, -5.259, -5.284, -5.268, -5.161, -4.982, -4.769, -4.286, -3.911, -3.213, -2.670, -2.176, -1.655]

# Perfil inferior posterior
xiE = [8.141, 8.473, 8.832, 9.337, 9.887, 10.572, 10.995, 11.501, 11.923, 12.364, 12.763, 13.300]
fiE = [-1.138, -0.434, -0.514, -0.494, -0.382, -0.005, -0.090, -0.085, -0.030, 0.093, 0.120, 0.250]

# Perfil ojo superior
xiF = [2.663, 2.700, 2.805, 2.886]
fiF = [2.202, 2.279, 2.293, 2.222]

# Perfil ojo inferior
xiG = [2.663, 2.720, 2.826, 2.886]
fiG = [2.202, 2.130, 2.143, 2.222]