数值计算方法 Chapter1. 插值插值问题的本质其实就是： - 给定一堆采样点，然后构造一个函数来对这堆采样点背后的

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

数值计算方法 Chapter1. 插值

1. 定义

插值问题的本质其实就是：

给定一堆采样点，然后构造一个函数来对这堆采样点背后的真实函数表达进行拟合。

即是说，找一条经过这一堆采样点的曲线来对这些采样点背后的真实函数曲线进行描述。

我们给出书中的定义如下：

$f(x)$ 为定义在区间 $[a,b]$ 上的函数， $x_0, ..., x_n$ 为 $[a,b]$ 上的 $n+1$ 个互不相同的点， $\Phi$ 为给定的某一函数类，若 $\Phi$ 上有函数 $\phi{x}$ ，满足 $\phi(x_i) = f(x_i), \ i=0,1,...,n$ 则称 $\phi(x)$ 为 $f(x)$ 关于节点 $x_0, ..., x_n$ 在 $\Phi$ 上的插值函数。称点 $x_0, ..., x_n$ 为插值节点；称 $(x_i, f(x_i)), \ i=0,1,...,n$ 为插值型值点，简称型值点或者插值点； $f(x)$ 称为被插函数。

2. Lagrange插值

1. 定义 & 实现

Lagrange插值公式本质上就是用一个 $n$ 阶函数来拟合这些采样点，因此，我们事实上就是要解如下方程组：

\left\{ \begin{aligned} y_0 &= a_0 + a_1 x_0 + a_2 x_0^2 + ... + a_n x_0^n \\ y_1 &= a_0 + a_1 x_1 + a_2 x_1^2 + ... + a_n x_1^n \\ ... \\ y_n &= a_0 + a_1 x_n + a_2 x_n^2 + ... + a_n x_n^n \end{aligned} \right.

当满足 $x_0, ..., x_n$ 互异时，上述方程的解是唯一的。

而Lagrange直接给出了一个插值函数，令其满足上述 $n+1$ 个方程，亦即Lagrange插值函数：

L_n(x) = \sum_{i=0}^{n} l_i(x) y_i

其中， $l_i(x)$ 满足条件：

\left\{ \begin{aligned} l_i &= 1 \text{ if } x = x_i \\ l_i &= 0 \text{ if } x = x_j, \ j \neq i \end{aligned} \right.

由是，可以直观地给出 $l_i(x)$ 的一种实现方式如下：

l_i(x) = \Pi_{j \neq i} \frac{x - x_j}{x_i - x_j}

2. 伪代码实现

我们在python当中给出Lagrange插值函数的具体代码实现如下：

def lagrange_fn(xlist, ylist):
    assert(len(xlist) == len(ylist))
    assert(len(set(xlist)) == len(xlist))
    n = len(xlist)
    llist = []
    for i in range(n):
        li = ylist[i]
        for j in range(n):
            if j == i:
                continue
            li /= (xlist[i] - xlist[j])
        llist.append(li)
    
    def fn(x):
        s = 1
        for xi, yi in zip(xlist, ylist):
            if x == xi:
                return yi
            s *= (x - xi)
        res = 0
        for xi, li in zip(xlist, llist):
            res += li * s / (x - xi)
        return res
    
    return fn

3. 误差分析

对于n次插值多项式的误差，有如下定理：

设 $L_n(x)$ 是 $[a,b]$ 上过 $(x_i, f(x_i)), i=0,1,...,n$ 的 $n$ 次插值多项式， $x_i \in [a, b]$ ， $x_i$ 互不相同，当 $f\in C^{n+1}[a,b]$ 时，则存在 $\xi \in [a,b]$ 使得插值多项式的误差为： $R_n(x) = \frac{f^{n+1}(\xi)}{(n+1)!}\Pi_{i=0}^{n}(x-x_i)$

3. Newton插值

1. 定义 & 实现

Newton插值公式和Lagrange插值公式本质上来说是一样的，都是用一个 $n$ 阶的多项式来对采样点进行拟合。

由前所述，由于n阶函数的解是唯一的，所以Newton插值公式本质上来说和Lagrange插值公式是完全等价的。

他们的区别在于具体的实现思路，Lagrange插值是平权地对每一个点进行描述，而Newton插值的思路则是通过残差的方式进行实现。

即是说，要对第 $k$ 个点进行拟合，只需要先拟合前 $k-1$ 个点进行拟合，得到一个 $k-1$ 维的函数，然后对于第k个点，就会有残差 $y_k - f_{k+1}(x_k)$ ，然后我们用一个k维的多项式来填补这个差距即可。

综上，我们可以写出Newton插值公式的表达式如下：

N(x) = f(x_0) + \sum_{i=1}^{n} f[x_0, x_1, ..., x_{i-1}]\Pi_{j=0}^{i-1}(x-x_{j})

其中， $f[x_0, x_1, ..., x_k]$ 就是 $k$ 阶差商，其定义可以通过迭代关系得到：

\left\{ \begin{aligned} &f[x_k] = f(x_k) \\ &... \\ &f[x_0, ..., x_k] = \frac{f[x_1, ..., x_k] - f[x_0, ..., x_{k-1}]}{x_k - x_0} \end{aligned} \right.

更进一步的，我们可以给出书中列举出的关于差商的几个性质如下：

性质1.1： $k$ 阶差商 $f[x_0, ..., x_k]$ 是由函数值 $f(x_0), f(x_1), ..., f(x_k)$ 的线性组合而成，即： $f[x_0, x_1, ... ,x_k] = \sum_{i=0}^{k} \frac{1}{\Pi_{j \neq i} (x_i - x_j)} f(x_i)$

性质1.2：若 $i_0, i_1, ..., i_k$ 为 $0, 1, ..., k$ 的任一排列，则有： $f[x_0, x_1, ..., x_k] = f[x_{i_0}, x_{i_1}, ..., x_{i_k}]$

性质1.3：若 $f(x)$ 为 $m$ 次多项式，则 $f[x_0, x_1, ..., x_{k-1}, x]$ 为 $m-k$ 次多项式。

2. 伪代码实现

同样的，我们给出python代码演示如下：

def newton_fn(xlist, ylist):
    assert(len(xlist) == len(ylist))
    assert(len(set(xlist)) == len(xlist))
    n = len(xlist)
    nmatrix = [[0 for _ in range(n)] for _ in range(n)]
    for i in range(n):
        nmatrix[i][0] = ylist[i]
        for j in range(1, i+1):
            nmatrix[i][j] = (nmatrix[i][j-1] - nmatrix[i-1][j-1]) / (xlist[i] - xlist[i-j])
    nlist = [nmatrix[i][i] for i in range(n)]
            
    def fn(x):
        s = 1
        res = 0
        for i in range(n):
            res += nlist[i] * s
            s *= (x-xlist[i])
        return res
    
    return fn

3. 误差分析

如前所述，由于Newton插值公式和Lagrange插值公式本质上是同一个多项式的两种构造方法，所以他们的误差分析是完全相同的，这里也就不再多做赘述了。

4. 分段插值

1. 定义 & 实现

上述Lagrange插值和Newton插值本质上都是用一个n阶方程来对 $n+1$ 个点进行描述，而高阶方程存在一个极大的问题就是过拟合问题，也就是说，可能为了对某些点进行拟合而导致函数去现在插值区间内剧烈震荡。

上述现象又称为Runge现象，一个典型的例子就是 $f(x) = \frac{1}{1+25x^2}, x \in [-1, 1]$ 。给出的拟合点越多，函数在边界附近的震荡越厉害。

因此，这里给出另外一种插值方法，即分段插值方法，其思路极其暴力，即首先对点进行排序处理，然后对每两个邻接点之间线性都采用线性连接。

此时，对于任何一个待求的点，我们只要找到其邻接的两个点 $x_i, x_{i+1}$ ，就可以快速地给出其函数估计值：

f(x) = \frac{x - x_{i+1}}{x_i - x_{i+1}}f(x_i) + \frac{x - x_{i}}{x_{i+1} - x_{i}}f(x_{i+1})

2. 伪代码实现

同样的，我们给出分段插值的python代码示例如下：

import bisect

def segment_fn(xlist, ylist):
    assert(len(xlist) == len(ylist))
    assert(len(set(xlist)) == len(xlist))
    xy = sorted([(x, y) for x, y in zip(xlist, ylist)])
    xlist = [x[0] for x in xy]
    ylist = [x[1] for x in xy]
            
    def fn(x):
        assert(xlist[0] <= x <= xlist[-1])
        i = bisect.bisect_left(xlist, x)
        y = (x - xlist[i])/(xlist[i+1]-xlist[i]) * ylist[i+1] + (x - xlist[i+1])/(xlist[i]-xlist[i+1]) * ylist[i]
        return y
    
    return fn

3. 误差分析

分段插值本质上来说就是在每一个小的区间上都采用线性拟合，因此，根据上述Lagrange插值的误差分析，其误差总可以表示为：

\begin{aligned} \delta(x) &= \frac{f^{(2)}(\xi)}{2!}(x-x_i)(x-x_{i+1}) \\ &\leq \frac{M_2}{2}|(x-x_i)(x-x_{i+1})| \\ &\leq \frac{M_2}{2} \cdot \frac{1}{4}(x_{i+1}-x_{i})^2 \\ &= \frac{M_2}{8}(x_{i+1}-x_{i})^2 \end{aligned}

其中， $M_2 = max|f^{(2)}(x)|, x \in [a, b]$ 。

5. 三次样条插值

1. 定义

如前所述，Lagrange插值和Newton插值平滑，但是容易过拟合，反之分段插值可以有效的防止过拟合，但是在连接处不够平滑，如果采样点不够充分，则拟合效果可能不太好。

而三次样条函数则是结合了上述几种方式的优点，它依然采用的是分段插值的方式，从而避免过拟合，但是，为了增加平滑性，他在两点之间不再使用线性连接，而是采用一个三次函数，然后限制连接处位置的一阶导数和二阶导数连续，从而保证了拟合曲线整体的平滑性。

给出书中的定义如下：

给定区间 $[a,b]$ 上的 $n+1$ 个节点 $a=x_0 < x_1 < ... < x_n = b$ 和这些点上的函数值 $f(x_i) = y_i, i = 0, 1, ..., n$ 。若 $S(x)$ 满足 $S(x_i) = y_i, i = 0,1,...,n$ $S(x)$ 在每个小区间 $[x_i, x_{i+1}]$ 上至多是一个三次多项式， $S(x)$ 在 $[a,b]$ 上有连续的二阶导数，则称 $S(x)$ 为 $f(x)$ 关于剖分 $a=x_0 < x_1 < ... < x_n = b$ 的三次样条插值函数，称 $x_0, x_1, ..., x_n$ 为样条节点。

但是，这里比较特殊的是，这里一共会有 $n$ 个区间， $4n$ 个待定参数，而通过 $n+1$ 个节点的值以及邻接点上的一阶导和二阶导连续条件，我们一共只有 $2 + 2(n-1) + 2(n-1) = 4n-2$ 个限制条件，不足以对全部的 $4n$ 个参数进行完全求解，因此，我们还需要额外提供两个边界条件，才能够完成所有的 $4n$ 个参数的求解。

这里的参数求解方式需要联合矩阵的求解，因此这里就不给出代码示例了。