牛顿差商插值

牛顿插值python实现

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)

结果：