Interpolaciones_metodos

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

#Importamos las librerías a utilizar
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO , Datos de prueba
xcoordenadas = np.array([0, 0.05, 0.1,0.1089])
funcioni = np.array([1, 1.1052, 1.2214,0.06328])
# PROCEDIMIENTO
# Polinomio de Lagrange
niteraciones = len(xcoordenadas)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(niteraciones, dtype = float)
for i in range(0,niteraciones,1):
 
 # Termino de Lagrange
 numerador = 1
 denominador = 1
 for j in range(0,niteraciones,1):
 if (j!=i):
 numerador = numerador*(x-xcoordenadas[j])
 denominador = denominador*(xcoordenadas[i]-xcoordenadas[j])
 terminoLi = numerador/denominador
 polinomio = polinomio + terminoLi*funcioni[i]
 divisorL[i] = denominador
 # simplifica el polinomio
polinomiosimple = polinomio.expand()
 # para evaluación numérica
evaluacion_num = sym.lambdify(x,polinomiosimple)
# Puntos para la gráfica
#Creamos las muestras que son necesarias para poder evaluar 
muestras = 101 #Entre cada par de puntos
punto_a = np.min(xcoordenadas)
punto_b = np.max(xcoordenadas)
punto_xi = np.linspace(punto_a,punto_b,muestras)
punto_fi = evaluacion_num(punto_xi)
#Imprimimos los resultados del polinomio de lagrange
print(' valores de funcion(i): ',funcioni)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polinomiosimple)
 # Gráfica
plt.plot(xcoordenadas,funcioni,'o', label = 'Puntos')
plt.plot(punto_xi,punto_fi, label = 'Polinomio')
plt.legend()
plt.xlabel('x coordenadas')
plt.ylabel('funcion(i)')
plt.title('Interpolación Lagrange')
plt.show()
#trazador lineal 
def trazaLineal(xi,fi):
 tramo=1
 px_tabla = []
 tramo = 1
 while not(tramo>=n):
 numerador = fi[tramo]-fi[tramo-1]
 denominador = xi[tramo]-xi[tramo-1]
 m = numerador/denominador
 pxtramo = fi[tramo-1] + m*(x-xi[tramo-1])
 px_tabla.append(pxtramo)
 tramo = tramo + 1
 print('/nInterpolacion lineal, polinomio por tramos: ')
 for tramo in range(1,n,1):
 print(' x = ['+str(xi[tramo-1])
 +','+str(xi[tramo])+']')
 print(str(px_tabla[tramo-1])) 
 
 return(px_tabla)
n=len(xcoordenadas)
px_tabla=trazaLineal(xcoordenadas,funcioni)
# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
 a = xcoordenadas[tramo-1]
 b = xcoordenadas[tramo]
 xtramo = np.linspace(a,b,muestras)
 # evalua polinomio del tramo
 pxtramo = px_tabla[tramo-1]
 pxt = sym.lambdify('x',pxtramo)
 ytramo = pxt(xtramo)
 # vectores de trazador en x,y
 xtraza = np.concatenate((xtraza,xtramo))
 ytraza = np.concatenate((ytraza,ytramo))
 tramo = tramo + 1
# Gráfica
plt.plot(xtraza,ytraza,label='trazador lineal')
plt.plot(xcoordenadas,funcioni,'o', label='puntos')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('trazador lineal')
plt.show()
#Colocación con funciones de base radial RBF
# Función de evaluación de FBR multicuádrica
def rbffunction(xev, xdat, c):
 rbfv = np.sqrt((xev - xdat)**2 + c**2)
 return rbfv
# Construcción de la matriz de interpolación
def interpmat(xdat, c):
 nd = len(xdat)
 mat1 = np.zeros((nd,nd),float)
 for i in range (nd):
 for j in range (nd):
 mat1[i,j] = rbffunction(xdat[i], xdat[j], c)
 return mat1
# Superposición de funciones de base radial
def rbfsuperposit(x, coef, xdat, c):
 y = np.zeros((len(x)))
 for i in range (len(x)):
 for j in range (len(xdat)):
 y[i] = y[i] + coef[j]*rbffunction(x[i], xdat[j], c)
 return y
# Información de entrada
xdat = xcoordenadas
ydat = funcioni
c = 0.5 # Parámetro de forma
# Matriz de interpolación
matint = interpmat(xdat, c)
# Coeficientes de la interpolación
coef = np.linalg.solve(matint, ydat)
# Evaluación de la superposición sobre un intervalo
xespacio = np.linspace(0,0.15,100)
yinterp = rbfsuperposit(xespacio, coef, xdat, c)
fbr1 = rbffunction(xespacio, xdat[0], c)
fbr2 = rbffunction(xespacio, xdat[1], c)
fbr3 = rbffunction(xespacio, xdat[2], c)
fbr4 = rbffunction(xespacio, xdat[3], c)
# Cálculo del error RMS entre la interpolación y la función dada
Err = np.sqrt((np.sum((yinterp - (np.cos(xespacio))**10))**2)/len(yinterp))
print('Parámetro de forma: ', c)
print('Error RMS de la aproximación: ', Err)
# Sumas parciales de la interpolación
yinterp6 = rbfsuperposit(xespacio, coef[0:6], xdat[0:6], c)
yinterp7 = rbfsuperposit(xespacio, coef[0:7], xdat[0:7], c)
yinterp8 = rbfsuperposit(xespacio, coef[0:8], xdat[0:8], c)
# Gráficas
plt.figure()
plt.plot(xdat, ydat, 'or', label = 'Datos')
plt.plot(xespacio, yinterp, label = 'Interpolación RBF')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('Interpolación con Funciones de Base Radial')
plt.show()
#Trazadores cuadraticos
def trazadorcuadratico(xi,yi):
 n = len(xi)
 a = np.zeros(n-1, dtype = float)
 b = np.zeros(n-1, dtype = float)
 c = np.zeros(n-1, dtype = float)
 for i in range(n-1):
 h = xi[i+1]-xi[i]
 a[i] = (yi[i+1]-yi[i])/h
 b[i] = yi[i]-a[i]*xi[i]
 c[i] = 0
 
 px_tabla = []
 for i in range(n-1):
 pxtramo = a[i]*(x-xi[i])**2 + b[i]*(x-xi[i]) + c[i]
 pxtramo = pxtramo.expand()
 px_tabla.append(pxtramo)
 
 return px_tabla
n = len(xcoordenadas)
px_tabla = trazadorcuadratico(xcoordenadas,funcioni)
print("Trazadores cuadraticos")
print('Polinomios por tramos: ')
for tram in range(1,n,1):
 print(' x = ['+str(xcoordenadas[tram-1])
 +','+str(xcoordenadas[tram])+']')
 print(str(px_tabla[tram-1]))
#Después graficamos los resultados
# Puntos para graficar cada tramo
xtraz = np.array([])
ytraz = np.array([])
tram = 1
while not(tram>=n):
 a = xcoordenadas[tram-1]
 b = xcoordenadas[tram]
 xtramo = np.linspace(a,b,muestras)
 
 # evalua polinomio del tramo
 pxtramo = px_tabla[tram-1]
 pxt = sym.lambdify('x',pxtramo)
 ytramo = pxt(xtramo)
 # vectores de trazador en x,y
 xtraz = np.concatenate((xtraz,xtramo))
 ytraz = np.concatenate((ytraz,ytramo))
 tram = tram + 1
# Graficar resultado de la función trazadorcubico
plt.plot(xcoordenadas, funcioni, 'o', label='Puntos')
plt.plot(xtraz, ytraz, label='Trazadores cuadraticos')
plt.title('Trazadores cuadraticos')
plt.xlabel('xi')
plt.ylabel('yi')
plt.legend()
plt.show()
#Trazadores cúbicos
def trazadorcubico(xi,yi):
 niteraciones = len(xi)
 
 # Valores h
 hiteraciones = np.zeros(niteraciones-1, dtype = float)
 for j in range(0,niteraciones-1,1):
 hiteraciones[j] = xi[j+1] - xi[j]
 
 # Sistema de ecuaciones
 A = np.zeros(shape=(niteraciones-2,niteraciones-2), dtype = float)
 B = np.zeros(niteraciones-2, dtype = float)
 S = np.zeros(niteraciones, dtype = float)
 A[0,0] = 2*(hiteraciones[0]+hiteraciones[1])
 A[0,1] = hiteraciones[1]
 B[0] = 6*((yi[2]-yi[1])/hiteraciones[1] - (yi[1]-yi[0])/hiteraciones[0])
#Se recorre el total de iteraciones que se van a hacer en el sistema de ecuaciones
 for i in range(1,niteraciones-3,1):
 A[i,i-1] = hiteraciones[i]
 A[i,i] = 2*(hiteraciones[i]+hiteraciones[i+1])
 A[i,i+1] = hiteraciones[i+1]
 factor21 = (yi[i+2]-yi[i+1])/hiteraciones[i+1]
 factor10 = (yi[i+1]-yi[i])/hiteraciones[i]
 B[i] = 6*(factor21 - factor10)
 #Después encontramos los factores del sistema de ecuaciones
 A[niteraciones-3,niteraciones-4] = hiteraciones[niteraciones-3]
 A[niteraciones-3,niteraciones-3] = 2*(hiteraciones[niteraciones-3]+hiteraciones[niteraciones-2])
factor12 = (yi[niteraciones-1]-yi[niteraciones-2])/hiteraciones[niteraciones-2]
 factor23 = (yi[niteraciones-2]-yi[niteraciones-3])/hiteraciones[niteraciones-3]
 B[niteraciones-3] = 6*(factor12 - factor23)
 
 # Resolver sistema de ecuaciones S
 r = np.linalg.solve(A,B)
 for j in range(1,niteraciones-1,1):
 S[j] = r[j-1]
 S[0] = 0
 S[niteraciones-1] = 0
 
 # Coeficientes
 a = np.zeros(niteraciones-1, dtype = float)
 b = np.zeros(niteraciones-1, dtype = float)
 c = np.zeros(niteraciones-1, dtype = float)
 d = np.zeros(niteraciones-1, dtype = float)
 for j in range(0,niteraciones-1,1):
 a[j] = (S[j+1]-S[j])/(6*hiteraciones[j])
 b[j] = S[j]/2
 factor10 = (yi[j+1]-yi[j])/hiteraciones[j]
 c[j] = factor10 - (2*hiteraciones[j]*S[j]+hiteraciones[j]*S[j+1])/6
 d[j] = yi[j]
 
 # Polinomio trazador
 x = sym.Symbol('x')
 px_tabla = []
 for j in range(0,niteraciones-1,1):
 pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
 pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
 
 pxtramo = pxtramo.expand()
 px_tabla.append(pxtramo)
 
 return(px_tabla)
 
 
# Tabla de polinomios por tramos
n = len(xcoordenadas)
px_tabla = trazadorcubico(xcoordenadas,funcioni)
#Imprimiendo los resultados de la ecuación encontrada
print("Trazadores Cúbicos")
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
 print(' x = ['+str(xcoordenadas[tramo-1])
 +','+str(xcoordenadas[tramo])+']')
 print(str(px_tabla[tramo-1]))
#Después graficamos los resultados
# Puntos para graficar cada tramo
xtrazador = np.array([])
ytrazador = np.array([])
tramo = 1
while not(tramo>=n):
 a = xcoordenadas[tramo-1]
 b = xcoordenadas[tramo]
 xtramo = np.linspace(a,b,muestras)
 
 # evalua polinomio del tramo
 pxtramo = px_tabla[tramo-1]
 pxt = sym.lambdify('x',pxtramo)
 ytramo = pxt(xtramo)
 # vectores de trazador en x,y
 xtrazador = np.concatenate((xtrazador,xtramo))
 ytrazador = np.concatenate((ytrazador,ytramo))
 tramo = tramo + 1
# Gráfica
plt.plot(xcoordenadas,funcioni,'ro', label='puntos')
plt.plot(xtrazador,ytrazador, label='trazador cúbico'
 , color='blue')
plt.title('Trazadores Cúbicos')
plt.xlabel('xi')
plt.ylabel('puntox(xi)')
plt.legend()
plt.show()
# Gráfica
plt.plot(xcoordenadas,funcioni,'o', label = 'Puntos',color='red')
plt.plot(punto_xi,punto_fi, label = 'Polinomio De Lagrange')
plt.plot(xespacio, yinterp, label = 'Interpolación RBF')
plt.plot(xtraza,ytraza,label='trazador lineal')
plt.plot(xtraz, ytraz, label='Trazadores cuadraticos')
plt.plot(xtrazador,ytrazador, label='trazador Cúbico')
plt.legend()
plt.xlabel('x coordenadas')
plt.ylabel('funcion(i)')
plt.title('Interpolaciones')
plt.show()