使用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) #調用函數得到牛頓插值多項式，類型是<class 'sympy.core.add.Add'>
    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

 各位大哥點個贊吶（卑微）