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. Por ejemplo si se toma el primer pungo con xi =0 y fi=1 se establece la ecuación:

a0(0.0)3 + a1(0.0)2 + a2(0.0)1 + a3 = 1

Luego se continúa con los otros puntos seleccionados hasta completar las ecuaciones necesarias para el grado del polinomio seleccionado.

El sistema obtenido 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()