【笔记】Lagrange 插值与 Newton 插值

· · 算法·理论

插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。

插值法常用于函数拟合中，也就是给定函数 f(x) 的图像上的一些点，来拟合 f(x) 的表达式。

通常我们尝试用多项式去拟合 f(x)，即多项式插值。

多项式插值的一般形式是对于已知的 n+1 个点 (x_0,y_0),\ldots,(x_n,y_n)，求形如 f(x)=\sum\limits_{i=0}^na_ix^i 且对于 \forall i =0,1,\ldots,n 满足 f(x_i)=y_i 的多项式 f(x)。

常见的多项式插值算法有 Lagrange 插值与 Newton 插值。

Lagrange 插值法

对 \forall i=0,1,\ldots,n，构造 L_i(x)=\prod\limits_{j=0,j \neq i}^n\dfrac{x-x_j}{x_i-x_j}。

容易看出，L_i(x_i)=1,L_i(x_j)=0。

于是所求的多项式 f(x)=\sum\limits_{i=0}^ny_iL_i(x)。

朴素实现的时间复杂度是 O(n^2)。可以利用一些技巧优化到 O(n\log^2n)。

特别地，如果已知点的横坐标是连续整数，可以做到 O(n)。

Newton 插值法

不同于 Lagrange 插值法一次性构造整个函数，Newton 插值法一项一项构造，支持 O(n) 插入新数据点。

具体地，Newton 插值多项式写成

每增加一个新点，只需要多增加一项，无需重算整个多项式，这也是它支持 $O(n)$ 插入的原因。 系数 $a_0,a_1,\ldots,a_n$ 用对应阶差商进行表示。 从 $x_0$ 开始的一阶差商 $f[x_0,x_1]=\dfrac{y_1-y_0}{x_1-x_0}$。 从 $x_0$ 开始的二阶差商 $f[x_0,x_1,x_2] = \dfrac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}$。 以此类推，从 $x_0$ 开始的 $n$ 阶差商 $f[x_0,\ldots,x_n]=\dfrac{f[x_1,\ldots,x_n]-f[x_0,\ldots,x_{n-1}]}{x_n-x_0}$。 对于 $\forall k=0,1,\ldots,n$，$a_k$ 即为从 $x_0$ 开始的 $k$ 阶差商。 那么所求的多项式 $f(x)=\sum\limits_{i=0}^na_i\prod\limits_{j=0}^{i-1}(x-x_j)$。 显然，为了计算从 $x_0$ 开始的 $k$ 阶差商，必须先计算出从 $x_0,\ldots,x_j$ 开始的 $k-j$ 阶的差商 $(j=0,\ldots,k-1)$，因此时间复杂度是 $O(n^2)$。 更好理解的方式是构建差商表，非常直观。 从 $x_j$ 开始的 $k$ 阶差商的计算方法为 $f[x_j,\ldots,x_{j+k}]=\dfrac{f[x_{j+1},\ldots,x_{j+k}]-f[x_j,\ldots,x_{j+k-1}]}{x_{j+k}-x_j}

误差分析

两种插值方式的误差是相同的。

考虑对于未知函数 f(x) 上的 n+1 个已知点进行多项式插值得到的拟合函数 P_n(x)，设截断误差为 R_n(x)，则有

R_n(x)=f(x)-P_n(x)=\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod\limits_{i=0}^n(x-x_i)

其中 \xi 为所有插值节点 x_i 和求值点 x 的某个区间内，但不能确定具体值。

通常采用如下估计： \left|R_n(x)\right|\leq \dfrac{\max_{\xi\in[a,b]}|f^{(n+1)}(\xi)|}{(n+1)!}\left|\prod_{i=0}^n(x-x_i)\right|，其中 [a,b] 为插值区间。

于是只需要求出 f^{(n+1)}(x) 在 [a,b] 上的最大值即可。

如何计算求出的近似值 P_n(x) 精确到小数点后几位：若 |R_n(x)|<0.5\times10^{-k}，则小数点后 k 位是精确的。

例题

已知函数 f(x)=\sin x,f(0.7)=0.6442,f(0.9)=0.7833,f(1.1)=0.8912，计算 f(1.0) 的近似值，并进行误差分析。

Lagrange 插值法

Lagrange 插值多项式为：

L(x) = f(x_0)L_0(x) + f(x_1)L_1(x) + f(x_2) L_2(x)

其中基函数为：

L_0(x) = \dfrac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)}=\dfrac{(1.0 - 0.9)(1.0 - 1.1)}{(0.7 - 0.9)(0.7 - 1.1)}=-0.125 L_1(x) = \dfrac{(x - x_0)(x - x_2)}{(x_1 - x_0)(x_1 - x_2)}=\dfrac{(1.0 - 0.7)(1.0 - 1.1)}{(0.9 - 0.7)(0.9 - 1.1)}=0.75 L_2(x) = \dfrac{(x - x_0)(x - x_1)}{(x_2 - x_0)(x_2 - x_1)}=\dfrac{(1.0 - 0.7)(1.0 - 0.9)}{(1.1 - 0.7)(1.1 - 0.9)}=0.375

综上有 L(1.0) = 0.6442 \times (-0.125) + 0.7833 \times 0.75 + 0.8912 \times 0.375 = 0.84115。

Newton 插值法

构造差商表： x f(x) 一阶差商 二阶差商
0.7 0.6442 0.6955 -0.39
0.9 0.7833 0.5395 -
1.1 0.8912 - -

从 x_0 开始的一阶差商 f[x_0,x_1]=\dfrac{f(x_1)-f(x_0)}{x_1-x_0}=0.6955

从 x_1 开始的一阶差商 f[x_1,x_2]=\dfrac{f(x_2)-f(x_1)}{x_2-x_1}=0.5395

从 x_0 开始的二阶差商 f[x_0,x_1,x_2]=\dfrac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}=-0.39

近似值 P_2(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)=0.6442+0.6955 \times (x-0.7)-0.39 \times (x-0.7)(x-0.9)。

代入 x=1.0 可知 P_2(1.0)=0.84115，故计算出 \sin(1.0) 的近似值为 0.84115。

误差分析

计算截断误差 R(x)=f(x)-P_2(x)=\dfrac{f^{(3)}(\xi)}{3!}(x-x_0)(x-x_1)(x-x_2)。

又 f^{(3)}(x)=-\cos x,|f^{(3)}(x)| \le 1，故 |R_2(x)| \le \dfrac{1}{3!}|(x-x_0)(x-x_1)(x-x_2)|，代入 x=1.0,x_0=0.7,x_1=0.9,x_2=1.1 可知 |R_2(1.0)| \le 0.5 \times 10^{-3}。

故计算结果精确到小数点后三位。