使用Python求牛顿插值多项式及其差商表（附加一个拉格朗日插值多项式）

闲话不多说，直接上代码。

import numpy as np
from sympy import *

# 定义一个求差商表的函数，使用递归求解差商表，返回值是差商的值
# x是数组，表示样本点的x
# f是数组，表示样本点的函数值f（x）
# start是int类型，表示x数组的起始下标，
# end是int类型，表示x数组的结束下标，
# res是数组类型，表示差商表，对角线以下为各阶差商表
def cs(x,f,start,end,res):
    # 当x中只有两个元素的时候，结束递归
    if((end-start)==1):
        # 将一阶差商填入差商表
        res[end-1][end-start-1]=(f[end]-f[start])/(x[end]-x[start])
        # 返回差商
        return res[end-1][end-start-1]
    # 当x中元素个数大于2时，根据公式递归调用此函数求差商，并将差商填入差商表
    res[end-1][end-start-1]=(cs(x,f,start+1,end,res)-cs(x,f,start,end-1,res))/(x[end]-x[start])
    # 返回差商
    return res[end-1][end-start-1]

# 定义一个求牛顿插值多项式的函数
# x是数组，表示样本点的x
# f是数组，表示样本点的函数值f（x）
def Newton(x,f):
    res = np.ones([x.size - 1, x.size - 1])*np.inf # 初始化差商表骨架结构
    cs(x, f, 0, x.size - 1, res) #调用差商表函数给差商表填值，对角线及以下才是差商表的值
    X=symbols("x") #定义x变量
    y=f[0] #初始化牛顿插值多项式，它的第一项是常数项，正好是f[0]
    for i in range(x.size-1):
        temp=1 #临时变量，保存  f[x0,x1,...,xn]*(x-x1)(x-x2)...(x-xn-1)
        for j in range(i+1):
            temp=temp*(X-x[j]) #(x-x1)(x-x2)...(x-xn-1)
        temp=res[i,i]*temp #将差商表对角线元素作为系数
        y=y+temp #将temp添加到多项式中
    # 返回牛顿插值多项式
    return y

if __name__=="__main__":
    # 样本点
    x=np.array([2.0,2.1,2.2,2.3,2.4])
    f=np.array([1.414214,1.449138,1.483240,1.516575,1.549193])

    ##### 可以直接得到差商表
    res = np.ones([x.size - 1, x.size - 1]) * np.inf  # 初始化差商表骨架结构
    # 调用差商表函数
    cs(x,f,0,x.size-1,res)
    print(res)

    #### 也可以直接得到牛顿插值多项式
    X=symbols("x") #定义自变量x
    y=Newton(x,f) #调用函数得到牛顿插值多项式，类型是
    print("N(x)=",y)
    # 给自变量x赋值，求出函数近似值
    print(y.evalf(subs={X:2.15})) #求给定自变量x值时函数f(x)的值 | 将表达式转化为函数

 得到的差商表：

 牛顿插值多项式（比较长，就截取了部分）：

拉格朗日插值多项式代码（使用方法很简单，和牛顿插值多项式一样）： 

#拉格朗日插值法
def L(x,f):
    X=symbols("x")
    m=x.size #x个数
    L=0
    for i in range(m):
        temp=1
        for j in range(m):
            if(i!=j):
                temp=temp*((X-x[j])/(x[i]-x[j]))
        L=L+temp*f[i]
    return L

 各位大哥点个赞呐（卑微）