插值思想

求已知的n个点对(x1,y1),(x2,y2),…,(xn,yn)的所有阶差商公式

联立以上差商公式建立如下插值多项式f(x)

其中P(x)是牛顿插值逼近函数，R(x)是误差函数

差商表

代码实现

import numpy as np
import matplotlib.pyplot as plt
import random

y_point = []
text_list = []
x_point = random.sample(range(0,12),10)#生成不重复的随机数
for num in range(10):
    y_point.append(random.randint(1, 12))#生成1-20随机整数
    text='('+str(x_point[num])+','+str(y_point[num])+')'#合并字符串
    text_list.append(text)
    plt.annotate(text_list[num],xy=(x_point[num],y_point[num]),xytext=(x_point[num]+0.5,y_point[num]+0.5))
    #annotate第一个参数是文本内容 第二个参数是要标记的位置 第三个参数是文本标记位置

x = np.arange(min(x_point)-0.01,max(x_point)+0.003,0.003)

def whole(n):#计算数值(x-x0)(x-x1)...(x-xn) n为项数次序
    if(n==1):
        return y_point[0]
    else:
        sum = 1
        for i in range(n-1):
            sum = sum * (x-x_point[i])
        return sum

def difference(n):#计算差分公式 f(x0,x1) n为项数次序
    if n==1:
        return 1
    else:
        _2Dlist = [([0] * (n - 1)) for i in range(n - 1)]#从f(x0,x1)开始的三角矩阵 #f(x0,x1),f(x1,x2),f(x2,x3)
        for i in range(n - 1):
            _2Dlist[i][0] = (y_point[i + 1] - y_point[i]) / (x_point[i + 1] - x_point[i]) #先计算三角矩阵第一列
        for a in range(1,n-1): #a代表列 n=4   a=2,b=23;a=3,b=3
            for b in range(a,n-1):  #b代表行
                _2Dlist[b][a]=(_2Dlist[b][a-1]-_2Dlist[b-1][a-1])/(x_point[b+1]-x_point[b-a]) #计算三角矩阵除第一列以外所有数据元素
    print(_2Dlist) #打印差商表
    return _2Dlist[n-2][n-2]


def newton(n):#n为数据元素个数
    sum = 0
    for i in range(n):
        sum = sum + whole(i+1)*difference(i+1)
    return sum

plt.scatter(x_point,y_point)
plt.plot(x,newton(10))
plt.show()

编写思想

牛顿插值代码编写难度主要在于差商公式，可以设计一个二维数组，观察每行列数据与其余元素对应关系，即可写出差商公式。

结果展示

作者：Atom_QQ2022313691

链接：https://atom2022313691.blog.csdn.net/article/details/107002299

来源：CSDN
著作权归作者所有，转载请联系作者获得授权，切勿私自转载。