【数值分析】插值方法
几种易于实现的插值方法
插值,即构造 f ( x ) f(x) f(x)的近似函数 p ( x ) p(x) p(x),若 p ( x ) p(x) p(x)取离散数据,称之为 f ( x ) f(x) f(x)的插值函数,使得在某些点的值或者加上导数值相等
根据插值函数类型分类:代数插值(代数多项式),样条插值
这里以王能超《数值分析简明教程》为基础,对泰勒插值、拉格朗日插值和牛顿插值等基础插值方式作介绍
泰勒插值
插值公式
取 f ( x ) f(x) f(x)在 x 0 x_0 x0处的泰勒展开(泰勒多项式) p n ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + ⋯ + f ( n ) n ! ( x − x 0 ) n p_n(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\cdots+\frac{f^{(n)}}{n!}(x-x_0)^n pn(x)=f(x0)+f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+⋯+n!f(n)(x−x0)n在某点 x 0 x_0 x0逼近 f ( x ) f(x) f(x),值相等,各阶导数值相等
余项分析
插值函数
p
n
(
x
)
p_n(x)
pn(x)的截断误差或称插值余项
R
(
x
)
=
f
(
x
)
−
p
n
(
x
)
R(x)=f(x)-p_n(x)
R(x)=f(x)−pn(x)
泰勒余项定理:含有
x
0
x_0
x0区间
[
a
,
b
]
[a,b]
[a,b]内有直到
n
+
1
n+1
n+1阶导数,区间内
f
(
x
)
−
p
n
(
x
)
=
f
(
n
+
1
)
(
ε
)
(
n
+
1
)
!
(
x
−
x
0
)
n
+
1
,
ε
∈
[
a
,
b
]
f(x)-p_n(x)=\frac{f^{(n+1)}(\varepsilon)}{(n+1)!}(x-x_0)^{n+1}, \varepsilon\in [a,b]
f(x)−pn(x)=(n+1)!f(n+1)(ε)(x−x0)n+1,ε∈[a,b]
拉格朗日插值
问题的描述
作
n
n
n次多项式
p
n
(
x
)
p_n(x)
pn(x),使得在各插值节点(
n
+
1
n+1
n+1个互不相同,
x
i
,
i
=
0
,
1
,
2
,
⋯
,
n
x_i,i=0,1,2,\cdots,n
xi,i=0,1,2,⋯,n)上的值与函数
f
(
x
)
f(x)
f(x)节点值相等
即构造
n
n
n次多项式
p
n
(
x
)
=
a
0
+
a
1
x
+
a
2
x
0
2
+
⋯
+
a
n
x
n
p_n(x)=a_0+a_1x+a_2x_0^2+\cdots+a_nx^n
pn(x)=a0+a1x+a2x02+⋯+anxn使得
p
n
(
x
i
)
=
y
i
=
f
(
x
i
)
,
i
=
0
,
1
,
⋯
,
n
p_n(x_i)=y_i=f(x_i),i=0,1,\cdots,n
pn(xi)=yi=f(xi),i=0,1,⋯,n
拉格朗日多插值多项式的系数的项数与插值点数相等,并对于给定的互不相同节点,插值多项式
p
n
(
x
)
p_n(x)
pn(x)一定存在:
求解关于系数
a
0
,
a
1
,
⋯
,
a
n
a_0,a_1,\cdots,a_n
a0,a1,⋯,an的线性方程组
{
a
0
+
a
1
x
0
+
a
2
x
0
2
+
⋯
+
a
n
x
0
n
=
y
0
a
0
+
a
1
x
1
+
a
2
x
1
2
+
⋯
+
a
n
x
1
n
=
y
1
⋯
a
0
+
a
1
x
n
+
a
2
x
n
2
+
⋯
+
a
n
x
n
n
=
y
n
\begin{cases}a_0+a_1x_0+a_2x_0^2+\cdots+a_nx_0^n=y_0\\a_0+a_1x_1+a_2x_1^2+\cdots+a_nx_1^n=y_1\\\cdots\\a_0+a_1x_n+a_2x_n^2+\cdots+a_nx_n^n=y_n\end{cases}
⎩
⎨
⎧a0+a1x0+a2x02+⋯+anx0n=y0a0+a1x1+a2x12+⋯+anx1n=y1⋯a0+a1xn+a2xn2+⋯+anxnn=yn
注意多项式的次数其实是度数要注意零次幂常数项的存在
其系数行列式为范德蒙(Vandermonde)行列式,性质即证
插值函数求解方法
直接求解上面的线性方程组过于复杂,尝试找出显式规律
构造的思路
利用插值节点处的值,每个节点
x
i
x_i
xi构造对应的插值基函数
l
i
(
x
)
l_i(x)
li(x)(关于
x
x
x,与系数组成插值函数),并将其与节点的值
y
i
y_i
yi(系数)分离,更好的体现出“基函数”性
核心思想:每个节点插值函数只在当前节点有贡献,其他节点为零,保证当前节点的值对应的同时避免对以同样方式构造的其他插值基函数负责的节点的值的干扰。注意到节点之间为各插值函数组合。
上面的构造思想道出了插值函数应满足的性质
拉格朗日插值公式的形式
根据构造思路的讨论,
y
=
f
(
x
)
y=f(x)
y=f(x)的插值函数
p
n
(
x
)
=
∑
i
=
0
n
y
i
l
i
(
x
)
p_n(x)=\sum_{i=0}^ny_il_i(x)
pn(x)=∑i=0nyili(x)
i
i
i节点
x
i
x_i
xi对应的插值基函数
l
i
(
x
)
l_i(x)
li(x)应满足:
l
i
(
x
j
)
=
δ
i
j
=
{
1
,
j
=
i
0
,
j
≠
i
l_i(x_j)=\delta_{ij}=\begin{cases}1,j=i\\0,j\ne i\end{cases}
li(xj)=δij={1,j=i0,j=i
给出其构造形式,考虑分式,除当前节点,节点均作为零点添加进式子,当前节点与其他节点之差作为分母,使当
i
=
j
i=j
i=j时分子和分母相等,且分母不为零
l
i
(
x
)
=
∏
j
=
0
,
j
≠
k
n
x
−
x
j
x
i
−
x
j
l_i(x)=\prod_{j=0,j\ne k}^n\frac{x-x_j}{x_i-x_j}
li(x)=j=0,j=k∏nxi−xjx−xj
i和j仅代表当前处理节点不同
按照此方法构造的插值函数即拉格朗日插值公式
p
n
(
x
)
=
∑
i
=
0
n
(
∏
j
=
0
,
j
≠
k
n
x
−
x
j
x
i
−
x
j
)
y
i
p_n(x)=\sum_{i=0}^n(\prod_{j=0,j\ne k}^n\frac{x-x_j}{x_i-x_j})y_i
pn(x)=i=0∑n(j=0,j=k∏nxi−xjx−xj)yi
拉格朗日余项定理
插值节点
x
0
,
x
1
,
⋯
,
x
n
x_0,x_1,\cdots,x_n
x0,x1,⋯,xn在
[
a
,
b
]
[a,b]
[a,b]内,
f
(
x
)
f(x)
f(x)在
[
a
,
b
]
[a,b]
[a,b]有连续的直到
n
+
1
n+1
n+1阶导数,在区间内有
f
(
x
)
−
p
n
(
x
)
=
f
(
n
+
1
)
(
x
)
(
n
+
1
)
!
∏
i
=
0
n
(
x
−
x
i
)
f(x)-p_n(x)=\frac{f^{(n+1)}(x)}{(n+1)!}\prod_{i=0}^n(x-x_i)
f(x)−pn(x)=(n+1)!f(n+1)(x)i=0∏n(x−xi)
推导过程设包含了当前项
p
n
(
t
)
p_n(t)
pn(t)和
∏
i
=
0
n
(
x
−
x
i
)
\prod_{i=0}^n(x-x_i)
∏i=0n(x−xi)形式余项形式的插值表示函数,并使得辅助函数在差值节点再加上当前点
x
x
x时值与
f
(
x
)
f(x)
f(x)相等,使用罗尔定理证明该形式合理,属于倒推
反映了若插值点落在区间 [ a , b ] [a,b] [a,b]内,属于内插;若落在区间外,属于外推,不适用余项定理,外推结果不可靠
拓展和不足
增加插值点的做法
拉格朗日插值增加一个插值点,需要重新计算公式中所有的系数(插值基函数),引出误差的事后估计法
事后估计法
从简单的两点线性插值引入
选取点
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)和点
(
x
1
,
y
1
)
(x_1,y_1)
(x1,y1)作线性插值,得到插值函数
p
1
(
x
)
p_1(x)
p1(x);选取点
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)和点
(
x
2
,
y
2
)
(x_2,y_2)
(x2,y2)作线性插值,得到插值函数
p
2
(
x
)
p_2(x)
p2(x)(注意不同于此前的记号,p的下标仅代表两个不同的插值函数;之所以这样记是使用
y
1
y_1
y1和
y
2
y_2
y2表示插值函数而不是
x
1
x_1
x1和
x
2
x_2
x2的值显得并不直观,但
y
y
y仍然代表
f
(
x
)
f(x)
f(x))
根据余项定理的估计
y
−
p
1
(
x
)
≈
f
’’
(
ε
)
2
(
x
−
x
0
)
(
x
−
x
1
)
y-p_1(x)\approx\frac{f’’(\varepsilon)}{2}(x-x_0)(x-x_1)
y−p1(x)≈2f’’(ε)(x−x0)(x−x1)
y
−
p
2
(
x
)
≈
f
’’
(
ε
)
2
(
x
−
x
0
)
(
x
−
x
2
)
y-p_2(x)\approx\frac{f’’(\varepsilon)}{2}(x-x_0)(x-x_2)
y−p2(x)≈2f’’(ε)(x−x0)(x−x2)
两式相比得
y
−
p
1
(
x
)
y
−
p
2
(
x
)
≈
x
−
x
1
x
−
x
2
\frac{y-p_1(x)}{y-p_2(x)}\approx\frac{x-x_1}{x-x_2}
y−p2(x)y−p1(x)≈x−x2x−x1,这个结果容易与斜率的含义混淆,因此要记住刚刚约定的下标的含义,但是类比直线的两点式用基函数的方法尝试将
y
y
y表示出来
y
=
x
−
x
2
x
1
−
x
2
p
1
(
x
)
+
x
−
x
1
x
−
x
1
p
2
(
x
)
y=\frac{x-x_2}{x_1-x_2}p_1(x)+\frac{x-x_1}{x-x_1}p_2(x)
y=x1−x2x−x2p1(x)+x−x1x−x1p2(x)
可以看到由两个线性插值得到了新的插值公式,其含义在于
y
−
p
1
(
x
)
=
x
−
x
1
x
2
−
x
1
(
p
2
(
x
)
−
p
1
(
x
)
)
y-p_1(x)=\frac{x-x_1}{x_2-x_1}(p_2(x)-p_1(x))
y−p1(x)=x2−x1x−x1(p2(x)−p1(x))
插值结果
p
1
(
x
)
p_1(x)
p1(x)的误差
y
−
p
1
(
x
)
y-p_1(x)
y−p1(x)可以由两个结果的偏差给出,经过实例验证后还发现这种添加一个新点计算新的插值函数估计出来的偏差估计出的结果用来修正旧点的插值结果与高一次拉格朗日插值的结果相同,引出下面的埃特金算法
埃特金(Aitken)算法
为了描述递推形式的埃特金算法引入专门的记号,
f
k
(
x
i
)
f_k(x_i)
fk(xi)表示在使用
x
0
,
x
1
,
⋯
,
x
k
−
1
x_0,x_1,\cdots,x_{k-1}
x0,x1,⋯,xk−1进行
k
−
1
k-1
k−1次拉格朗日插值的基础上,加上新的插值节点
x
i
(
i
≥
k
)
x_i(i\ge k)
xi(i≥k)所得的插值结果,在下面的逐阶递推表达时
i
i
i就取
k
k
k,
f
k
(
x
k
)
f_k(x_k)
fk(xk)就是
k
k
k次拉格朗日插值多项式的结果,但仍然是关于
x
x
x的函数,不代表取值
从事后估计法的两个线性插值得到的新插值公式其实就是抛物插值的结果,归纳的有(按照基函数的构造去记忆)
f
k
(
x
i
)
=
x
−
x
i
x
k
−
1
−
x
i
f
k
−
1
(
x
k
−
1
)
+
x
−
x
k
−
1
x
i
−
x
k
−
1
f
k
−
1
(
x
i
)
f_k(x_i)=\frac{x-x_i}{x_{k-1}-x_i}f_{k-1}(x_{k-1})+\frac{x-x_{k-1}}{x_i-x_{k-1}}f_{k-1}(x_i)
fk(xi)=xk−1−xix−xifk−1(xk−1)+xi−xk−1x−xk−1fk−1(xi)
含义上不是由
k
−
1
k-1
k−1次插值新添加的
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)点得到
k
k
k次插值下
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)的结果,要从递推的角度去理解,0次插值就是
f
(
x
0
)
f(x_0)
f(x0)
因此无法给出显示形式,需要建立具有承袭性的显示的插值公式
牛顿插值
想要显示地承袭性插值公式,在
k
−
1
k-1
k−1次插值的基础上,就只能根据已知插值节点去计算新的插值节点
x
i
x_i
xi带来的项的系数(借用插值基函数的思想,不影响其它插值节点
x
0
,
x
1
,
⋯
,
x
k
−
1
x_0,x_1,\cdots,x_{k-1}
x0,x1,⋯,xk−1,故应是
c
∏
j
=
0
,
j
≠
i
k
−
1
(
x
−
x
j
)
c\prod_{j=0,j\ne i}^{k-1}(x-x_j)
c∏j=0,j=ik−1(x−xj)的形式),求解出新的系数值
对比拉格朗日插值公式的插值基函数设法,其实牛顿插值是拉格朗日插值的另一表现形式,直接在上一次插值基础上递推求出当前次的拉格朗日插值(埃特金算法需要两个上一次插值且非标准化仅能由递推给出)
差商
为了表示方便引入差商的概念(P24给出了不用差商时的形式),这里的常变量(指不含
x
x
x)
n
n
n阶差商就对应
n
n
n次插值的
n
n
n次项
∏
i
=
0
n
(
x
−
x
i
)
\prod_{i=0}^n(x-x_i)
∏i=0n(x−xi)项的系数,而其余次项具有承袭性与前面得到的结果相同(引入
n
n
n阶差商时
n
n
n次项对应的式就永久成立)
一阶差商:
f
(
x
0
,
x
1
)
=
f
(
x
1
)
−
f
(
x
0
)
x
1
−
x
0
f(x_0,x_1)=\frac{f(x_1)-f(x_0)}{x_1-x_0}
f(x0,x1)=x1−x0f(x1)−f(x0)
二阶差商:
f
(
x
0
,
x
1
,
x
2
)
=
f
(
x
1
,
x
2
)
−
f
(
x
0
,
x
1
)
x
2
−
x
0
f(x_0,x_1,x_2)=\frac{f(x_1,x_2)-f(x_0,x_1)}{x_2-x_0}
f(x0,x1,x2)=x2−x0f(x1,x2)−f(x0,x1)
n
n
n阶差商:
f
(
x
0
,
x
1
,
⋯
,
x
n
)
=
f
(
x
1
,
x
2
,
⋯
,
x
n
)
−
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
)
x
n
−
x
0
f(x_0,x_1,\cdots,x_n)=\frac{f(x_1,x_2,\cdots,x_n)-f(x_0,x_1,\cdots,x_{n-1})}{x_n-x_0}
f(x0,x1,⋯,xn)=xn−x0f(x1,x2,⋯,xn)−f(x0,x1,⋯,xn−1)
可以看到差商也使用了类似递推的定义,补充零阶差商为
f
(
x
i
)
f(x_i)
f(xi)
将差商展开成节点的形式可得到差商的对称性,任意交换节点的次序不对差商值产生影响,差商与节点的排列顺序无关
f ( x 0 , x 1 , ⋯ , x n ) = ∑ i = 0 n f ( x i ) ∏ j = 0 , j ≠ i n ( x i − x j ) f(x_0,x_1,\cdots,x_n)=\sum_{i=0}^n\frac{f(x_i)}{\prod_{j=0,j\ne i}^n(x_i-x_j)} f(x0,x1,⋯,xn)=i=0∑n∏j=0,j=in(xi−xj)f(xi)
差商形式的插值公式(牛顿插值公式)
逆推差商定义
推导过程
根据差商的定义式尝试逆推,需要把新引进的变量
x
x
x当作此前的常数
x
i
x_i
xi看,我们的目标就是求出
f
(
x
)
f(x)
f(x)
f
(
x
)
=
f
(
x
0
)
+
f
(
x
0
,
x
)
(
x
−
x
0
)
f(x)=f(x_0)+f(x_0,x)(x-x_0)
f(x)=f(x0)+f(x0,x)(x−x0)
继续利用高阶差商项逆推下去,这时要引进的新变量不是
x
x
x了,先参照差商定义式对节点排列
f
(
x
0
,
x
)
=
f
(
x
1
,
x
0
)
+
f
(
x
1
,
x
0
,
x
)
(
x
−
x
1
)
f(x_0,x)=f(x_1,x_0)+f(x_1,x_0,x)(x-x_1)
f(x0,x)=f(x1,x0)+f(x1,x0,x)(x−x1)
根据对称性写成
f
(
x
0
,
x
)
=
f
(
x
0
,
x
1
)
+
f
(
x
0
,
x
1
,
x
)
(
x
−
x
1
)
f(x_0,x)=f(x_0,x_1)+f(x_0,x_1,x)(x-x_1)
f(x0,x)=f(x0,x1)+f(x0,x1,x)(x−x1)
继续逆推
f
(
x
0
,
x
1
,
x
)
=
f
(
x
2
,
x
0
,
x
1
)
+
f
(
x
2
,
x
0
,
x
1
,
x
)
(
x
−
x
2
)
f(x_0,x_1,x)=f(x_2,x_0,x_1)+f(x_2,x_0,x_1,x)(x-x_2)
f(x0,x1,x)=f(x2,x0,x1)+f(x2,x0,x1,x)(x−x2)
即
f
(
x
0
,
x
1
,
x
)
=
f
(
x
0
,
x
1
,
x
2
)
+
f
(
x
0
,
x
1
,
x
2
,
x
)
(
x
−
x
2
)
f(x_0,x_1,x)=f(x_0,x_1,x_2)+f(x_0,x_1,x_2,x)(x-x_2)
f(x0,x1,x)=f(x0,x1,x2)+f(x0,x1,x2,x)(x−x2)
总结规律:展开含有
x
x
x的
n
−
1
n-1
n−1阶差商
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
)
f(x_0,x_1,\cdots,x_{n-1},x)
f(x0,x1,⋯,xn−1,x),得到在原有常量节点上加上新节点
x
n
x_n
xn的常数差商加上在新的常数差商上引入变量
x
x
x的差商与
x
−
x
n
x-x_n
x−xn之积
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
)
=
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
n
)
+
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
n
,
x
)
(
x
−
x
n
)
f(x_0,x_1,\cdots,x_{n-1},x)=f(x_0,x_1,\cdots,x_{n-1},x_n)+f(x_0,x_1,\cdots,x_{n-1},x_n,x)(x-x_n)
f(x0,x1,⋯,xn−1,x)=f(x0,x1,⋯,xn−1,xn)+f(x0,x1,⋯,xn−1,xn,x)(x−xn)
结论
同样地递推得到插值公式,为了得到公式形式的一致性,遵循特定的排列顺序,运用对称性进行验证
f
(
x
)
=
f
(
x
0
)
+
f
(
x
0
,
x
)
(
x
−
x
0
)
f(x)=f(x_0)+f(x_0,x)(x-x_0)
f(x)=f(x0)+f(x0,x)(x−x0)
f
(
x
0
,
x
)
=
f
(
x
0
,
x
1
)
+
f
(
x
0
,
x
1
,
x
)
(
x
−
x
1
)
f(x_0,x)=f(x_0,x_1)+f(x_0,x_1,x)(x-x_1)
f(x0,x)=f(x0,x1)+f(x0,x1,x)(x−x1)
f
(
x
0
,
x
1
,
x
)
=
f
(
x
0
,
x
1
,
x
2
)
+
f
(
x
0
,
x
1
,
x
2
,
x
)
(
x
−
x
2
)
f(x_0,x_1,x)=f(x_0,x_1,x_2)+f(x_0,x_1,x_2,x)(x-x_2)
f(x0,x1,x)=f(x0,x1,x2)+f(x0,x1,x2,x)(x−x2)
f
(
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
)
=
f
(
x
0
,
x
1
,
⋯
,
x
n
)
+
f
(
x
0
,
x
1
,
⋯
,
x
n
,
x
)
(
x
−
x
n
)
f(x_0,x_1,\cdots,x_{n-1},x)=f(x_0,x_1,\cdots,x_n)+f(x_0,x_1,\cdots,x_n,x)(x-x_n)
f(x0,x1,⋯,xn−1,x)=f(x0,x1,⋯,xn)+f(x0,x1,⋯,xn,x)(x−xn)
牛顿插值公式
逐层代入含有
x
x
x的差商使得其只在末项作为余项
R
(
x
)
R(x)
R(x)存在
得到
p
n
(
x
)
p_n(x)
pn(x)即牛顿插值公式
f
(
x
)
=
p
n
(
x
)
+
R
(
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
)
+
⋯
+
f
(
x
0
,
x
1
,
⋯
,
x
n
)
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
−
1
)
+
f
(
x
0
,
x
1
,
⋯
,
x
n
,
x
)
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
)
⏟
R
(
x
)
\begin{aligned}f(x)=&p_n(x)+R(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})\\&+\begin{matrix}\underbrace{f(x_0,x_1,\cdots,x_n,x)(x-x_0)(x-x_1)\cdots(x-x_n)}\\R(x)\end{matrix}\end{aligned}
f(x)==pn(x)+R(x)f(x0)+f(x0,x1)(x−x0)+f(x0,x1,x2)(x−x0)(x−x1)+⋯+f(x0,x1,⋯,xn)(x−x0)(x−x1)⋯(x−xn−1)+
f(x0,x1,⋯,xn,x)(x−x0)(x−x1)⋯(x−xn)R(x)
p
n
(
x
)
p_n(x)
pn(x)的每一项的系数其实是由其下一项在代入时“下放”得来,因此系数是
n
n
n阶差商但该项只有
n
−
1
n-1
n−1次
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)