import numpy as np
import sympy # 符号运算库
class LagrangeInterpolation:
'''
拉格朗日插值
'''
def __init__(self,x,y):
'''
拉格朗日必要参数的初始化,及各健壮性的检测;健壮性:系统在不正常输入情况下仍能表现正常的程度
:param x: 已知数据x的坐标点
:param y: 已知数据y的坐标点
'''
#构造一下作为类的属性,防止后续进行数组的运算,做一个类型转换
self.x = np.asarray(x, dtype=np.float)
self.y = np.asarray(y, dtype=np.float) #类型转换,数据结构采用array
if len(self.x)>1 and len(self.x)==len(self.y): #如果x,y是一个数据点,就不用做了,且个数应相同
self.n=len(self.x) #有多少个数据点,即已知离散数据点的个数
else:
raise ValueError("插值数据(x,y)维度不匹配")
#插值最终结果是一个多项式的形式,可作为一个类的属性由用户进行返回
self.polynomial = None #最终的插值多项式,符号表示
self.poly_coefficient = None # 最终插值多项式的系数向量,幂次从高到低
self.coefficient_order = None # 对应多项式系数的阶次(所有项中最高次的幂)
self.y0 = None # 所求插值点的值,单个值或向量
def fit_interp(self): # 拟合
'''
核心算法:生成拉格朗日插值多项式
:return:
'''
t = sympy.Symbol("t") # 定义符号变量,因为参数已经用过x
self.polynomial = 0.0 # 插值多项式实例化
for i in range(self.n): # 由多少个样本点,n个数据生成n-1次
# 根据每个数据带点,构造插值基函数
basis_fun = self.y[i] # 插值基函数
# 计算每个 yi*li(x)
for j in range(i):
basis_fun *= (t-self.x[j]) / (self.x[i]-self.x[j])
for j in range(i + 1 , self.n):
basis_fun *= (t - self.x[j]) / (self.x[i] - self.x[j])
# 下面开始累加
self.polynomial +=basis_fun # 插值多项式累加
self.polynomial = sympy.expand(self.polynomial) #多项式的展开
print(self.polynomial)
def cal_interp_x0(self, x0): # 计算插值的值
'''
计算所给定的插值点的值,即插值
:param x0: 所求插值的x坐标
:return:
'''
pass
测试
from lagelangri.lagrange_interpolation import LagrangeInterpolation
import numpy as np
if __name__ =='__main__':
x = np.linspace(0, 2*np.pi, 10 ,endpoint=True) #0到2*pi,等分十个,最后一个值包含
y=np.sin(x)
"""print(x)
print(y)"""
lag_interp = LagrangeInterpolation(x=x,y=y)
lag_interp.fit_interp()
未简化的结果
0.920725428958529*t*(1.125 - 0.179049310978382*t)*(1.14285714285714 - 0.204627783975294*t)*(1.16666666666667 - 0.238732414637843*t)*(1.2 - 0.286478897565412*t)*(1.25 - 0.358098621956765*t)*(1.33333333333333 - 0.477464829275686*t)*(1.5 - 0.716197243913529*t)*(2.0 - 1.43239448782706*t) + 0.705316598492019*t*(1.28571428571429 - 0.204627783975294*t)*(1.33333333333333 - 0.238732414637843*t)*(1.4 - 0.286478897565412*t)*(1.5 - 0.358098621956765*t)*(1.66666666666667 - 0.477464829275686*t)*(2.0 - 0.716197243913529*t)*(3.0 - 1.43239448782706*t)*(1.43239448782706*t - 1.0) + 0.413496671566344*t*(1.5 - 0.238732414637843*t)*(1.6 - 0.286478897565412*t)*(1.75 - 0.358098621956765*t)*(2.0 - 0.477464829275686*t)*(2.5 - 0.716197243913529*t)*(4.0 - 1.43239448782706*t)*(0.716197243913529*t - 0.5)*(1.43239448782706*t - 2.0) + 0.122476942006377*t*(1.8 - 0.286478897565412*t)*(2.0 - 0.358098621956765*t)*(2.33333333333333 - 0.477464829275686*t)*(3.0 - 0.716197243913529*t)*(5.0 - 1.43239448782706*t)*(0.477464829275686*t - 0.333333333333333)*(0.716197243913529*t - 1.0)*(1.43239448782706*t - 3.0) - 0.0979815536051016*t*(2.25 - 0.358098621956765*t)*(2.66666666666667 - 0.477464829275686*t)*(3.5 - 0.716197243913529*t)*(6.0 - 1.43239448782706*t)*(0.358098621956765*t - 0.25)*(0.477464829275686*t - 0.666666666666667)*(0.716197243913529*t - 1.5)*(1.43239448782706*t - 4.0) - 0.206748335783172*t*(3.0 - 0.477464829275686*t)*(4.0 - 0.716197243913529*t)*(7.0 - 1.43239448782706*t)*(0.286478897565412*t - 0.2)*(0.358098621956765*t - 0.5)*(0.477464829275686*t - 1.0)*(0.716197243913529*t - 2.0)*(1.43239448782706*t - 5.0) - 0.201519028140577*t*(4.5 - 0.716197243913529*t)*(8.0 - 1.43239448782706*t)*(0.238732414637843*t - 0.166666666666667)*(0.286478897565412*t - 0.4)*(0.358098621956765*t - 0.75)*(0.477464829275686*t - 1.33333333333333)*(0.716197243913529*t - 2.5)*(1.43239448782706*t - 6.0) - 0.115090678619816*t*(9.0 - 1.43239448782706*t)*(0.204627783975294*t - 0.142857142857143)*(0.238732414637843*t - 0.333333333333333)*(0.286478897565412*t - 0.6)*(0.358098621956765*t - 1.0)*(0.477464829275686*t - 1.66666666666667)*(0.716197243913529*t - 3.0)*(1.43239448782706*t - 7.0) - 3.89817183251938e-17*t*(0.179049310978382*t - 0.125)*(0.204627783975294*t - 0.285714285714286)*(0.238732414637843*t - 0.5)*(0.286478897565412*t - 0.8)*(0.358098621956765*t - 1.25)*(0.477464829275686*t - 2.0)*(0.716197243913529*t - 3.5)*(1.43239448782706*t - 8.0)
self.polynomial = sympy.expand(self.polynomial)多项式展开后
-2.29369601281319e-6*t**9 + 6.48527268907878e-5*t**8 - 0.00061946678451133*t**7 + 0.00167479859899278*t**6 + 0.0040407632285091*t**5 + 0.00707566335693954*t**4 - 0.173843738744177*t**3 + 0.00401585715826069*t**2 + 0.999072579744674*t