插值法是什么
插值,就是给定一定的离散数据点,范围内估计新数据点的过程或方法。在这个过程中,我们当然希望得到一个连续的光滑曲线同时经过所有的
(xi,yi)
(
x
i
,
y
i
)
,并求得该曲线在需要求值的点上的值。
具体定义如下:
给定
n
n
个离散数据点 ,对于
x (x≠xi,i∈[1,n])
x
(
x
≠
x
i
,
i
∈
[
1
,
n
]
)
,求
x
x
所对应的 值称之为 内插。
f(x)
f
(
x
)
为定义在区间
[a,b]
[
a
,
b
]
上的函数。
x1,x2,⋯,xn
x
1
,
x
2
,
⋯
,
x
n
为
[a,b]
[
a
,
b
]
上
n
n
个互不相同的点, 为给定的某一函数类。若
G
G
上有函数 满足:
(摘自维基百科)
牛顿插值法
差商
给定函数 和插值节点
x0,x1,⋯,xn
x
0
,
x
1
,
⋯
,
x
n
,用
f[x0,x1,⋯,xk]
f
[
x
0
,
x
1
,
⋯
,
x
k
]
表示
f(x)
f
(
x
)
关于节点
x0,x1,⋯,xk
x
0
,
x
1
,
⋯
,
x
k
的
k
k
阶差商(k-th Difference Quotient) ,它们可递归定义为
其中 f(x) f ( x ) 关于 xi x i 的 0 0 阶差商为 。
根据差商的定义,差商具有如下性质:
1、
证明的话可以用数学归纳法直接爆肝就行了,此处省略(其实这个性质下面并不会用到,主要是引入对称这个性质)。
2、由上面那个公式可以显然看出, xi x i 顺序的改变并不会影响差商的值,因此我们称差商具有对称性。
牛顿插值公式
于是我们利用差商的性质,可以推出牛顿插值公式。方法如下:
首先由差商的定义,以及其对称的性质,可以得到:
我们发现,对于式子 (b) ( b ) ,我们将其乘上 (x−x0) ( x − x 0 ) ,那么就会有
递推地看,我们发现一个规律,每个式子 (c) ( c ) 如果乘上一个 (x−x0)(x−x1)⋯(x−xn−1) ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) 以后可以上下相加来抵消一些项。
因此对于每一个 (c) ( c ) 式,我们将其乘上 (x−x0)(x−x1)⋯(x−xn−1) ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) 后求和,可以得到这样的式子
这就是 牛顿插值公式。
形式化地,有 f(x)=Pn(x)+Rn(x) f ( x ) = P n ( x ) + R n ( x ) ,并记
并称 Pn(x) P n ( x ) 为 牛顿插值多项式, Rn(x) R n ( x ) 为 牛顿插值余项。
复杂度
显然我们预处理差商的时候需要 O(n2) O ( n 2 ) 的复杂度,而单次计算的时候只需要用秦九韶算法 O(n) O ( n ) 求出。
优缺点
牛顿插值法可以在
O(n)
O
(
n
)
的复杂度计算单次的值, 并且当插值点增加的时候,我们也只需要
O(n)
O
(
n
)
的复杂度刷新一遍差商即可,这就是牛顿插值法的优越性。缺点的话,就是计算较为复杂,式子也不好背。
(其实缺点是这个算法无法在特殊情况下优化到
O(n)
O
(
n
)
,如何这样请看下文)
实现
在插值这个过程中,我们只需要求出
P(x)
P
(
x
)
的值即可,余项只是用来判误差的。
观察这个式子,我们发现对于
1
1
这个因式,所有的项都有(这不是废话);对于 这个因式,从第二项到最后一项都有;对于
(x−x0)(x−x1)
(
x
−
x
0
)
(
x
−
x
1
)
这个因式,从第三项到最后一项都有
⋯
⋯
发现什么没有?这玩意非常像那个叫秦九韶算法的玩意。因此我们考虑预处理差商,然后利用秦九韶算法从内往外求出
P(x)
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
)
需要满足
L(xi)=yi
L
(
x
i
)
=
y
i
,因此我们考虑构造这样一个式子,使得
fi(xi)=yi
f
i
(
x
i
)
=
y
i
且其余的
fi(xj) (j≠i)
f
i
(
x
j
)
(
j
≠
i
)
都为
0
0
。
显然如果取 有
(x−xj) (j≠i)
(
x
−
x
j
)
(
j
≠
i
)
这些因式,那么就可以保证这一点。
因此取
那么我们要构造的多项式即是
那么 L(x) L ( x ) 称为 拉格朗日插值多项式, fi(x) f i ( x ) 称为 拉格朗日基本多项式(或者插值基函数)。
复杂度
显然求每一个 fi(x) f i ( x ) 是 O(n) O ( n ) 的,于是求 L(x) L ( x ) 是 O(n×n)=O(n2) O ( n × n ) = O ( n 2 ) 的。因此对于每个插值点,拉格朗日插值法的复杂度是 O(n2) O ( n 2 ) 的。
优缺点
拉格朗日插值多项式结构整齐,计算十分方便(比牛顿插值法好记多了)。但是在实际计算中,当插值点有所修改时,所对应的整个多项式就需要全部重新计算,复杂度是 O(n2) O ( n 2 ) 的(相比之下牛顿插值法只需要 O(n) O ( n ) ),因此我们就要用重心拉格朗日插值法或牛顿插值法来代替。此外,当插值点比较多的时候,拉格朗日插值多项式的次数可能会很高,因此具有数值不稳定的特点,也就是说尽管在已知的几个点取到给定的数值,但在附近却会和“实际上”的值之间有很大的偏差。这类现象也被称为龙格现象,相同地,牛顿插值法由于也是一个 n n 次的多项式,也会出现同样的问题。解决的办法是分段用较低次数的插值多项式。
重心插值法
再观察一下这个式子
我们会发现我们每个 fi(x) f i ( x ) 都重复计算了这样一个东西
于是乎,我们记
来解决这个问题,这样一来 L(x) L ( x ) 就会变成这样:
为了解决拉格朗日插值法每次加入插值点都需要重新计算 fi(x) f i ( x ) 的问题,对于每一个 fi( f i ( 我们记重心权
则
这样,求预处理 wi w i 的复杂度为 O(n2) O ( n 2 ) ,单次求值 L(x) L ( x ) 的复杂度为 O(n) O ( n ) 。
而对于动态加点的问题,我们只需要 O(n) O ( n ) 地修改 wi w i 就可以了,避免了 O(n2) O ( n 2 ) 的复杂运算,具体做法呢,就是把所有的 wi w i 除以 xi−xn+1 x i − x n + 1 得到新的 wn+1 w n + 1 ,这样一来我们的添加插值点的复杂度就变成 O(n) O ( n ) 的了。
线性插值法
所谓线性插值法,并不是指你把插值点两两直线连起来的那个插值法,而是能在保证插值点 xi=x0+i x i = x 0 + i 的情况下将复杂度优化到 O(n) O ( n ) 的重心拉格朗日插值法。
推导如下
由上面讲的重心插值法,我们有
其中的 g(x) g ( x ) 显然可以 O(n) O ( n ) 预处理,因此我们要想办法 O(n) O ( n ) 地计算
显然前面那个 ∑ ∑ 是省不掉的啊,所以我们考虑如何 O(1) O ( 1 ) 求
我们考虑把下面的 ∏ ∏ 分成 j<i j < i 与 j>i j > i 两类,那么
由于 xi=x0+i x i = x 0 + i ,那么
这样一来,我们可以预处理
那么
于是我们就可以 O(n) O ( n ) 地预处理 vi v i ,然后 O(1) O ( 1 ) 地求出 wi w i ,然后 O(n) O ( n ) 地计算
于是此时,我们的插值法就变成线性的啦(耶)
(代码与例题未完待续…
(可能会没空写…