Language interpolation and Aitken interpolation

两种插值方法(拉格朗日法和艾特肯法)

先复制一段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!(xx0)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!(xx0)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+θ(xx0)])k!(xx0)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=0kaixi,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 x00x10x20xk0x01x11x21xk1x0kx1kx2kxkk=0<=i<j<=k(xixj)=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=ikxixjxxj得到拉格朗日内插函数 L ( x ) = ∑ i = 0 k y i l i ( x ) L(x)=\sum_{i=0}^{k}y_il_i(x) L(x)=i=0kyili(x)利用拉格朗日内插函数求得内插值。

拉格朗日函数的性质

  1. 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带入基函数即可得。

  2. 拉格朗日函数本质上就是基函数根据权进行叠加。感觉这条性质到处都是(说出来只是为了说明我好像有在好好看书)

  3. 拉格朗日内插法本质上与上述直接求解线性方程组得到的函数是相同的函数。将 L ( x ) L(x) L(x)展开后不难得到这一结论。

  4. 假如原函数 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+θ(xx0)])k!(xx0)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)=xnxk1p0,1,,k1,n(x)p0,1,,k(x)xkxxnx
利用递推公式不难得出艾特肯内插法其实就是缓存了各个基函数的拉格朗日内插法,其优点在于每次更新内插多项式时,时间复杂度为 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
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值