5.2 El polinomio de interpolación

Referencia: Rodriguez 6.1 pdf.191

Unir n puntos en el plano usando un polinomios de interpolación se puede realizar con un polinomio de grado (n-1) puntos de datos 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) = a0x3 + a1x2 + a2x1 + a3

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 los puntos, luego se resuelve con alguno de los métodos conocidos.

a0(0.0)3 + a1(0.0)2 + a2(0.0)1 + a3 = 1
a0(0.2)3 + a1(0.2)2 + a2(0.2)1 + a3 = 1.6
a0(0.3)3 + a1(0.3)2 + a2(0.3)1 + a3 = 1.7
a0(0.4)3 + a1(0.4)2 + a2(0.4)1 + a3 = 2.0

Para la solución se propone un algoritmo realizado en python, teniendo como resultado:

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

# Interpolación por polinomio
import numpy as np
import matplotlib.pyplot as plt

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

# PROCEDIMIENTO
# Arreglos numpy 
xi = np.array(xi)
fi = np.array(fi)
n = len(xi)

# Vector B en columna
B = np.array(fi)
B = fi[:,np.newaxis]

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for f in range(0,n,1):
    for c in range(0,n,1):
        potencia = (n-1)-c
        D[f,c] = xi[f]**potencia

# Resolver matriz aumentada
coeficiente =  np.linalg.solve(D,B)

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

# Convierte polinomio a funcion para evaluación más rápida
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio
k = 100
inicio = np.min(xi)
fin = np.max(xi)
puntosx = np.linspace(inicio,fin,k)
puntosy = px(puntosx)

# Usando evaluación simbólica
##puntosy = np.zeros(k,dtype=float)
##for j in range(0,k,1):
##    puntosy[j] = polinomio.subs(x,puntosx[j])
    
# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)

# Grafica
plt.plot(xi,fi,'o')
plt.plot(puntosx,puntosy)
plt.show()