1. n阶差商实现
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
|
def diff(xi,yi,n):
"""
param xi:插值节点xi
param yi:插值节点yi
param n: 求几阶差商
return: n阶差商
"""
if len(xi) != len(yi): #xi和yi必须保证长度一致
return
else:
diff_quot = [[] for i in range(n)]
for j in range(1,n+1):
if j == 1:
for i in range(n+1-j):
diff_quot[j-1].append((yi[i]-yi[i+1]) / (xi[i] - xi[i + 1]))
else:
for i in range(n+1-j):
diff_quot[j-1].append((diff_quot[j-2][i]-diff_quot[j-2][i+1]) / (xi[i] - xi[i + j]))
return diff_quot
|
测试一下：
|
1
2
3
4
|
xi = [1.615,1.634,1.702,1.828]
yi = [2.41450,2.46259,2.65271,3.03035]
n = 3
print(diff(xi,yi,n))
|
返回的差商结果为：
[[2.53105263157897, 2.7958823529411716, 2.997142857142854], [3.0440197857724347, 1.0374252793901158], [-9.420631485362996]]
2. 牛顿插值实现
|
1
2
3
4
5
6
7
8
9
|
def Newton(x):
f = yi[0]
v = []
r = 1
for i in range(n):
r *= (x - xi[i])
v.append(r)
f += diff_quot[i][0] * v[i]
return f
|
测试一下：
|
1
2
|
x = 1.682
print(Newton(x))
|
结果为：
2.5944760289639732
3.完整Python代码
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
|
def Newton(xi,yi,n,x):
"""
param xi:插值节点xi
param yi:插值节点yi
param n: 求几阶差商
param x: 代求近似值
return: n阶差商
"""
if len(xi) != len(yi): #xi和yi必须保证长度一致
return
else:
diff_quot = [[] for i in range(n)]
for j in range(1,n+1):
if j == 1:
for i in range(n+1-j):
diff_quot[j-1].append((yi[i]-yi[i+1]) / (xi[i] - xi[i + 1]))
else:
for i in range(n+1-j):
diff_quot[j-1].append((diff_quot[j-2][i]-diff_quot[j-2][i+1]) / (xi[i] - xi[i + j]))
print(diff_quot)
f = yi[0]
v = []
r = 1
for i in range(n):
r *= (x - xi[i])
v.append(r)
f += diff_quot[i][0] * v[i]
return f
|
到此这篇关于用Python实现牛顿插值法的文章就介绍到这了,更多相关python牛顿插值法内容请搜索服务器之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持服务器之家！
原文链接：https://blog.csdn.net/weixin_45812669/article/details/115694394