数计算数学

牛顿插值的原理与实现

字数统计: 2.7k字   |   阅读时长: 12分
, , , ,

摘要: 牛顿插值的原理与应用

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

在业务场景中，我们经常要找到某个指标与某个因素之间的变化关系，这种变化关系有时可以用函数来表示。

对于这种函数，我们一般不会知道其具体形式，只能通过数据分析拿到一些散点。也就是说我们只知道函数在一些点的位置，却不知道函数的具体的表达式。此时使用代数插值的方法给出函数的近似形式是一种可以考虑的方案。

众多代数插值算法中，牛顿插值法是比较经典的一种方法，是在《自然哲学的数学原理》（1687年出版）中提出的，计算简单，逻辑清晰，适合编程实现。

本文首先看一下牛顿提出的插值过程中引入的差商的概念，然后推导出牛顿插值公式，最后我们编程实现这个算法并进行一些插值实验。

差商的定义

设函数 $f(x)$ 有 $n+1$ 个已知的点 $(x_{i}, y_{i}), i=0,1,\cdots,n$，定义 $f(x)$ 的各阶差商如下：

  • $f(x)$ 在 $x_{i}$ 的零阶差商 $f[x_{i}]$：
  • $f(x)$ 在 $x_{i}, x_{j}$ 的二阶差商 $f[x_{i}, x{j}]$：
  • $f(x)$ 在 $x_{i}, x_{j}, x_{k}$ 的三阶差商 $f[x_{i}, x{j}, x_{k}]$：

后面我们在实现牛顿插值的时候，以上定义是核心的计算逻辑。

一般地，$f(x)$ 在点 $x_{0}, x_{1}, \cdots, x_{k}$ 的 k 阶差商为：

经过推导，$f(x)$ 在点 $x_{0}, x_{1}, \cdots, x_{k}$ 的 k 阶差商 $f[x_{0}, x_{1}, \cdots, x_{k}]$ 可以用函数值 $f(x_{0}), f(x_{1}), \cdots, f(x_{k})$ 表示：

牛顿插值公式推导

按照前面的定义，当我们有了 $n+1$ 个已知点 $(x_{i}, y_{i}), i=0,1,\cdots,n$，可以按照以下步骤推导出牛顿插值公式：

首先按定义写出 $f(x)$ 的 $1 \sim n+1$ 阶差商：

分别变形，得：

依次代入，得：

记：

于是得到最终的插值公式：

其中 $R_{n}(x)$ 为牛顿插值公式的余项，也称为截断误差，当 n 趋于无穷大时为 0。

代码模板

前面我们已经推导出了牛顿插值公式：

下面我们看一下如何实现这个公式。

当我们有了 $n+1$ 个已知点 $(x_{i}, y_{i}), i=0,1,\cdots,n$，要得出 $P_{n}(x)$ 我们需要计算 $f(x)$ 在各个已知点的各阶差商。我们将需要计算的差商整理成以下表格：

$x_{i}$ $f(x_{i})$ 1 阶差商 2 阶差商 3 阶差商 n 阶差商
$x_{0}$ $f(x_{0})$
$x_{1}$ $f(x_{1})$ $f[x_{0}, x_{1}]$
$x_{2}$ $f(x_{2})$ $f[x_{1}, x_{2}]$ $f[x_{0}, x_{1}, x_{2}]$
$x_{3}$ $f(x_{3})$ $f[x_{2}, x_{3}]$ $f[x_{1}, x_{2}, x_{3}]$ $f[x_{0}, x_{1}, x_{2}, x_{3}]$
$\cdots$
$x_{n}$ $f(x_{n})$ $f[x_{n-1}, x_{n}]$ $f[x_{n-2}, x_{n-1}, x_{n}]$ $f[x_{n-3}, x_{n-2}, x_{n-1}, x_{n}]$ $f[x_{0}, x_{1},\cdots,x_{n}]$

按照插值公式，我们需要的是上表中对角线的这些。根据前面的各阶差商的定义公式，在编程实现时可以从左到右一列一列的计算。具体地，在上表中计算某个位置的差商时，按照定义，需要左边和左上这两个位置的已经算好的差商的值，类似于动态规划的过程。

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
import numpy as np
import matplotlib.pyplot as plt
import random

class NewtonInterpolation:
    def __init__(self, ps):
        """
        输入 ps 为 np.ndarray，n 个点
        ps[i] = [xi, yi] 为 第 i 个已知点，i = 0,1,...,n
        计算牛顿插值公式中的各个参数，并保存
        """
        self.n = len(ps[0]) # 已知点个数
        self.range = [min(ps[0]), max(ps[0])]
        self.ps = ps
        # self.diff_coef[i][j] := xi 的 j 阶差商
        self.diff_coef = np.zeros((self.n, self.n)) # 按照表格顺序保存差商
        self.get_diff_coef()

    def get_diff_coef(self):
        """
        计算差商
        self.diff_coef[i][j] := xi 的 j 阶差商
        """
        # 先计算三角矩阵第 0 列, 即 xi 的 0 阶差商
        for i in range(self.n):
            self.diff_coef[i][0] = self.ps[1][i]
        # 计算三角矩阵第 1 到第 n - 1 列
        for j in range(1, self.n):
            # 第 i 行
            for i in range(j, self.n):
                self.diff_coef[i][j] = (self.diff_coef[i][j - 1] - self.diff_coef[i - 1][j - 1]) / (self.ps[0][i] - self.ps[0][i - j])

    def plot(self):
        """
        画出插值后的函数图像
        """
        text_list = []
        for i in range(self.n):
            # 画 n 个已知点
            text = "({}, {})".format(self.ps[0][i], self.ps[1][i])
            text_list.append(text)
            # annotate 第一个参数是文本内容 第二个参数是要标记的位置 第三个参数是文本标记位置
            plt.annotate(text_list[i]
                        ,xy=(self.ps[0][i], self.ps[1][i])
                        ,xytext=(self.ps[0][i]+0.5, self.ps[1][i]+0.5)
                        )
        plt.scatter(self.ps[0], self.ps[1])

        # 画图的点的横坐标
        x_list = np.arange(self.range[0] - 0.01, self.range[1] + 0.003, 0.003)
        y_list = np.zeros(x_list.shape)
        for i, x in enumerate(x_list):
            y_list[i] = self.__call__(x)
        plt.plot(x_list, y_list)
        plt.show()

    def __call__(self, x):
        """
        对于插值后的函数，输入 x，返回 y
        """
        # 需要遍历 self.diff_coef[i][i]，同时乘以一个系数 tmp
        # 计算 tmp = (x-x0)(x-x1)...
        tmp = [1] * self.n
        for i in range(1, self.n):
            tmp[i] = tmp[i - 1] * (x - self.ps[0][i - 1])
        y = 0
        for i in range(self.n):
            y += tmp[i] * self.diff_coef[i][i]
        return y

def main():
    # 生成已知点的横坐标，不重复
    x_point = random.sample(range(0, 12), 10)
    # 生成已知点的纵坐标，并在图上标出
    y_point = random.sample(range(0, 12), 10)
    newton_interpolation = NewtonInterpolation([x_point, y_point])
    newton_interpolation.plot()

if __name__ == "__main__":
    main()

代码中 NewtonInterpolation 在创建实例时接收已知点，__init__ 的时候根据已知点计算前面的差商表格，具体计算逻辑在 get_diff_coef 方法中，动态规划算法，时间复杂度 $O(n^{2})$。

初始化完成后，牛顿插值公式所需的各个系数分别存在 self.ps 和 self.diff_coef 中。

__call__ 方法接收一个新的 x，根据插值后的函数（背后是已知点 ps 和差商表格 diff_coef）计算对应的 y 并返回，时间复杂度 $O(n)$。

plot 方法画牛顿插值后的函数图像。代码中的 10 个已知点的例子的插值结果如下：

Scipy 组件

在 Scipy 中进行一维多项式插值需要用到 BarycentricInterpolator 类和 barycentric_interpolate 函数，格式如下：

1
2
BarycentricInterpolator(xi, yi=None, axis=0)
barcentric_interpolate(xi, yi, x, axis=0)

其中，xi 和 yi 是已知的数据，多项式会通过这些数据点，xi 的数据通常按升序方式排列；axis 指定 yi 中数据所在的轴；x 是要进行插值计算的数据点。BarycentricInterpolator 类有方法 add_xi(xi[，yi]) 和 set_yi(yi[，axis])，用于添加 xi 的数据和更新 yi 的数据。

以下代码中用与前面例子相同的离散值，结果与前面的图完全一样。

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BarycentricInterpolator, barycentric_interpolate

x_point = np.array([0, 1, 2, 4, 5, 6, 7, 8, 10, 11])
y_point = np.array([1, 10, 2, 7, 11, 3, 6, 9, 4, 0])
n = x_point.shape[0]

# 画出插值后的函数图像
text_list = []
for i in range(n):
    # 画 n 个已知点
    text = "({}, {})".format(x_point[i], y_point[i])
    text_list.append(text)
    # annotate 第一个参数是文本内容 第二个参数是要标记的位置 第三个参数是文本标记位置
    plt.annotate(text_list[i]
                ,xy=(x_point[i], y_point[i])
                ,xytext=(x_point[i]+0.5, y_point[i]+0.5)
                )
    plt.scatter(x_point, y_point)

# 画图的点的横坐标
x_list = np.arange(min(x_point) - 0.01, max(x_point) + 0.003, 0.003)

# 用 BarycentricInterpolator 类的方式
fun_inter = BarycentricInterpolator(x_point) # 创建插值函数
fun_inter.set_yi(y_point)
y_list1 = fun_inter(x_list)

# 用 barycentric_interpolate 函数的方式
y_list2 = barycentric_interpolate(x_point, y_point, x_list)

plt.plot(x_list, y_list1)
plt.plot(x_list, y_list2)
plt.legend()
plt.show()

总结

我们熟知的泰勒公式其实就是对牛顿插值算法进行了改进：

$$ f(x) = f(x_{0}) + f'(x_{0})(x - x_{0}) + \frac{f''(x_{0})}{2!}(x - x_{0})^{2} + \frac{f'''(x_{0})}{3!}(x - x_{0})^{3} + \cdots $$

英国人沃利斯在研究插值的过程中发现了沃利斯公式、之后牛顿在沃利斯的基础上，研究插值的过程中推导出了广义二项式定理、之后同是英国人的泰勒在牛顿的基础上，也是在研究插值的过程中提出了泰勒公式。

三个英国人在 17 ~ 18 世纪的重大成果都是在研究插值问题时发现的，都是证明不严格但是先取得了广泛应用，且之后都由后世数学家陆续完成了严格证明。可见插值问题除了实用性，还具有理论价值。

Share