计算方法：牛顿插值及其代码实现

Posted on Edited on In , Views:

这篇博文记录计算方法中牛顿插值法及其Python实现，并且记录差商的相关知识。

牛顿插值法的引入

拉格朗日插值可以解决一般的插值问题，但是问题是，每次添加一个新的插值点，就需要重新计算一次插值基函数，这会浪费大量的算力、降低效率，我们希望找到一种方法，能够在加入新的插值点时，不重新计算原先的部分，这就使用了牛顿插值法。

差商

定义

给定函数在互不相同点处的函数值，定义为在处的一阶差商

接着我们定义为在处的二阶差商

一般地，如果我们有阶差商，那么我们可以定义阶差商为

性质

证明可参考csdn文章

牛顿插值公式

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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import numpy as np
import sympy

class NewtonInterpolation:
    def __init__(self, x, y):
        self.x = np.asarray(x, dtype=float)
        self.y = np.asarray(y, dtype=float)
        self.n = len(x)
        self.diff_mat = np.zeros((self.n, self.n))
        self.polynomial = None  # 最终的插值多项式
        self.biuld_diff_mat()
        self.build_polynomial()

    def biuld_diff_mat(self):
        for i in range(self.n):
            self.diff_mat[i][0] = self.y[i]
        for i in range(1,self.n):
            for j in range(i,self.n):
                self.diff_mat[j][i] = (self.diff_mat[j][i-1]-self.diff_mat[j-1][i-1])/(self.x[j]-self.x[j-i])
    
    def build_polynomial(self):
        t = sympy.Symbol("t")
        self.polynomial = self.diff_mat[0][0]
        for i in range(1,self.n):
            temp = self.diff_mat[i][i]
            for j in range(i):
                temp *= (t-self.x[j])
            self.polynomial += temp
        self.polynomial = sympy.expand(self.polynomial)
    
    def calc_interp(self, x0):
        if self.polynomial == None:
            return ValueError("未生成插值多项式")
        
        x0 = np.asarray(x0, dtype=float)
        n0 = len(x0)
        y_0 = np.zeros(n0)
        t = self.polynomial.free_symbols.pop()
        for i in range(n0):
            y_0[i] = self.polynomial.evalf(subs={t: x0[i]})
        return y_0

if __name__ == "__main__":
    # x = [0.4, 0.55, 0.65, 0.8, 0.9]#, 1.05
    # y = [0.41075, 0.57815, 0.69675, 0.88811, 1.02652]#, 1.25382
    x = [1,2,3,4,5]
    y = [1,2,3,4,5]
    newton = NewtonInterpolation(x, y)
    # print(newton.diff_mat)
    print(newton.polynomial)
    print(newton.calc_interp([0.596]))