两种插值方法(拉格朗日法和艾特肯法)
先复制一段wiki上对插值的描述:数学的数值分析领域中,内插或称插值(英语:interpolation)是一种通过已知的、离散的数据点,在范围内推求新数据点的过程或方法。
本文以多项式内插值为主,介绍拉格朗日法和艾特肯法构造内插函数。
内插原理
由Taylor展开,有
f
(
x
)
=
∑
i
=
0
i
n
f
(
f
(
i
)
(
x
0
)
(
x
−
x
0
)
i
i
!
)
f(x)=\sum_{i=0}^{inf}(f^{(i)}(x_0)\frac{(x-x_0)^i}{i!})
f(x)=∑i=0inf(f(i)(x0)i!(x−x0)i),我们将之在
i
=
k
i=k
i=k处截断,得到一个近似函数
f
k
(
x
)
=
∑
i
=
0
k
(
f
(
i
)
(
x
0
)
(
x
−
x
0
)
i
i
!
)
f_{k}(x)=\sum_{i=0}^{k}(f^{(i)}(x_0)\frac{(x-x_0)^i}{i!})
fk(x)=∑i=0k(f(i)(x0)i!(x−x0)i),误差项为
R
(
x
)
=
f
(
k
+
1
)
[
x
0
+
θ
(
x
−
x
0
)
]
)
(
x
−
x
0
)
k
+
1
k
!
R(x)=f^{(k+1)}[x_0+\theta(x-x_0)])\frac{(x-x_0)^{k+1}}{k!}
R(x)=f(k+1)[x0+θ(x−x0)])k!(x−x0)k+1,该误差项即高数中的拉格朗日余项。由高等数学我们得知当k充分大时,拉格朗日余项趋近于0,这给我们用近似函数
f
k
(
x
)
f_k(x)
fk(x)代替
f
(
x
)
f(x)
f(x)提供理论基础。在给定已知点
(
x
i
,
y
i
)
,
i
∈
[
0
,
k
]
(x_i,y_i),i\in[0,k]
(xi,yi),i∈[0,k]的基础上,我们总可以通过求解线性方程组,得到多项式各系数,进而求得内插函数。
证明如下:
以
f
(
x
)
=
∑
i
=
0
k
a
i
x
i
,
a
i
=
f
(
i
)
(
x
0
)
i
!
f(x)=\sum_{i=0}^{k}a_{i}x^{i},a_i=\frac{f^{(i)}(x_0)}{i!}
f(x)=i=0∑kaixi,ai=i!f(i)(x0)代替上述
f
k
(
x
)
f_k(x)
fk(x),带入
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi),得到关于
a
i
a_i
ai的k+1元一次方程组,由于系数矩阵
∣
x
0
0
x
0
1
…
x
0
k
x
1
0
x
1
1
…
x
1
k
x
2
0
x
2
1
…
x
2
k
⋮
⋮
⋱
x
k
0
x
k
1
…
x
k
k
∣
=
∏
0
<
=
i
<
j
<
=
k
(
x
i
−
x
j
)
≠
0
\begin{vmatrix}x_0^0 & x_0^1 & \dots &x_0^k\\ x_1^0 & x_1^1 & \dots &x_1^k\\ x_2^0 & x_ 2^1& \dots &x_2^k\\ \vdots & \vdots & \ddots\\ x_k^0 & x_ k^1& \dots &x_k^k\\ \end{vmatrix}=\prod_{0<=i<j<=k}(x_i-x_j)\quad \neq0
∣∣∣∣∣∣∣∣∣∣∣x00x10x20⋮xk0x01x11x21⋮xk1………⋱…x0kx1kx2kxkk∣∣∣∣∣∣∣∣∣∣∣=0<=i<j<=k∏(xi−xj)=0所以该线性方程组必有唯一解,可以利用克莱默法则求得对应解。
我们可以粗略估计求行列式的值的时间复杂度,LU分解法求行列式时间复杂度
O
(
K
3
)
O(K^3)
O(K3),共K个行列式求值
O
(
K
4
)
O(K^4)
O(K4)。
拉格朗日内插法
先上结论:
拉格朗日内插法的做法是取K个基函数
l
i
(
x
)
=
∏
j
=
0
,
j
≠
i
k
x
−
x
j
x
i
−
x
j
l_i(x)=\prod_{j=0,j\neq i}^{k}\quad \frac{x-x_j}{x_i-xj}
li(x)=j=0,j=i∏kxi−xjx−xj得到拉格朗日内插函数
L
(
x
)
=
∑
i
=
0
k
y
i
l
i
(
x
)
L(x)=\sum_{i=0}^{k}y_il_i(x)
L(x)=i=0∑kyili(x)利用拉格朗日内插函数求得内插值。
拉格朗日函数的性质
-
在 x = x i x=x_i x=xi处, l i ( x ) = 1 l_i(x)=1 li(x)=1,在 x = x j , j ≠ i x=x_j,j\neq i x=xj,j=i处, l i ( x ) = 0 l_i(x)=0 li(x)=0。将 x = x i x=x_i x=xi带入基函数即可得。
-
拉格朗日函数本质上就是基函数根据权进行叠加。感觉这条性质到处都是(说出来只是为了说明我好像有在好好看书)
-
拉格朗日内插法本质上与上述直接求解线性方程组得到的函数是相同的函数。将 L ( x ) L(x) L(x)展开后不难得到这一结论。
-
假如原函数 f ( x ) f(x) f(x)是 k 1 k_1 k1次多项式函数,且 k 1 < = k k_1<=k k1<=k,由拉格朗日余项公式可以得知 R ( x ) = f ( k + 1 ) [ x 0 + θ ( x − x 0 ) ] ) ( x − x 0 ) k + 1 k ! = 0 R(x)=f^{(k+1)}[x_0+\theta(x-x_0)])\frac{(x-x_0)^{k+1}}{k!}=0 R(x)=f(k+1)[x0+θ(x−x0)])k!(x−x0)k+1=0,即有拉格朗日插值误差始终为0.
一些没必要放的代码
%丁度觉得有点没必要放MATLAB代码
function main()
x=1:0.5:3;
y=2*sin(x);
a=x(1)-1;
b=x(end)+1
x1=a:0.1:b;
for i =1:length(x)
[yk(i),y1(:,i)]=languagei(x,y,1.6,x1,i);
end
yk
double(sum(yk))
2*sin(1.6)
y1=sum(y1');
figure
plot(x,y,'b*',x1,y1,'ro')
hold on
[f,fk,f1]=Aitken(x,y,1.6,x1);
double(fk)
plot(x1,f1,'g')
legend('数据点','拉格朗日插值','Aitken插值')
end
%实际使用的时候可以只保留yk返回值,
%这里使用y1是为了作图方便
%各个输入参数的含义:
%x,y为给定数据点
%x1与xk都为待插值的数据点 差别为x1为数组 xk为数据点
%i为第i个基函数
%返回值f1为数组,fk为数据点
function [fk,f1]=languagei(x,y,xk,x1,i)
syms t
l=y(i);
for k =1:length(x)
if k~=i
l=l*(t-x(k))/(x(i)-x(k));
end
end
l;
f1=subs(l,'t',x1);
fk=subs(l,'t',xk);
end
艾特肯内插法
虽然拉格朗日内插法将时间复杂度降低至
O
(
N
2
)
O(N^2)
O(N2),但是其存在的一个问题是,一旦我们需要增加数据点,其所有基函数将需要重新计算,更新时间很慢。在其基础上,我们引入了艾特肯内插法。
艾特肯内插法的核心思想是使用之前建立好的k次内插函数,线性内插得到k+1次插值多项式。其推导公式为:
p
0
,
1
,
…
,
k
,
n
(
x
)
=
1
x
n
−
x
k
∣
p
0
,
1
,
…
,
k
−
1
,
n
(
x
)
x
k
−
x
p
0
,
1
,
…
,
k
(
x
)
x
n
−
x
∣
p_{0,1,\dots,k,n}(x)=\frac{1}{x_n-x_k}\begin{vmatrix} p_{0,1,\dots,k-1,n}(x) & x_k-x \\ p_{0,1,\dots,k}(x) & x_n-x \end{vmatrix}
p0,1,…,k,n(x)=xn−xk1∣∣∣∣p0,1,…,k−1,n(x)p0,1,…,k(x)xk−xxn−x∣∣∣∣
利用递推公式不难得出艾特肯内插法其实就是缓存了各个基函数的拉格朗日内插法,其优点在于每次更新内插多项式时,时间复杂度为
O
(
N
)
O(N)
O(N)
一些没必要的代码
%用python可以通过缺省输入参数可以直观体现初次建立以及后续更新的差异
%main函数还是用的上面的调用方式
%各个输入参数的含义:
%x,y为给定数据点
%x1与xk都为待插值的数据点 差别为x1为数组 xk为数据点
%返回值f1为数组,fk为数据点
function [f,fk,f1]=Aitken(x,yaitken,xk,x1)
syms t
n=length(x)
yaitken1(1:n)=t;
for i=1:n-1
for j =i+1:n
yaitken1(j)=(yaitken(j)*(t-x(i))-yaitken(i)*(t-x(j)))/(x(i)-x(j));
end
yaitken=yaitken1;
end
f=yaitken1(n)
fk=subs(yaitken1(n),'t',xk)
f1=subs(yaitken1(n),'t',x1)
end