数据分析（Python） 插值——牛顿插值法
插值思想
求已知的n个点对(x1,y1),(x2,y2),…,(xn,yn)的所有阶差商公式
联立以上差商公式建立如下插值多项式f(x)
其中P(x)是牛顿插值逼近函数，R(x)是误差函数
差商表
| xk | f(xk) | 一阶差商 | 二阶差商 | 三阶差商 |
|---|---|---|---|---|
| x0 | f(x0) | |||
| x1 | f(x1) | f[x0,x] | ||
| x2 | f(x2) | f[x1,x2] | f[x0,x1,x2] | |
| x3 | f(x3) | f[x2,x3] | f[x1,x2,x3] | f[x0,x1,x2,x3] |
| … | … | … | … | … |
代码实现
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()
编写思想
牛顿插值代码编写难度主要在于差商公式，可以设计一个二维数组，观察每行列数据与其余元素对应关系，即可写出差商公式。
结果展示
