matlab的插值方法

1. Lagrange插值

插值是在已知数据之间寻找估计值的过程。在信号处理和图像处理中，插值极其常用。

类型很多：比如多项式插值，一、二、三维插值，样条插值等。

方法介绍：
对给定的n个插值点x1,x2,⋯,xnx1,x2,⋯,xn{x_1},{x_2}, \cdots ,{x_n}及对应的函数值y1,y2,⋯,yny1,y2,⋯,yn{y_1},{y_2}, \cdots ,{y_n}，利用构造的n-1次Lagrange插值多项式，则对插值区间内任意x的函数值y可通过下式求的：

MATLAB实现:

function y=lagrange(x0,y0,x)
ii=1:length(x0);   y=zeros(size(x));
for i=iiij=find(ii~=i);   y1=1;for j=1:length(ij),  y1=y1.*(x-x0(ij(j)));  endy=y+y1*y0(i)/prod(x0(i)-x0(ij));
end

例1：给出f(x)=ln(x)的数值表，用Lagrange计算ln(0.54)的近似值。

>> x=[0.4:0.1:0.8]；
>> y=[-0.916291,-0.693147,-0.510826,-0.356675,-0.223144];
>> lagrange(x,y,[0.54,0.55,0.78])
ans =-0.6161   -0.5978   -0.2484

2. Hermite插值

方法介绍：
不少实际问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。
已知n个插值点x1,x2,⋯,xnx1,x2,⋯,xn{x_1},{x_2}, \cdots ,{x_n}及对应的函数值y1,y2,⋯,yny1,y2,⋯,yn{y_1},{y_2}, \cdots ,{y_n}和一阶导数值y1′,y2′,⋯,y′ny1′,y2′,⋯,y′n{y_1}',{y_2}', \cdots ,{y'}_n。则对插值区间内任意x的函数值y的Hermite插值公式：

MATLAB实现:

hermite.m
function y=hermite(x0,y0,y1,x)
n=length(x0);  m=length(x);
for k=1:m    yy=0.0;for i=1:n        h=1.0;  a=0.0;for j=1:nif j~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;
end

例2:对给定数据,试构造Hermite多项式求出sin0.34的近似值。

>> x0=[0.3,0.32,0.35];
>> y0=[0.29552,0.31457,0.34290];
>> y1=[0.95534,0.94924,0.93937];
>> y=hermite(x0,y0,y1,0.34)
y =0.3335
>> sin(0.34)               ％与精确值比较
ans =0.3335

例3：

>> x=[0.3:0.005:0.35];
>>y=hermite(x0,y0,y1,x);
>> plot(x,y)
>> y2=sin(x);
>>hold on
>> plot(x,y2,'--r')

3. Runge现象

问题的提出：
根据区间[a,b]上给出的节点做插值多项式p(x)的近似值，一般总认为p(x)的次数越高则逼近f(x)的精度就越好，但事实并非如此。
反例：

f(x)=11+x2f(x)=11+x2

f(x) = \frac{1}{{1 + {x^2}}}
在区间[-5,5]上的各阶导数存在，但在此区间上取n个节点所构成的Lagrange插值多项式在全区间内并非都收敛。
例4：

x=[-5:1:5];
y=1./(1+x.^2);
x0=[-5:0.1:5];
y0=lagrange(x,y,x0);
y1=1./(1+x0.^2);%绘制图形
plot(x0,y0,'--r')%插值曲线
hold on
plot(x0,y1,'-b')%原曲线

为解决Rung问题，引入分段插值。

4.分段插值

算法分析：
所谓分段插值就是通过插值点用折线或低次曲线连接起来逼近原曲线。
MATLAB实现可调用内部函数。
interp1:
功能： 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。
格式1 yi = interp1(x,Y,xi)
%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。

例5：对于t，beta 、alpha分别有两组数据与之对应，用分段线性插值法计算当t=321, 440, 571时beta 、alpha的值。

>> temp=[300,400,500,600]';
>> beta=1000*[3.33,2.50,2.00,1.67]';
>> alpha=10000*[0.2128,0.3605,0.5324,0.7190]';
>> ti=[321,400,571]';
>> propty=interp1(temp,[beta,alpha],ti);
propty=interp1(temp,[beta,alpha],ti ,‘linear’);（默认）
>> [ti,propty]
ans =1.0e+003 *0.3210    3.1557    2.43820.4000    2.5000    3.60500.5710    1.7657    6.6489

格式2 yi = interp1(Y,xi)
%默认x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。
格式3 yi = interp1(x,Y,xi,method) %用指定的算法计算插值：
’nearest’：最近邻点插值，直接完成计算；
’linear’：线性插值（缺省方式），直接完成计算；
’spline’：三次样条函数插值。
’cubic’： 分段三次Hermite插值。
其它，如’v5cubic’ 。（三次多项式插值）
对于超出x范围的xi的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。
yi = interp1(x,Y,xi,method,’extrap’)
yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。

如何选择？

考虑方法的执行速度，占用内存，获得数据的平滑度。
’nearest’：速度最快，但数据平滑度最差，一般变化不连续。
’linear’：与上相比占内存更多，执行速度稍慢，但数据平滑度更优，数据是连续变化的。
’spline’：速度最慢，占用内存小，数据光滑，常用。
’cubic’： 速度和占内存都比线性差，但数据和一阶导数都是连续的。
例6：

>> year=1900:10:2010;
>> product=[79.995,91.972,109.711,123.203,131.669,150.697,
179.323,203.212,226.505,249.633,256.344,267.893];
>> p1995 = interp1(year,product,1995)
p1995 =252.9885
>> x = 1900:1:2010;
>> y = interp1(year,product,x,'cubic');
>> plot(year,product,'o',x,y)

例7:已知的数据点来自函数

f(x)=(x2−3x+5)e−5xsinxf(x)=(x2−3x+5)e−5xsin⁡x

f(x) = ({x^2} - 3x + 5){e^{ - 5x}}\sin x
根据生成的数据进行插值处理，得出较平滑的曲线直接生成数据。

>> x=0:.12:1;
>> y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> plot(x,y,x,y,'o')

调用interp1( )函数：

***未运行出***
x1=0:02:1;
y0=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1);
y1=interp1(x,y,x1);
y2=interp1(x,y,x1,'cubic');
y3=interp1(x,y,x1,'spline');
y4=interp1(x,y,x1,'nearest');
plot(x1,[y1',y2',y3',y4'],':',x,y,'o',x1,y0)

误差分析：

未运行出
[max(abs(y0(1:49)-y2(1:49))),max(abs(y0-y3)),max(abs(y0-y4))]

例6：对f(x)=1(1+25x2),−1≤x≤1f(x)=1(1+25x2),−1≤x≤1f(x) = \frac{1}{{(1 + 25{x^2})}}, - 1 \le x \le 1进行Lagrange插值

>> x0=-1+2*[0:10]/10;
>> y0=1./(1+25*x0.^2);
>> x=-1:.01:1;
>> y=lagrange(x0,y0,x);  % Lagrange 插值
>> ya=1./(1+25*x.^2);
>> plot(x,ya,x,y,':')

例7（基于例6数据）：分段插值

y1=interp1(x0,y0,x,'cubic');
y2=interp1(x0,y0,x,'spline');
plot(x,ya,x,y1,':',x,y2,'--')

interp2 ：
功能：二维数据内插值（表格查找）（用于图像处理和数据的可视化）
格式1 ZI = interp2(X,Y,Z,XI,YI)
%返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回NaN。
格式2 ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。
格式3 ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值：
’linear’：双线性插值算法（缺省算法）；
’nearest’：最临近插值；
’spline’：三次样条插值；
’cubic’：双三次插值。
例8：

>> years=1950:10:1990;
>> service=10:10:30;
>> wage = [150.697 199.592 187.625179.323 199.072 250.287203.212 179.092 322.767226.505 153.706 426.730249.633 120.281 599.243];
>> w = interp2(service,years,wage,15,1975)
w =190.6288

例9：由z=f(x,y)=(x2−2x)e−x2−y2−xyz=f(x,y)=(x2−2x)e−x2−y2−xyz = f(x,y) = ({x^2} - 2x){e^{ - {x^2} - {y^2} - xy}}可计算出一些较稀疏的网格数据，对整个函数曲面进行各种插值拟合，并比较拟合结果，绘制已知数据的网格图：

>> [x,y]=meshgrid(-3:.6:3,-2:.4:2);%.6、。4代表x,y轴单位间距
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,z)%三维着色表面图
axis([-3,3,-2,2,-0.7,1.5])%标注输出的图线的最大值最小值。

选较密的插值点，用默认的线性插值算法进行插值

>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=interp2(x,y,z,x1,y1);
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])

立方和样条插值：

>> z1=interp2(x,y,z,x1,y1,'cubic');
>> z2=interp2(x,y,z,x1,y1,'spline');
>> surf(x1,y1,z1)
>>axis([-3,3,-2,2,-0.7,1.5])
>> figure
>>surf(x1,y1,z2)
>>axis([-3,3,-2,2,-0.7,1.5])

算法误差的比较:

> z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>>  surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.08])
>>  figure;surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.025])%abs:数值的绝对值和复数的幅值

二维一般分布数据的插值
功能:可对非网格数据进行插值
**格式:**z=griddata(x0,y0,z0,x,y,’method’)
’ v4 ’:MATLAB4.0提供的插值算法，公认效果较好；
’linear’：双线性插值算法（缺省算法）；
’nearest’：最临近插值；
’spline’：三次样条插值；
’cubic’：双三次插值。
例10：z=f(x,y)=(x2−2x)e−x2−y2−xyz=f(x,y)=(x2−2x)e−x2−y2−xyz = f(x,y) = ({x^2} - 2x){e^{ - {x^2} - {y^2} - xy}}
在x为［-3,3］，y为［-2,2］矩形区域随机选择一组坐标，用’ v4 ’与’cubic’插值法进行处理，并对误差进行比较。

>> x=-3+6*rand(200,1);
>>y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'cubic');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
>> z2=griddata(x,y,z,x1,y1,'v4');
>> figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])

误差分析

>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>>  surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])
>>  figure;surf(x1,y1,abs(z0-z2)),axis([-3,3,-2,2,0,0.15])

例11:z=f(x,y)=(x2−2x)e−x2−y2−xyz=f(x,y)=(x2−2x)e−x2−y2−xyz = f(x,y) = ({x^2} - 2x){e^{ - {x^2} - {y^2} - xy}}在x为［3,3],y为［-2,2］矩形区域随机选择一组坐标中，对分布不均匀数据，进行插值分析。

>> x=-3+6*rand(200,1); %x=rand(m,n)产生m行n列的位于（0，1）区间的随机数
>>y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);   % 生成已知数据
>> plot(x,y,'x')   % 样本点的二维分布
>> figure, plot3(x,y,z,'x'), axis([-3,3,-2,2,-0.7,1.5])，grid %grid打开网格

去除在(-1,-1/2)点为圆心，以0.5为半径的圆内的点

>> x=-3+6*rand(200,1); y=-2+4*rand(200,1);  % 重新生成样本点
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> ii=find((x+1).^2+(y+0.5).^2>0.5^2);  % 找出满足条件的点坐标
>> x=x(ii); y=y(ii); z=z(ii); plot(x,y,'x')
>> t=[0:.1:2*pi,2*pi]; x0=-1+0.5*cos(t); y0=-0.5+0.5*sin(t);
>> line(x0,y0)   % 在图形上叠印该圆，可见，圆内样本点均已剔除

用新样本点拟合出曲面

>> [x1,y1]=meshgrid(-3:.2:3, -2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'v4');
>> surf(x1,y1,z1), axis([-3,3,-2,2,-0.7,1.5])

误差分析

>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z0-z1)), axis([-3,3,-2,2,0,0.1])
>> contour(x1,y1,abs(z0-z1),30); hold on, plot(x,y,‘x’); line(x0,y0)  ％误差的二维等高线图

interp3：
三维网格生成用meshgrid( )函数,
调用格式: [x,y,z]=meshgrid(x1,y1,z1)
其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。
griddata3( ) 三维非网格形式的插值拟合
interpn：
n维网格生成用ndgrid( )函数,
调用格式: [x1,x2,…,xn]=ndgrid[v1,v2,…,vn]
griddatan( ) n维非网格形式的插值拟合interp3 ( )、 interpn( )调用格式同interp2( )函数一致；
griddata3( )、 griddatan( )调用格式同griddata( )函数一致。
例12：V(x,y,z)=ex2z+y2x+z2ycos(x2yz+z2yx)V(x,y,z)=ex2z+y2x+z2ycos⁡(x2yz+z2yx)V(x,y,z) = {e^{{x^2}z + {y^2}x + {z^2}y}}\cos ({x^2}yz + {z^2}yx)
通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。

>> [x,y,z]=meshgrid(-1:0.2:1);  [x0,y0,z0]=meshgrid(-1:0.05:1);
>> V=exp(x.^2.*z+y.^2.*x+z.^2.*y).*cos(x.^2.*y.*z+z.^2.*y.*x);>> V0=exp(x0.^2.*z0+y0.^2.*x0+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0);>> V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:))
ans =0.0419

5. 样条插值的MATLAB表示

spline:
功能:三次样条数据插值
**格式:**yy = spline(x,y,xx)
例13:对离散分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算：

>>x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; >> y = exp(x).*sin(x);>>xx = 0:.25:20;>>yy = spline(x,y,xx);>>plot(x,y,'o',xx,yy)