【数值分析】插值

发表于 | 分类于 |

本篇作为数值分析的入门教程，介绍一下插值相关的常用算法。

概述

插值，是通过已知的数据点来估算未知数据点的方法。

它基于这样一个假设：数据点之间存在某种内在关系或规律，通过这些已知点，我们可以推测出未知点的值。

插值在数学、统计学、计算机图形学等领域都有广泛的应用:

  • 科学研究，插值可用于处理实验数据，提高数据质量
  • 图像处理，插值可用于图像缩放、旋转等操作中，实现图像的平滑变换
  • 机器学习，插值可用于数据预处理，提高模型的预测性能
  • …

本文将对常用的插值算法进行介绍。

插值定义

设 $y=f(x)$ 定义于 $[a,b]$, 并测得点 $a \leq x_0 \leq x_1… \leq x_n \leq b$ 上的值 $y_0,y_1…y_n$,

若存在 $p(x)$ 使得对任意 $x_i$, 有 $p(x_i)=y_i$, 则 $p(x)$ 为插值函数。

插值分类

根据插值函数的形式，插值函数大致可以分为以下几类：

  • 多项式插值：$p(x)$ 是多项式
  • 分段多项式插值：$p(x)$ 是分段多项式
  • 三角插值：$p(x)$ 是三角函数
  • 有理插值：$p(x)$ 是有理函数

多项式插值

设 $y=f(x)$ 定义于 $[a,b]$, 有 $n+1$ 个点 $a < x_0 < x_1… < x_n < b$ 上的值 $y_0,y_1…y_n$,

若存在次数不超过n的多项式 $p(x) = a_0 + a_1 x + … +a_n x^n$ 使得对任意 $x_i$, 有 $p(x_i)=y_i$, 则 $p(x)$ 为插值多项式。

插值多项式唯一性

定理：插值多项式 存在 且 唯一。

对于多项式基,$1,x,x^2,x^3…x^n$

假设插值多项式存在，那么应满足

$ a_0 + a_1 x_0 + … a_n x_0^n = y_0$

$ a_0 + a_1 x_1 + … a_n x_1^n = y_1$

…

$ a_0 + a_1 x_n + … a_n x_n^n = y_n$

显然方程组系数是范德蒙矩阵。

由范德蒙行列式不等于0，我们立刻得到：

对于一组不重复的点，插值多项式存在且唯一。插值多项式唯一性得证。

线性插值

给定两个点$(x_0,y_0)(x_1,y_1)$, 两点间的插值函数显然是一条直线, 斜率为$\frac{y_1-y_0} {x_1-x_0}$

$ p(x) = y_0 + \frac{y_1-y_0} {x_1-x_0} (x-x_0)$

$ p(x) = y_1 + \frac{y_1-y_0} {x_1-x_0} (x-x_1)$

求平均得到两点式： $p(x)=y_0\frac{x-x_1} {x_0-x_1}+y_1\frac{x-x_0} {x_1-x_0}$

若记作$p(x)=y_0 l_0(x)+y_1 l_1(x)$

有$l_0(x) + l_1(x) = 1$

两个点各自对应

  • $l_0(x_0)=1, l_0(x_1)=0$

  • $l_1(x_0)=0, l_1(x_1)=1$

因此我们临时得到对线性插值的一种描述方式：

  • 对于一组点，每个点对应一个自己的权重函数$l(x)$
  • 整体插值函数就是这些点的函数值对$l(x)$加权平均

抛物线插值

给定两个点$(x_0,y_0)(x_1,y_1)(x_2,y_2)$, 三点间的插值函数是二次曲线。

按照我们刚才的理解，我们对三个点有三个权重函数：

记作$p(x)=y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x)$

假设三个点各自对应如下方程，这样的函数组是否存在呢？

  • $l_0(x_0)=1, l_0(x_1)=0, l_0(x_2)=0$
  • $l_1(x_0)=0, l_1(x_1)=1, l_0(x_2)=0$
  • $l_2(x_0)=0, l_2(x_1)=0, l_2(x_2)=1$

显然在这三个点中选择任意一个点做分析是等价的。这种对称性说明只要$l_0$存在，$l_1,l_2$都存在。

对于多项式函数$l_0$, 我们观察到

$l_0(x)$ 在$x_1,x_2$都为0, 因此我们假设$l_0 = k (x-x_1)(x-x_2)$

$l_0(x)$ 在$x_0$ 为1, 有$l_0 = k (x_0-x_1)(x_0-x_2) = 1, k = \frac{1}{(x_0-x_1)(x_0-x_2)}$

因此$l_0 = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}$,存在且唯一(插值多项式唯一性)

拉格朗日插值

根据线性插值和抛物线插值的算法，我们对n次做推广，即为拉格朗日插值。

拉格朗日插值多项式

拉格朗日插值多项式 $ L_n = \sum _{k=0} ^n y_k l_k(x) $

其中基函数 $ l_k(x) = \frac{(x-x_0)(x-x_1)…(x-x_{k-1})(x-x_{k+1})…(x-x_n)} {(x_k-x_0)(x_k-x_1)…(x_k-x_{k-1})(x_k-x_{k+1})…(x_k-x_n)} = \prod _{i=0,i \neq k} ^n \frac{x-x_i} {x_k-x_i} $

引入记号 $ \omega _{n+1} (x) = \prod _{i=0} ^n {(x-x_i)} $

显然 $x = x_k$ 时 $ \omega _{n+1}’ (x_k) = ((x-x_k) \prod _{i=0,i \neq k} ^n {(x-x_i)})’ $

$= \prod _{i=0,i \neq k} ^n {(x-x_i)} + (x-x_k)(…) = \prod _{i=0,i \neq k} ^n {(x_k-x_i)} $

因此 $ \omega _{n+1}’ (x_k) = \prod _{i=0,i \neq k} ^n {(x_k-x_i)} $

拉格朗日插值余项

若在$[a,b]$上用$L(x)$近似$f(x)$，则其截断误差为 $R(x)=f(x)-L(x)$ ,也称为插值多项式的余项。

关于插值余项估计有以下定理：

定理：设$f^{(n)}(x)$在$[a，b]$上连续，$f^{(n+1)}(x)$在$(a，b)$内存在，$L(x)$是的插值多项式,

则对任何$x \in [a，b]$，插值余项 $R(x)=f(x)-L(x) = \frac{f^{(n+1)}(\xi)} {(n+1)!}\omega _{n+1} (x) , \xi \in (a,b) $

在所有节点上，$R_n(x_k)=0$,我们设 $R_n(x)=k(x) \omega _{n+1}(x)$

$f(x)-L(x)-R(x)=0$ 显然在$[a,b]$连续，在$[a,b]$上有n+1个零点， 假设$f^{(n+1)}(x)$存在

由罗尔定理，其一阶导至少有n个零点; 二阶导至少有n-1个零点…

n+1阶导至少有1个零点, 存在$\xi \in [a,b]$, 使得 $f^{(n+1)}(\xi) - (n+1)!k(x) = 0$

$k(x)$带入 $R_n(x)=k(x) \omega _{n+1}(x)$即可

如果 $f^{(n+1)}(x)$不存在呢？我们可以想办法求一个最大值 $M_{N+1} =max \lvert f^{(n+1)}(x) \rvert $ 来代替:

有 $\lvert R(x) \rvert \leq \frac{M_{N+1}} {(n+1)!} \lvert \omega _{n+1} (x)\rvert $

牛顿插值

我们发现对于拉格朗日插值，如果我们已经计算了n个点，我们再添加一个点进去，所有的系数都需要重新计算。

牛顿插值就解决了这个问题。

还记得我们开始做线性插值，使用了两点式来构造插值函数么？我们现在改用点斜式, 并记$ p_0(x) = y_0$

$ p_1(x) = p_0 + \frac{y_1-y_0} {x_1-x_0} (x-x_0)$

那么抛物线插值可以表示为

$ p_2(x) = p_1(x) +\frac{\frac{y_2-y_0}{x_2-x_0} - \frac{y_1-y_0}{x_1-x_0}} {x_2-x_1} (x-x_0)(x-x_1) $

我们定义差商为$\frac{f(x_1)-f(x_0)}{x_1-x_0}$, 我们发现抛物线的系数是差商的差商。

均差：

  • $f[x_0,x_k] = \frac{f(x_k)-f(x_0)}{x_k-x_0}$为函数$f(x)$关于$x_0,x_k$的一阶均差
  • $f[x_0,x_1,x_k] = \frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}$为函数$f(x)$关于$x_0,x_k$的二阶均差
  • $f[x_0,x_1,…,x_k] = \frac{f[x_0,…x_{k-2},x_k]-f[x_0,x_1,…,x_{k-1}]}{x_k-x_{k-1}}$为函数$f(x)$关于$x_0,x_k$的k阶均差

借助均差，我们有

  • $ p_1(x) = p_0 + f[x_0,x_1] (x-x_0)$
  • $ p_2(x) = p_1 + f[x_0,x_1,x_2] (x-x_0)(x-x_1)$

我们发现后一项可由前一项迭代出来，因此添加新点可以免去一部分重复的计算。

牛顿均差多项式

$ p_n(x) = a_0 + a_1 (x-x_0) + … + a_n (x-x_0)(x-x_1)…(x-x_n)$

其中系数为 $a_k = f[x_0,x_1,…x_k]$

余项满足 $f(x) = p(x)+R(x), R(x) = f[x,x_0,x_1,…x_n]\omega _{n+1}(x)$

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
import numpy as np
import matplotlib.pyplot as plt

def divided_diff(x, y):
    """计算差商"""
    n = len(y)
    coef = np.zeros([n, n])
    coef[:,0] = y
    
    for j in range(1, n):
        for i in range(n-j):
            coef[i][j] = (coef[i+1][j-1] - coef[i][j-1]) / (x[i+j] - x[i])
    
    return coef

def newton_poly(coef, x_data, x):
    """计算牛顿插值多项式的值"""
    result = coef[0][0]
    for k in range(1, len(x_data)):
        prod = coef[0][k]
        for j in range(k):
            prod *= (x - x_data[j])
        result += prod
    return result

def newton_interpolation(x_data, y_data, x):
    """牛顿插值法"""
    coef = divided_diff(x_data, y_data)
    return newton_poly(coef, x_data, x)

# 示例数据
x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([1, 3, 2, 5, 7])

# 执行牛顿插值
x = np.linspace(1,5,100)
y = [newton_interpolation(x_data, y_data, xi) for xi in x]
plt.plot(x,y)
牛顿前插公式

上面给出节点任意分布的情况，实际中等距节点十分普遍。

对于等距节点$x_k=x_0+kh$，h称为步长，公式可以简化。

我们对步长h的情况做如下定义：

  • 一阶差分算子 $\Delta f_k = f_{k+1} - f_k$
  • 二阶差分算子 $\Delta^2 f_k = \Delta f_{k+1} -\Delta f_k$
  • n阶差分算子 $\Delta^n f_k = \Delta^{n-1} f_{k+1} -\Delta^{n-1} f_k$

为了方便，我们定义两个常用算子：

  • 不变算子： $If_k = f_k$
  • 位移算子： $Ef_k = f_{k+1}$

于是差分算子可以表示为：

  • $\Delta f_k = E f_k - I f_k = (E-I)f_k$
  • $\Delta ^n f_k = (E-I)^n f_k = \sum _{j=0} ^n (-1)^j C_n ^j E ^{n-j} f_k$

$f[x_k,x_{k+1}] = \frac {f_{k+1}-f_k}{x_{k+1}-x_k} = \frac {\Delta f_{k}}{h}$

$f[x_k,…,x_{k+m}] = \frac {f_{k+1}-f_k}{x_{k+1}-x_k} = \frac {\Delta ^m f_{k}}{m! h^m}$

用上式代替均差，令$x=x_0+th$,得到牛顿前插公式：

$p_n(x_0+th) = f_0 + t\Delta h f_0 + \frac{t(t-1)}{2!}\Delta ^2 f_0 + … + \frac{t(t-1)…(t-n+1)}{n!}\Delta ^n f_0$

余项为 $R_n(x_0+th) = \frac{t(t-1)…(t-n)}{(N+1)!} h ^{n+1} f^{n+1}(\xi), \xi \in (x_0,x_n) $

埃尔米特插值

我们之前讨论的都是满足再节点上函数中相等，实际中我们有时需要一些节点导数/二阶导数/…等于特定数值。

为此我们引入埃尔米特插值。

重节点均差

假设$f\in C^1 [a,b]$, $ \lim_{x \to x_0} \frac{f(x) - f(x_0)}{x-x_0} = f^{‘}(x_0) $

因此定义重节点:

  • 一阶均差 $ f[x_0,x_0] = \lim_{x \to x_0} f[x_0,x] = f^{‘}(x_0) $
  • 二阶均差 $ f[x_0,x_0,x_0] = \lim_{x_1 \to x_0, x_2 \to x_0} f[x_0,x_1,x_2] = \frac{1}{2} f^{‘’}(x_0) $
  • n阶均差 $ f[x_0,x_0,…,x_0] = \lim_{x_i \to x_0} f[x_0,x_1,…,x_n] = \frac{1}{n!} f^{(n)}(x_0) $

由重节点均差带入牛顿均差插值多项式，得到泰勒多项式：

$ p_n(x) = f(x_0) + f^{‘}(x_0)(x-x_0) +\frac{f^{(n)}(x_0)}{n!} (x-x_0)^n $

它实际是在点$x_0$附近逼近函数的带导数的插值多项式，满足$p^{(k)}(x_0) = f^{(k)}(x_0)$

余项为$R_n(x) =\frac{f^{(n+1)}(\xi)}{(n+1)!} (x-x_0)^{(n+1)}, \xi \in (a,b)$

给出m+1个插值条件，可以得到不超过m次的埃尔米特插值多项式。

由于导数条件不同，这里不给出通用公式。可灵活使用基函数法或列方程组用矩阵求插值多项式。

分段低次插值

龙格现象

有一些人猜测n的次数越高，逼近函数效果就越好，是这样么？

一个典型的反例由龙格提出：在$[-5,5]$上$f(x) = \frac {1}{1+x^2}$

插值点越多，在两边震荡越明显。

为了解决高次插值的龙格现象，我们提出分段低次插值。

分段线性插值

分段线性插值其实就是每段两点之间做线性插值，(折线图)。 误差估计可直接由每段线性插值的余项的最大值得出。 这种插值比较少用，在此不详细叙述了。

分段三次埃尔米特插值

分段三次埃尔米特插值已知各点的函数值和导数值，其实就是每段两点之间做三次埃尔米特插值。 误差估计可直接由每段三次埃尔米特插值的余项的最大值得出。 这种插值需求的信息很多，但只有一阶光滑，因此也比较少用，在此不详细叙述了。 为改进其缺点，提出三次样条插值。

三次样条插值

边界条件

三次样条函数对于不同边界条件的算法是不同的。

  • 第一类边界条件：给定边界点一阶导数
  • 第二类边界条件：给定边界点二阶导数
    • 自然边界条件：两端点二阶导数为0，是第二类边界条件特例
  • 周期边界条件：在一个周期两端点函数值，一阶导数，二阶导数相等。
三次样条求解

设子区间 $[x_i,x_{i+1}]$ 上插值函数$S_i$为$ S_i $, 令 $h_i(x) = x-x_{i}$

$ S_i = a_i(x-x_i)^3 + b_i(x-x_i)^2 + c_i(x-x_i) + d_i $

$ S_i = a_i(h_i)^3 + b_i(h_i)^2 + c_i(h_i) + d_i $

$ S_i’ = 3a_i(h_i)^2 + 2b_i(h_i) + c_i $

$ S_i’’ = 6a_i(h_i) + 2b_i $

在起点处，有

$ S_{0}(x_{0}) = y_0$

$ d_0 = y_0 $

在子区间右端点处，有

$ S_{i+1}(x_{i+1}) = S_i(x_{i+1})$

$ d_{i+1} = a_i(x_{i+1}-x_{i})^3 + b_i(x_{i+1}-x_{i})^2 + c_i(x_{i+1}-x_{i}) + d_{i} $

在当前子区间终点和下一个子区间的起点上，一阶导数和二阶导数连续，因此有

$ S_i’(x_{i+1}) = S_{i+1}’(x_{i+1})$

$ c_{i+1} = 3a_i(x_{i+1}-x_{i})^2 + 2b_i(x_{i+1}-x_{i}) + c_i $

$ S_i’‘(x_{i+1}) = S_{i+1}’‘(x_{i+1})$

$ 2b_{i+1} = 6a_i(x_{i+1}-x_{i}) + 2b_i $

假设已经知道 $b_{i}, c_{i}, d_{i}$ 则可以根据 $S_{i}(x_{i+1}) = y_{i}$ 求出 $a_{i}$

对于 n+1 个点，存在 n 个多项式，需要 4n 个方程； 由上文，递推方程有 3n-3 个，过 n+1 个点对应 n+1 个方程， 共 4n-2 个方程。 因此只需两个方程，即上面三种边界条件。

若取 $ M_i = S’‘_{i}(x_i) = 2b_i $ 对于第一类边界条件，有三对角方程 $ AM = D $ ，M 为 二阶导数的列向量，

列向量 D 中元素

\[d_0 = \frac{6(-f'_0+f[x_{0}, x_1])}{h_{0}} , d_n = \frac{6(f'_n-f[x_{n-1}, x_n])}{h_{n-1}}\]

系数矩阵 A 如下，

\[\begin{pmatrix} 2 & \lambda _0 & 0 & 0 & \cdots & 0 \\ \mu _1 & 2 & \lambda _1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \mu _{n-1} & 2 & \lambda _{n-1} \\ 0 & 0 & 0 & \cdots & \mu _n & 2 \end{pmatrix}\]

其中参数

\[\mu _n = \frac{h_{j-1}}{h_{j-1}+h_j} , \mu _j = \frac{h_{j-1}}{h_{j-1}+h_j}\] \[\lambda _0 = 1 , \lambda _j = \frac{h_j}{h_{j-1}+h_j}\]

第二类第三类同理，这里就不展开了。

需要注意的是，三对角矩阵可以用追赶法来快速求解，即先消元成上/下三角矩阵，再从最后开始依次求解其他(矩阵中相邻的非0元素)的解。