Categoría: Unidad 4 Interpolación polinómica

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] [ Dif_Divididas_Newton ]
..

---

1. Diferencias finitas

Referencia: Rodríguez 6.5 p211, Chapra 18.1.3 p509

Establece relaciones simples entre los puntos dados que describen la función 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 casilla de diferencia finita Δ^jf_i 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 Δ¹f₃ no es posible calcular.

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

\Delta ^2 f_0 = 0.1-0.15 = -0.05 \Delta ^2 f_1 = 0.3-0.1 = 0.2 \Delta ^3 f_0 = 0.2-(-0.05) = 0.25

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

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] [ Dif_Divididas_Newton ]

---

2. Algoritmo en Python

Para crear la tabla de diferencias finitas, las primeras columnas requieren concatenar los valores de los índice i,  xi 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.  ]]

>>> 

---

3. 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 instrucciones 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()

[ Dif_Finitas ] [ Dif_Finitas Avanzadas] [ Dif_Divididas_Newton ]

Interpolación: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

1. Ejercicio

Referencia: Rodríguez 6.1 p191, Chapra 18.3 p520

Unir n puntos en el plano usando un polinomios de interpolación se puede realizar con un polinomio de grado (n-1).

xi	0	0.2	0.3	0.4
fi	1	1.6	1.7	2.0

Por ejemplo, dados los 4 puntos en la tabla [xi,fi] se requiere generar un polinomio de grado 3 de la forma:

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3

Interpolación: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

2. Desarrollo Analítico

Es posible generar un sistema de ecuaciones para p(x) haciendo que pase por cada uno de los puntos o coordenadas.
Por ejemplo si se toma el primer punto con x_i=0 y f_i=1 se establece la ecuación:

p(0) = 1 a_0 (0)^3 + a_1 (0)^2 + a_2 (0)^1 + a_3 = 1

Note que ahora las incógnitas son los coeficientes a_i. Luego se continúa con los otros puntos seleccionados hasta completar las ecuaciones necesarias para el grado del polinomio seleccionado.

a_0 (0.0)^3 + a_1 (0.0)^2 + a_2 (0.0)^1 + a_3 = 1.0 a_0 (0.2)^3 + a_1 (0.2)^2 + a_2 (0.2)^1 + a_3 = 1.6 a_0 (0.3)^3 + a_1 (0.3)^2 + a_2 (0.3)^1 + a_3 = 1.7 a_0 (0.4)^3 + a_1 (0.4)^2 + a_2 (0.4)^1 + a_3 = 2.0

El sistema obtenido se resuelve con alguno de los métodos conocidos para la Solución a sistemas de ecuaciones, que requieren escribir las ecuaciones en la forma de matriz A y vector B, desarrollar la matriz aumentada, pivotear por filas, etc.

\begin{pmatrix} 0.0^3 & 0.0^2 & 0.0^1 & 1 \\ 0.2^3 & 0.2^2 & 0.2^1 & 1 \\ 0.3^3 & 0.3^2 & 0.3^1 & 1 \\ 0.4^3 & 0.4^2 & 0.4^1 & 1 \end{pmatrix} . \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} = \begin{pmatrix} 1.0 \\ 1.6 \\ 1.7 \\ 2.0 \end{pmatrix}

La matriz A también se conoce como Matriz Vandermonde D, que se construye observando que los coeficientes se elevan al exponente referenciado al índice columna pero de derecha a izquierda. La última columna tiene valores 1 por tener como coeficiente el valor de xi⁰.

Para enfocarnos en la interpolación, en la solución se propone usar un algoritmo o función en Python, obteniendo el siguiente resultado:

los coeficientes del polinomio: 
[ 41.66666667 -27.5          6.83333333   1.        ]

a partir del cual que se construye el polinomio con los valores obtenidos.

Polinomio de interpolación p(x): 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

que re-escrito es,

p(x) = 41.6666 x^3 - 27.5 x^2 + 6.8333 x + 1

Se realiza  la gráfica de p(x) dentro del intervalo de los datos del ejercicio y se obtiene la línea contínua que pasa por cada uno de los puntos.

Interpolación: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

3. Algoritmo en Python

Para desarrollar el ejercicio, se realiza un bloque para construir la matriz A y vector B, usando los vectores que representan los puntos de muestra de la función o experimento.

xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]

Observe que en éste caso los datos se dan en forma de lista y se deben convertir hacia arreglos.

La matriz A o Matriz Vandermonde D, se construye observando que los coeficientes del vector x_i se elevan al exponente referenciado al índice columna pero de derecha a izquierda.

Construir el polinomio consiste en resolver el sistema de ecuaciones indicado en la sección anterior. Para simplificar la solución del sistema, se usa Numpy, que entrega el vector solución que representan los coeficientes del polinomio de interpolación.

Para construir la expresión del polinomio, se usa la forma simbólica con Sympy, de forma semejante a la usada para construir el polinomio de Taylor.
Para mostrar el polinomio de una manera más fácil de interpretar se usa la instrucción sym.pprint(), usada al final del algoritmo.

# El polinomio de interpolación
import numpy as np
import sympy as sym

# INGRESO
xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)

B = fi
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype=float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

con lo que se obtiene el siguiente resultado:

Matriz Vandermonde: 
[[0.    0.    0.    1.   ]
 [0.008 0.04  0.2   1.   ]
 [0.027 0.09  0.3   1.   ]
 [0.064 0.16  0.4   1.   ]]
los coeficientes del polinomio: 
[ 41.66666667 -27.5 6.83333333 1. ]
Polinomio de interpolación: 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

formato pprint
                  3         2                           
41.6666666666667*x  - 27.5*x  + 6.83333333333333*x + 1.0

Gráfica del polinomio

Para facilitar la evaluación numérica del polinomio, se convierte el polinomio a la forma Lambda px. La gráfica se realiza con un número de muestras suficientes para suavizar la curva dentro del intervalo [a,b], por ejemplo 21, con lo que se comprueba que la curva pase por todos los puntos de la tabla xi,fi dados en el ejercicio.

Instrucciones adicionales al algoritmo para la gráfica:

# GRAFICA
import matplotlib.pyplot as plt
# Polinomio a forma Lambda x:
# para evaluación con vectores de datos xin
# muestras = tramos+1
muestras = 21

px = sym.lambdify(x,polinomio)

a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])

plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()

Interpolación: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

4. Función en Scipy

>>> import scipy as sp
>>> xi = [0,0.2,0.3,0.4]
>>> fi = [1,1.6,1.7,2.0]
>>> sp.interpolate.lagrange(xi,fi)
poly1d([ 41.66666667, -27.5       ,   6.83333333,   1.        ])
>>> 

Interpolación: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]