牛顿插值python实现

import  numpy as np #数值运算
import sympy
import matplotlib.pyplot as plt
import pandas as pdclass 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)# 类型转换，数据结构采用arrayif 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_quotdef 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替换多项式中的tself.y0 = y_0return  y_0def 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)

结果：