1. Interpolación polinómica de Lagrange
Referencia: Chapra 18.2 p516, Burden 3.1 p80, Rodríguez 6.2 p195
El polinomio de interpolación de Lagrange reformula el polinomio de interpolación de Newton evitando el cálculo de la tabla de diferencias divididas. El método de Lagrange tolera las diferencias entre distancias x de los puntos de muestra.
El polinomio de Lagrange se construye a partir de las fórmulas:
f_{n} (x) = \sum_{i=0}^{n} L_{i} (x) f(x_{i}) L_{i} (x) = \prod_{j=0, j \neq i}^{n} \frac{x-x_j}{x_i - x_j}Donde una vez que se han seleccionado los puntos a usar que generan la misma cantidad de términos que puntos.
2. Ejercicio
Dados los 4 puntos en la tabla se requiere un polinomio de grado 3.
| xi | 0 | 0.2 | 0.3 | 0.4 |
|---|---|---|---|---|
| fi | 1 | 1.6 | 1.7 | 2.0 |
3. Desarrollo analítico
Se calculan cuatro términos de Lagrange según las fórmulas,
término 1
L_{0} (x) = \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)}término 2
L_{1} (x) = \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)}término 3
L_{2} (x) = \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)}término 4
L_{3} (x) = \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}se construye el polinomio usando la fórmula para fn(x) para cada valor fi,
p_3(x) = 1 L_{0} (x) + 1.6 L_{1} (x) + 1.7 L_{2} (x) + 2 L_{3} (x) p_3(x) = 1 \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)} + + 1.6 \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)} + 1.7 \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)} + 2 \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}La simplificación de la expresión del polinomio se puede dejar como tarea o realizarla con Sympy con la instrucción sym.expand().
Una forma de verificar que el polinomio es correcto es usar un punto original xi[i] y comprobar que la evaluación tiene como resultado fi[i].
3. Algoritmo en Python
En el algoritmo en Python, para construir las expresiones de cada término de Lagrange se usa la forma simbólica con Sympy, con los datos:
xi = [0, 0.2, 0.3, 0.4]
fi = [1, 1.6, 1.7, 2.0]En cada término Li(x) se usan todos los elementos i , excepto el mismo elemento i, en el numerador y denominador de la expresión.
En polinomio se agrupan todos los términos multiplicados por su respectivo valor fi[i].
valores de fi: [1. 1.6 1.7 2. ]
divisores en L[i]: [-0.024 0.004 -0.003 0.008]
Polinomio de Lagrange, expresiones
250.0*x*(x - 0.3)*(x - 0.2) +
400.0*x*(x - 0.4)*(x - 0.3) +
-566.666666666667*x*(x - 0.4)*(x - 0.2) +
-41.6666666666667*(x - 0.4)*(x - 0.3)*(x - 0.2)
Polinomio de Lagrange:
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333334*x + 1.0
3 2
41.6666666666667⋅x - 27.5⋅x + 6.83333333333334⋅x + 1.0Las expresiones del polinomio contiene los binomios en su forma básica, para resolver y simplificar las ecuaciones se usa polinomio.expand().
Para realizar la gráfica del polinomio es conveniente convertirlo a la forma lambda con Numpy, de esta forma se evalúa en una sola linea todos los puntos para el intervalo [a,b] en x.
# Interpolacion de Lagrange
# divisoresL solo para mostrar valores denominador
import numpy as np
import sympy as sym
# INGRESO , Datos de prueba
xi = [0, 0.2, 0.3, 0.4]
fi = [1, 1.6, 1.7, 2.0]
# PROCEDIMIENTO
# Vectores como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Polinomio de Lagrange
x = sym.Symbol('x')
polinomio = 0*x # sym.S.Zero en Sympy
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
# Termino de Lagrange
numerador = 1
denominador = 1
for j in range(0,n,1):
if (j!=i):
numerador = numerador*(x-xi[j])
denominador = denominador*(xi[i]-xi[j])
terminoLi = numerador/denominador
polinomio = polinomio + terminoLi*fi[i]
divisorL[i] = denominador
polisimple = polinomio.expand() # simplifica los (x-xi)
px = sym.lambdify(x,polisimple) # evaluación numérica
# SALIDA
print(' valores de fi: ',fi)
print('divisores en L[i]: ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
#print(polinomio)
terminos = sym.Add.make_args(polinomio)
n_term = len(terminos)
for i in range(0,n_term,1):
if i<(n_term-1):
print(terminos[i],'+')
else:
print(terminos[i])
print()
print('Polinomio de Lagrange: ')
print(polisimple)
print()
sym.pprint(polisimple)
4. Gráfica del polinomio
Para la gráfica se añaden las siguientes instrucciones
# Gráfica --------------
import matplotlib.pyplot as plt
titulo = 'Interpolación Lagrange'
muestras = 21 # resolución gráfica
a = np.min(xi) # intervalo [a,b]
b = np.max(xi)
xk = np.linspace(a,b,muestras)
yk = px(xk)
plt.plot(xi,fi,'o', label = '[xi,fi]')
plt.plot(xk,yk, label = 'p(x)')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title(titulo)
plt.show()
5. Algoritmo como función
Se convierte el algoritmo a una función para usar en ejercicios donde se requiera mas de un proceso de interpolación:
Se expresa el resultado en mas detalle en las operaciones,
Interpola con Lagrange
1.0 * (x - 0.4)*(x - 0.3)*(x - 0.2) /( (0.0 - 0.2)*(0.0 - 0.3)*(0.0 - 0.4) )
+ 1.6 * x*(x - 0.4)*(x - 0.3) /( (-0.0 + 0.2)*(0.2 - 0.3)*(0.2 - 0.4) )
+ 1.7 * x*(x - 0.4)*(x - 0.2) /( (-0.0 + 0.3)*(-0.2 + 0.3)*(0.3 - 0.4) )
+ 2.0 * x*(x - 0.3)*(x - 0.2) /( (-0.0 + 0.4)*(-0.2 + 0.4)*(-0.3 + 0.4) )
polinomio simplificado
41.666667*x**3 - 27.5*x**2 + 6.833333*x + 1
Polinomio de Lagrange:
41.666667*x**3 - 27.5*x**2 + 6.833333*x + 1
3 2
41.666667⋅x - 27.5⋅x + 6.833333⋅x + 1con lo que las instrucciones quedan como:
# Interpolacion de Lagrange
# divisoresL solo para mostrar valores denominador
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]
# Algoritmos como funcion
def interpola_Lagrange(xi,fi,vertabla=False,
precision=6, casicero = 1e-15):
'''
Interpolación con método de Lagrange
resultado: polinomio en forma simbólica
'''
# Matrices como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
if vertabla==True:
print('Interpola con Lagrange')
# Polinomio con Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = sym.S.Zero
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
termino = 1 # Termino de Lagrange
numerador = 1
denominador = 1
for j in range(0,n,1):
if (j!=i):
numerador = numerador*(x-xi[j])
denominador = denominador*(sym.UnevaluatedExpr(xi[i])-sym.UnevaluatedExpr(xi[j]))
terminoLi = numerador/denominador
if vertabla==True:
txt0='' ; txt1='('; txt2=')'
if i>0:
txt0='+'
if n>2:
txt1=''; txt2=''
print(txt0,fi[i],'*'+txt1,numerador,txt2+'/('+ txt1,
denominador,txt2,')')
polinomio = polinomio + (terminoLi.doit())*fi[i]
denominador = denominador.doit()
divisorL[i] = denominador
# simplifica el polinomio
polisimple = polinomio.expand()
try: # intenta redondear coeficientes a precisiion
polisimple = redondea_coef(polisimple, precision,casicero)
except NameError:
print("redondea_coef() no se encuentra definida.")
if vertabla==True:
print('polinomio simplificado')
print(polisimple)
return(polisimple)
def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
''' redondea coeficientes de términos suma de una ecuacion
'''
tipo = type(ecuacion)
tipo_eq = False
if tipo == sym.core.relational.Equality:
RHS = ecuacion.rhs
ecuacion = ecuacion.lhs
tipo = type(ecuacion)
tipo_eq = True
if tipo == sym.core.add.Add: # términos suma de ecuacion
term_sum = sym.Add.make_args(ecuacion)
ecuacion = sym.S.Zero
for term_k in term_sum:
# factor multiplicativo de termino suma
term_mul = sym.Mul.make_args(term_k)
producto = sym.S.One
for factor in term_mul:
if not(factor.has(sym.Symbol)):
factor = np.around(float(factor),precision)
if (abs(factor)%1)<casicero: # si es entero
factor = int(factor)
producto = producto*factor
ecuacion = ecuacion + producto
if tipo == sym.core.mul.Mul: # termino único, busca factores
term_mul = sym.Mul.make_args(ecuacion)
producto = sym.S.One
for factor in term_mul:
if not(factor.has(sym.Symbol)):
factor = np.around(float(factor),precision)
if (abs(factor)%1)<casicero: # si es entero
factor = int(factor)
producto = producto*factor
ecuacion = producto
if tipo == float: # si es entero
if (abs(ecuacion)%1)<casicero:
ecuacion = int(ecuacion)
if tipo_eq:
ecuacion = sym.Eq(ecuacion,RHS)
return(ecuacion)
# PROGRAMA DE PRUEBA -----------------------------
# PROCEDIMIENTO
polisimple = interpola_Lagrange(xi,fi,vertabla=True)
x = sym.Symbol('x')
px = sym.lambdify(x,polisimple)
# SALIDA
print('Polinomio de Lagrange: ')
print(polisimple)
print()
sym.pprint(polisimple)
La gráfica se puede realizar con las instrucciones anteriores.
6.Función en librería Scipy.interpolate.lagrange
>>> 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. ])
>>>