数值计算之牛顿插值

import numpy as np

def diff_quot(xi, yi, n):
    """
    通过定义计算差商
    输出为差商表主对角线, type为list
    如输入[1,2,3]和[1,3,8]，输出[2, 1.5]
    """
    quot=[]
    temp=yi
    for i in range(1, n+1):
        # 计算i阶差商
        tem = np.diff(temp) / [xi[_+i] - xi[_] for _ in range(n-i+1)]
        # 分母为i-1阶差商的逐项差值，分子为xi间隔i项的差值，/为向量逐元素相除(numpy特性)
        quot.append(tem[0]) # 将i阶差商的第一个元素添加到结果
        temp=tem
    return(quot)

def Newton(xi, yi, n, x):
    Dt = diff_quot(xi, yi, n)
    y = yi[0]
    for i in range(n):
        omiga = 1 # 计算连乘
        for j in range(i+1):
            omiga *= (x - xi[j])
        y += Dt[i] * omiga # 相加
    return y

xi = np.array([0.4, 0.55, 0.65, 0.80, 0.90])
yi = np.array([0.41075, 0.57815, 0.69675, 0.88811, 1.02652])
n = 4

print("差商计算结果：", diff_quot(xi, yi, n))
print("牛顿差值结果：", Newton(xi, yi, n, 0.596))