通常我们在求插值节点的开头部分插值点附近函数值时,使用Newton前插公式;求插值节点的末尾部分插值点附近函数值时,使用Newton后插公式。
代码:
1 # -*- coding: utf-8 -*-
2 """ 3 Created on Wed Mar 25 15:43:42 2020
4
5 @author: 35035
6 """ 7
8
9 import numpy as np 10
11 # 等距节点的Newton向前插值(输入的x向量和y向量注意保证从小到大顺序) 12 def Newton_iplt_forword(x, y, xi): 13 """x,y是插值节点,xi是一个值"""
14 n = len(x) 15 m = len(y) 16 if n != m: 17 print("Error!") 18 return None 19 h = x[1] - x[0] 20 t = (xi - x[0]) / h 21
22 # 先计算差分表(cf) 23 cf = [] 24 temp = y.copy() 25 for i in range(n): 26 if i != 0: 27 iv_1 = temp[i - 1] 28 for j in range(i, n): 29 iv_2 = temp[j] 30 temp[j] = iv_2 - iv_1 31 iv_1 = iv_2 32 cf.append(temp[i]) 33 # 再计算Newton插值 34 ans = 0
35 for i in range(n): 36 w = 1
37 # 计算多项式部分,让差分作为其系数 38 for j in range(i): 39 w *= ((t - j) / (j + 1)) 40 ans += w*cf[i] 41 return ans 42
43 # 等距节点的Newton向后插值(输入的x向量和y向量保证从小到大顺序) 44 def Newton_iplt_backword(x, y, xi): 45 """x,y是插值节点ndarray,xi是一个值"""
46 n = len(x) 47 m = len(y) 48 if n != m: 49 print("Error!") 50 return None 51 h = x[1] - x[0] 52 t = (xi - x[n - 1]) / h 53
54 # 先计算差分表(cf) 55 cf = [] 56 temp = y.copy() 57 for i in range(n): 58 if i != 0: 59 iv_1 = temp[i - 1] 60 for j in range(i, n): 61 iv_2 = temp[j] 62 temp[j] = iv_2 - iv_1 63 iv_1 = iv_2 64 cf.append(temp[n - 1]) 65 # 再计算Newton插值 66 ans = 0
67 for i in range(n): 68 w = 1
69 # 计算多项式部分,让差分作为其系数 70 for j in range(i): 71 w *= ((t + j) / (j + 1)) 72 ans += w*cf[i] 73 return ans 74
75 # 当对多个值使用Newton插值时,使用map()建立映射: 76 # Iterator = map(Newton, Iterable) 77
78 # 数值运算时使用float参与运算,dtype定为内置float 79
80 x = np.array((1,2,3,4,5,6), dtype=float) 81 y = np.array((1.0, 1.2599, 1.4422, 1.5874, 1.71, 1.8171), dtype=float) 82 print(Newton_iplt_forword(x, y, 5.6)) 83 print(Newton_iplt_backword(x, y, 5.6)) 84 # 结果:1.775416 测试成功!