Night关于数学的杂谈-插值法

插值法是什么

插值，就是给定一定的离散数据点，范围内估计新数据点的过程或方法。在这个过程中，我们当然希望得到一个连续的光滑曲线同时经过所有的 $(x_i,y_i)$ ，并求得该曲线在需要求值的点上的值。
具体定义如下：
给定 $n$ 个离散数据点 $(x_i,y_i) \ (i \in [1,n])$，对于 $x \ ( x \ne x_i , i \in [1,n] )$，求 $x$ 所对应的 $y$ 值称之为 内插。
$f(x)$ 为定义在区间 $[a,b]$ 上的函数。$x_1,x_2,\cdots ,x_n$ 为 $[a,b]$ 上 $n$ 个互不相同的点，$G$ 为给定的某一函数类。若 $G$ 上有函数 $g(x)$ 满足：

$$g(x_i)=f(x_i)\ , \ i \in [1,n]$$ 则称 $g(x)$ 为 $f(x)$ 关于节点 $x_1,x_2,\cdots ,x_n$ 在 $G$ 上的插值函数。
（摘自维基百科）

---

牛顿插值法

差商

给定函数 $f(x)$ 和插值节点 $x_0,x_1,\cdots ,x_n$，用 $f[x_0,x_1,\cdots ,x_k]$ 表示 $f(x)$ 关于节点 $x_0,x_1,\cdots ,x_k$ 的 $k$ 阶差商(k-th Difference Quotient) $(k=1,2,\cdots,n)$ ，它们可递归定义为

$$f[x_0,x_1,\cdots ,x_k]=\frac{f[x_1,x_2,\cdots ,x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_0}$$
其中 $f(x)$ 关于 $x_i$ 的 $0$ 阶差商为 $f(x_i)$。
根据差商的定义，差商具有如下性质：
1、 $$f[x_0,x_1,\cdots,x_k]=\sum_{i=0}^{k}{\frac{f(x_i)}{\prod_{j=0,j \ne i}^{k}{(x_i-x_j)}}}$$
证明的话可以用数学归纳法直接爆肝就行了，此处省略（其实这个性质下面并不会用到，主要是引入对称这个性质）。

2、由上面那个公式可以显然看出，$x_i$ 顺序的改变并不会影响差商的值，因此我们称差商具有对称性。

牛顿插值公式

于是我们利用差商的性质，可以推出牛顿插值公式。方法如下：
首先由差商的定义，以及其对称的性质，可以得到：

$$f[x,x_0]=\frac{f(x)-f(x_0)}{x-x_0} \ \ \rightarrow \ \ f(x)=f(x_0)+f[x,x_0](x-x_0) \ \ \ \ (a)$$ $$f[x,x_0,x_1]=\frac{f[x,x_0]-f[x_0,x_1]}{x-x_1}$$ $$\rightarrow \ \ f[x,x_0]=f[x_0,x_1]+f[x,x_0,x_1](x-x_1) \ \ \ \ (b)$$ $$\cdots$$ $$f[x,x_0,\cdots,x_n]=\frac{f[x,x_0,\cdots ,x_{n-1}]-f[x_0,x_1,\cdots ,x_n]}{x-x_n}\ \ \rightarrow \ \ f[x,x_0,\cdots ,x_{n-1}]=f[x_0,x_1,\cdots ,x_n]+f[x,x_0,\cdots,x_n](x-x_n) \ \ \ \ (c)$$
我们发现，对于式子 $(b)$，我们将其乘上 $(x-x_0)$，那么就会有
$$(a)+(b)\times(x-x_0)\ ：\ \ \ f(x)+f[x,x_0](x-x_0)=f(x_0)+f[x,x_0](x-x_0)+f[x_0,x_1](x-x_0)+f[x,x_0,x_1](x-x_0)(x-x_1)$$
递推地看，我们发现一个规律，每个式子 $(c)$ 如果乘上一个 $(x-x_0)(x-x_1)\cdots(x-x_{n-1})$ 以后可以上下相加来抵消一些项。
因此对于每一个 $(c)$ 式，我们将其乘上 $(x-x_0)(x-x_1)\cdots(x-x_{n-1})$ 后求和，可以得到这样的式子 $$f(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) + \cdots +f[x_0,x_1,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})+f[x,x_0,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})(x-x_n)$$
这就是 牛顿插值公式。
形式化地，有 $f(x)=P_n(x) + R_n(x)$ ，并记 $$P_n(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) + \cdots +f[x_0,x_1,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})$$ $$R_n(x)=f[x,x_0,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})(x-x_n)$$
并称 $P_n(x)$ 为 牛顿插值多项式， $R_n(x)$ 为 牛顿插值余项。

复杂度

显然我们预处理差商的时候需要 $O(n^2)$ 的复杂度，而单次计算的时候只需要用秦九韶算法 $O(n)$ 求出。

优缺点

牛顿插值法可以在 $O(n)$ 的复杂度计算单次的值， 并且当插值点增加的时候，我们也只需要 $O(n)$ 的复杂度刷新一遍差商即可，这就是牛顿插值法的优越性。缺点的话，就是计算较为复杂，式子也不好背。
（其实缺点是这个算法无法在特殊情况下优化到 $O(n)$，如何这样请看下文）

实现

在插值这个过程中，我们只需要求出 $P(x)$ 的值即可，余项只是用来判误差的。
观察这个式子，我们发现对于 $1$ 这个因式，所有的项都有（这不是废话）；对于 $(x-x_0)$ 这个因式，从第二项到最后一项都有；对于 $(x-x_0)(x-x_1)$ 这个因式，从第三项到最后一项都有 $\cdots$
发现什么没有？这玩意非常像那个叫秦九韶算法的玩意。因此我们考虑预处理差商，然后利用秦九韶算法从内往外求出 $P(x)$ 的值。代码如下：

//Night's template
#include <bits/stdc++.h>
#define R register
#define LL long long
template<class TT>inline TT Max(R TT a,R TT b){return a<b?b:a;}
template<class TT>inline TT Min(R TT a,R TT b){return a<b?a:b;}
using namespace std;
template<class TT>inline void read(R TT &x){
    x=0;R bool f=false;R char c=getchar();
    for(;c<48||c>57;c=getchar())f|=(c=='-');
    for(;c>47&&c<58;c=getchar())x=(x<<1)+(x<<3)+(c^48);
    (f)&&(x=-x);
}
//end template
const int maxn = 2020;
struct node{
    double x,y;
}P[maxn];//储存点的结构题体
int n;
double a[maxn];//差商，作为系数
double x;//求值点
int main(){
    read(n);
    for(R int i=0;i<n;++i){
        scanf("%lf%lf",&P[i].x,&P[i].y);
    }//输入数据点
    scanf("%lf",&x);//求值点
    for(R int i=0;i<n;++i)a[i]=P[i].y;
    for(R int i=0;i<n;++i){
        for(R int j=n-1;j>i;j--){
            a[j]=(a[j]-a[j-1])/(P[j].x-P[j-1-i].x);
        }
    }//处理差商

    R double tmp=1,ans=a[0];
    for(R int i=0;i<n;++i){
        tmp*=(x-P[i].x);
        ans+=tmp*a[i+1];
    }
    printf("%.5lf\n",ans);
    return 0;
}

---

拉格朗日插值法

推导分析

由于我们要构造的多项式 $L(x)$ 需要满足 $L(x_i)=y_i$，因此我们考虑构造这样一个式子，使得 $f_i(x_i)=y_i$ 且其余的 $f_i(x_j)\ (j \ne i)$ 都为 $0$。
显然如果取 $f_i(x)$ 有 $(x-x_j)\ (j \ne i)$ 这些因式，那么就可以保证这一点。
因此取

$$f_i(x) = \prod_{j=0,\ j \ne i}^{n}{\frac{x-x_j}{x_i-x_j}}$$
那么我们要构造的多项式即是
$$L(x)=\sum_{i=0}^{n}{f_i(x)\ y_i}$$
那么 $L(x)$ 称为 拉格朗日插值多项式， $f_i(x)$ 称为 拉格朗日基本多项式（或者插值基函数）。

复杂度

显然求每一个 $f_i(x)$ 是 $O(n)$ 的，于是求 $L(x)$ 是 $O(n\times n)=O(n^2)$ 的。因此对于每个插值点，拉格朗日插值法的复杂度是 $O(n^2)$ 的。

优缺点

拉格朗日插值多项式结构整齐，计算十分方便（比牛顿插值法好记多了）。但是在实际计算中，当插值点有所修改时，所对应的整个多项式就需要全部重新计算，复杂度是 $O(n^2)$ 的（相比之下牛顿插值法只需要 $O(n)$），因此我们就要用重心拉格朗日插值法或牛顿插值法来代替。此外，当插值点比较多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差。这类现象也被称为龙格现象，相同地，牛顿插值法由于也是一个 $n$ 次的多项式，也会出现同样的问题。解决的办法是分段用较低次数的插值多项式。

重心插值法

再观察一下这个式子

$$f_i(x) = \prod_{j=0,\ j \ne i}^{n}{\frac{x-x_j}{x_i-x_j}}$$

$$L(x)=\sum_{i=0}^{n}{(f_i(x)\times y_i)}$$
我们会发现我们每个 $f_i(x)$ 都重复计算了这样一个东西
$$\prod_{j=0,j\ne i}^{n}{(x-x_j)}$$
于是乎，我们记
$$g(x)=\prod_{i=0}^{n}{(x-x_i)}$$
来解决这个问题，这样一来 $L(x)$ 就会变成这样：
$$L(x)=g(x)\sum_{i=0}^{n}{\frac{y_i}{\left (\prod_{j\ne i}{(x_i-x_j)}\right)\times (x-x_i)}}$$
为了解决拉格朗日插值法每次加入插值点都需要重新计算 $f_i(x)$ 的问题，对于每一个 $f_i($ 我们记重心权
$$w_i=\frac{y_i}{\prod_{j\ne i}{(x_i-x_j)}}$$
则
$$L(x)=g(x)\sum_{i=0}^{n}{\frac{w_i}{x-x_i}}$$
这样，求预处理 $w_i$ 的复杂度为 $O(n^2)$，单次求值 $L(x)$ 的复杂度为 $O(n)$。
而对于动态加点的问题，我们只需要 $O(n)$ 地修改 $w_i$ 就可以了，避免了 $O(n^2)$ 的复杂运算，具体做法呢，就是把所有的 $w_i$ 除以 $x_i-x_{n+1}$ 得到新的 $w_{n+1}$，这样一来我们的添加插值点的复杂度就变成 $O(n)$ 的了。

---

线性插值法

所谓线性插值法，并不是指你把插值点两两直线连起来的那个插值法，而是能在保证插值点 $x_i=x_0+i$ 的情况下将复杂度优化到 $O(n)$ 的重心拉格朗日插值法。

推导如下

由上面讲的重心插值法，我们有

$$L(x)=g(x)\sum_{i=0}^{n}{\frac{w_i}{x-x_i}}$$
其中的 $g(x)$ 显然可以 $O(n)$ 预处理，因此我们要想办法 $O(n)$ 地计算 $$\sum_{i=0}^{n}{\frac{w_i}{x-x_i}}$$
显然前面那个 $\sum$ 是省不掉的啊，所以我们考虑如何 $O(1)$ 求
$$w_i=\frac{y_i}{\prod_{j\ne i}{(x_i-x_j)}}$$
我们考虑把下面的 $\prod$ 分成 $j<i$ 与 $j>i$ 两类，那么
$$w_i=\frac{y_i}{\prod_{j=0}^{i-1}{(x_i-x_j)}\times \prod_{j=i+1}^{n}{(x_i-x_j)}}$$
由于 $x_i=x_0+i$，那么
$$w_i=\frac{y_i}{\prod_{j=1}^{i}{j}\times \prod_{j=1}^{n-i}{-j}}$$
这样一来，我们可以预处理
$$v_i=\prod_{j=1}^{i}{j}$$
那么
$$w_i=\frac{y_i}{v_j \times -v_{n-i}}$$
于是我们就可以 $O(n)$ 地预处理 $v_i$，然后 $O(1)$ 地求出 $w_i$，然后 $O(n)$ 地计算
$$L(x)=g(x)\sum_{i=0}^{n}{\frac{w_i}{x-x_i}}$$
于是此时，我们的插值法就变成线性的啦（耶）
（代码与例题未完待续…
（可能会没空写…