import numpy as np #数值运算
import sympy
import matplotlib.pyplot as plt
import pandas as pd
class NewtonDifferenceQuotient:
"""
牛顿差商插值算法
"""
def __init__(self,x,y):
"""
牛顿差商插值必要参数的初始化以及健壮性的检测
:param x: 已知数据x坐标点
:param y: yi知数据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):
self.n = len(self.x) #已知离散数据点的个数
else:
raise ValueError("插值数据（x,y）维度不匹配！")
self.polynomial = None #最终的插值多项式，符号表示
self.poly_coefficient = None#最终插值多项式的系数向量，幂次从高到底
self.coefficient_order =None#对应多项式系数的阶次
self.y0 =None# 所求插值点的值，单个值或向量
self.diff_quot =None #存储离散数据点的差商
def __diff_quotient__(self):
"""
计算牛顿差商表
:return:
"""
diff_quot = np.zeros((self.n,self.n)) #存储差商
diff_quot[:,0] = self.y #差商表中第一列存储 y值
for j in range(1,self.n):#按列计算
for i in range(j,self.n): #行 对角线
diff_quot[i,j] = (diff_quot[i,j-1]-diff_quot[i-1,j-1])
/(self.x[i]-self.x[i-j])
self.diff_quot = pd.Dataframe(diff_quot)
return diff_quot
def fit_interp(self):
"""
核心算法：生成牛顿差商插值多项式
:return:
"""
diff_quot = self.__diff_quotient__()#计算差商
d_q = np.diag(diff_quot) #构造牛顿差商插值只需要对角线的值即可
t =sympy.Symbol("t")#定义符号变量
self.polynomial=d_q[0]#插值多项式实例化
term_poly = t-self.x[0]# 初始牛顿每项
for i in range(1,self.n):
#针对每个数据点构造
self.polynomial+=d_q[i]*term_poly;
term_poly*= (t-self.x[i])
#插值多项式特征
self.polynomial=sympy.expand(self.polynomial)#多项式展开
polynomial = sympy.Poly(self.polynomial,t)#根据多项式构造多项式对象
self.poly_coefficient = polynomial.coeffs() #获取多项式系数
self.coefficient_order=polynomial.monoms() #获取多项式系数对应阶次
# print(self.polynomial)
def cal_interp_x0(self,x0):
"""
计算所给定的插值点的数值，即插值
:param x0: 所求插值的x坐标
:return:
"""
x0 = np.asarray(x0,dtype=np.float)
n0 = len(x0)# 插值点个数
y_0=np.zeros(n0) #存储插值点x0所对应的插值
t = self.polynomial.free_symbols.pop() #返回值是集合，获取插值多项式的自由变量
for i in range(n0):
y_0[i] = self.polynomial.evalf(subs={t:x0[i]})# 用x0替换多项式中的t
self.y0 = y_0
return y_0
def plt_interpolation(self,x0=None,y0=None):
"""
可视化插值图像和所求的插值点
:return:
"""
plt.figure(figsize=(8,6))
plt.plot(self.x,self.y,"ro",label = "Interpolation base points")
xi = np.linspace(min(self.x),max(self.x),100) #模拟100个值
yi = self.cal_interp_x0(xi)
plt.plot(xi,yi,"b--",label = "Interpolation polynomisl")
if x0 is not None and y0 is not None:
plt.plot(x0, y0, "g*", label="Interpolation point values")
plt.legend()
plt.xlabel("x",fontdict={"fontsize":12})
plt.ylabel("y", fontdict={"fontsize": 12})
plt.title("Lagrange interpolation polynomial and values",fontdict={"fontsize":14})
plt.grid(ls=":")
plt.show()
主函数
if __name__ == "__main__":
x = np.linspace(0,2*np.pi,10,endpoint=True)
y = np.sin(x)
x0 = np.array([np.pi/2,2.158,3.58,4.784])
ndq = NewtonDifferenceQuotient(x=x,y=y)
ndq.fit_interp()
ndq.__diff_quotient__()
print("牛顿插值差商表：")
print(ndq.diff_quot)
print("牛顿差商插值多项式如下：")
print(ndq.polynomial)
print("牛顿差商插值多项式系数和对应阶次：")
print(ndq.poly_coefficient)
print(ndq.coefficient_order)
y0 = ndq.cal_interp_x0(x0)
print("所求插值点的值是:",y0,"精确值为：",np.sin(x0))
print(ndq.y0)
ndq.plt_interpolation(x0,y0)
结果：
