4.1 El polinomio de interpolación con Sympy-Python

Referencia: Rodriguez 6.1 p191, Chapra 18.3 p520 pdf544

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

Por ejemplo, dados los 4 puntos en la tabla 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
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

Es posible generar un sistema de ecuaciones en base p(x) que pase por cada uno de los puntos. Por ejemplo si se toma el primer punto con xi=0 y fi=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 ai. 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.

Los métodos requieren usar la matriz A y vector B:

\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, se construye observando que los coeficientes se elevan al exponente referenciado al índice columna pero de derecha a izquierda. La última columa da 1 por tener como coeficiente el valor de xi0.

Por los métodos para resolver sistemas de ecuaciones, se observa que hay que realizar la matriz aumentada, pivotear por filas, etc.

Para la solución se propone un algoritmo realizado en Python, que para enfocarnos en la interpolación, se obtiene el siguente resultado:

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

con lo que se contruye el polinomio con los valores ai.

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

Polinomio que se grafica dentro del intervalo de los datos del problema y se obtiene la línea que pasa por los puntos.


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 xi 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 para algebra lineal, que entrega el vector solución que representan los coeficientes del polinomio de interpolación.

Para contruir 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 facilitar la evación del polinomio de forma numérica, se convierte el polinomio a la forma Lambda. Con un número de muestras suficientes para suavizar la curva dentro del intervalo [a,b] del problema se puede realizar la gráfica del polinomio.

Finalmente, con la gráfica se comprueba que la curva pase por todos los puntos del intervalo.

# El polinomio de interpolación
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

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

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi)
B = np.array(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

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
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])
    
# 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)

# Grafica
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()

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

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.