插值算法 - 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  测试成功！