实验三插值.docx

1、实验三插值试验三Matlab插值与拟合一. 实验目的了解插值的基本内容和原理；掌握用matlab求解插值问题,包括一维插值和二维插值的各种常用方法；二. 实验原理和方法插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根据观测点的值，构造一个比较简单的函数y=(x)，使函数在观测点的值等于已知的数值或导数值，进而用简单函数y=(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数(x)，办

2、法是很多的。(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数；函数类的不同，自然地有不同的逼近效果。一维插值已知个数据节点:构造一个(相对简单)函数(称为插值函数),通过全部结点即（j0,1,n）再用计算插值，即数学上插值方法非常多，只介绍几种常用方法：（1）拉格朗日(Lagrange)插值已知函数f(x)（称为被插值函数）在n+1个点x0,x1,xn处的函数值为 y0,y1,yn 。求一n次多项式函数Pn(x)（称为插值函数），使其满足： Pn(xi)=yi,i=0,1,n.解决此问题的拉格朗日插值多项式公式如下：其中为次

3、多项式：称为拉格朗日插值基函数,可以验证该多项式通过所给数据点。特别地:两点一次(线性)插值多项式为：三点二次(抛物)插值多项式为:matlab没有Lagrange插值的函数，需要自己编写。Lagrange插值插值的最大优点是函数具有很好的解析性质,无穷次可微。缺点是：（1）可能出现严重的振荡现象，在头尾附近会出现一些不好的振荡现象（龙格现象）（2）多项式的系数对数据非常敏感。导致Lagrange插值如上缺点的原因在于：多项式次数过高不一定是好事。而克服这一缺点的办法之一就是采用分段低次插值，其中最简单的是分段线性插值。（2）分段线性插值分段线性插值就是在每个子区间上做通过两个端点的直线。该方

4、法具有良好的收敛性，即，其中为被插值函数，可以通过增加插值节点的个数来控制插值误差，是一种较实用的插值方法，Matlab有分段线性插值函数，但总体光滑程度不够。（3）分段三次埃尔米特（Hermite）插值对插值函数的要求可以根据需要进行改变，由此产生了不同的插值方法，分段三次埃尔米特插值就是其一。在插值问题中，如果除了指定在插值节点的函数值以外，同时指定插值节点处的导数值，这时插值问题变为：给定一组观察数据及，要求一个分段多项式函数，满足：这就相当于所求分段多项式在每个小段上满足四个条件：可以确定四个待定参数，三次多项式恰好有四个系数，所以可以考虑用三次多项式作为插值函数，这就是分段三次Her

5、mite插值。（4）分段三次样条插值在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则曲线光滑程度越好。分段线性插值具有零阶光滑性，也就是不光滑；分段三次埃尔米特插值具有一阶光滑性，但是这样的光滑程度在工程设计和机械加工等实际应用中是远远不够的。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。什么是样条？所谓样条，就是一个细的、可弯曲的木制或塑料条，在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线(放样)所用的工具。1946年,Schoenberg将样条引入数学,即所谓的样条函数。三次样条函数本

6、质上是一段一段的三次多项式拼合而成的曲线，在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的。满足以上条件的分段函数称为三次样条插值函数。可以看出，三次样条插值是分段三次多项式插值，但是不仅要求过所给数据节点，而且要求在所给节点处的二阶导数也是连续的，从而构造出来的插值函数比较光滑。在n个小区间构造S(x)，共有n个三次多项式，每个多项式有4个系数需要确定，总共需确定4n个参数才能给出S(x)（1）在n+1个节点上要求S(x)过点：， 该条件给出n+1个方程。（2）在n-1个内部节点上要求连续，一阶导数和二阶导数连续：（），该条件给出3（n-1）个方程。 由（1）和（2）总共得到4n-2个

7、方程，还差两个。为此常用的方法是在两个边界节点和附加要求，这就是所谓的边界条件。根据实际问题的不同，三次样条插值常用到下列三类边界条件。（1）边界条件：。即指定两个边界节点的一阶导数值为给定值（2）边界条件：。即指定两个边界节点的二阶导数值为给定值特别地，当和都为零时，称为自然边界条件。(3)周期性边界条件：总而言之，理论上三次样条插值是确定的。matlab 插值命令命令1 interp1功能 一维数据插值。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图如下图图 数据点与插值点关系示意图格式：（1）yi = inte

8、rp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。（2）yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。（3）yi = interp1(x,Y,xi,method) %用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear： 线性插值（缺省方式），直接完成计算；spline：三次样条函数插值，使用命令spline的默认的边界条件。即要求第一个点的三次导数和第

9、二点的三次导数一样；最后一个点的三次导数和倒数第二个点一样。pchip：分段三次Hermite插值。对于该方法，命令interp1调用函数pchip，用于对向量x与y执行分段三次内插值。该方法保留单调性与数据的外形； cubic：与pchip操作相同；v5cubic：在MATLAB 5.0中的三次插值。对于超出x范围的xi的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。（4）yi = interp1(x,Y,xi,method,extrap) %对于超出x范围的xi中的分量将执行特殊的外插值法

10、extrap。(5)yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。（6）pp= interp1(x,Y,method,pp) %使用指定方法生成pp形式的分片多项式，可以使用格式（3）中所列举的除v5cubic外任何一种方法。PP是一种结构体类型的数据，pp中存储了各分段三次样条插值多项式的系数和其他相关的信息。pp.breaks 插值节点pp.coefs 插值多项式系数pp.pieces 多项式个数 pp.order 每个多项式系数的个数pp.dim 插值维数欲求该分片多项式在xi处

11、的值，可使用yi=ppval(pp,xi)来计算。ppval(pp,xi)等价于interp1(x,Y,xi,method,extrap).例：已知标准正态分布函数值。利用分段线性插值计算x=2.34 2.35;y=0.99036 0.99061;x0=2.3456789;y0=interp1(x,y,2.345678,linear)最后结果为：y0=0.9905例：x = 0:10; y=sin(x); xi=0:.25:10; yi=interp1(x,y,xi); plot(x,y,o,xi,yi)例：x = 0:10; y = x.*sin(x); xx = 0:.25:10; yy

12、= interp1(x,y,xx); plot(x,y,kd,xx,yy)例year=1900:10:2010; %年份从1900到2010，每间隔为10年product=75.995 91.972 105.711 123.203 131.669 150.697 .179.323 203.212 226.505 249.633 256.344 267.893 ; %产量 p1995=interp1(year,product,1995) %插值计算1995年产量，采用%默认的线性差值方法。x = 1900:1:2010; %1900到2010年，每间隔1年%采用分段三次Hermite插值y =

13、interp1(year,product,x,pchip); plot(year,product,o,x,y)插值结果为：p1995=252.9885例：x=0:.2:pi;y= sin(x);pp=interp1(x,y,cubic,pp);xi=0:.1:pi;yi=ppval(pp,xi);plot(x,y,ko), hold on, plot(xi,yi,r:), hold off例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。解：hours=1:12;temps=5

14、 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); % (直接输出数据将是很多的)plot(hours,temps,+,h,t,r:) %作图xlabel(Hour),ylabel(Degrees Celsius)命令2 spline功能：三次样条数据插值格式（1）：pp=spline(x,y)说明：返回分段多项式形式的三次样条插值函数供函数ppval使用，Pp为x必须是一个向量.1.若y也是一个向量，而且length(x)=length(y)则所求的样条函数S（x(i)）=y(i),此时用非扭

15、结(not-a-knot)条件,即一个点的三次导数和第二点的三次导数一样，最后一个点的三次导数和倒数第二个点一样.2.如果向量y的元素个数比x多两个，此时Y的第一个和最后一个元素的值用作所求三次样条函数的S(x)的m边界条件：S(x)=Y(2:end-1)dS(min(x)=Y(1) % Y(1)和Y(end)用作导数值。dS(max(x)=Y(end)3.如果y是矩阵，则以y的每一行和x配对构造三次样条插值函数。注意：interp1是按列做插值，而 spline是按行做插值格式（2）：yy=spline(x,y,xx)说明：返回由向量x和y确定的三次样条函数在xx处的函数值yy。x和y的要求

16、同格式（1）。例：对离散地分布在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)插值图形结果为下图。例：x = -4:4;y = 0 .15 1.12 2.36 2.36 1.46 .49 .06 0;cs = spline(x,0 y 0); % 数组0 y 0所含分两个数比x多两个，% 使用指定一阶导数的边界条件xx = linspace(-4,4,101);plot(x,y

17、,o,xx,ppval(cs,xx),-);例：x = 0:.25:1;Y = sin(x); cos(x);% 数据Y为含两行的矩阵.xx = 0:.1:1;YY = spline(x,Y,xx);% YY的第一行为正弦数据插在xx上的插值 % YY的第二行为余弦数据插在xx上的插值 plot(x,Y(1,:),o,xx,YY(1,:),-); hold on;plot(x,Y(2,:),o,xx,YY(2,:),:); hold off;命令3 csape功能：使用各种边界条件的三次样条插值，返回PP形式的样条函数。基本使用格式：（1）pp=csape(x,y)；%使用默认的边界条件，即n

18、ot-a-knot条件。（2）pp=csape(x,y,conds,valconds)说明：conds指定边界类型，有如下几种：complete : 边界条件,即给定边界一阶导数.一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。也即 not-a-knot：非扭结条件，就是第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第一个点一样。periodic：周期性边界条件second： 边界条件，指定定边界二阶导数，在参数valconds给出。二阶导数的缺省值为0, 0variational：自然样条，边界二阶导数为0，如果选择该方法，则忽略

19、给定的二阶导数对于一些特殊的边界条件，可以通过conds的一个1x2矩阵来表示，conds元素的取值可为0，1，2。conds(i)=j的含义是给定端点i的j 阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。详细情况请使用帮助help csape。例： 考虑数据x1245y1342边界条件S(1)=2.5,S(5)=-3,程序如下：x=1 2 4 5;y=1 3 4 2;pp=csape(x,y,second,2.5,-3);pp.coefsxi=1:0.1:5;yi=ppv

20、al(pp,xi);plot(x,y,o,xi,yi);二维网格数据插值已知个节点其中互不相同，不妨设：构造一个二元函数，通过全部已知节点，即 再用计算在处的插值，（1）最临近点插值二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。（2）分片线性插值将四个插值点（矩形的四个顶点）处的函数值依次简记为：f (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4第一片(下三角形区域)：插值函数为：第二片(上三角形区域)插值函数为：注意：(x, y

21、)当然应该是在插值节点所形成的矩形区域内；分片线性插值函数是连续的；（3）双线性插值双线性插值是一片一片的空间二次曲面构成，双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数，但分片双线性插值曲面在所给区域上不连续。（4）分片双三次样条插值分片双三次样条插值函数在每一片（即小矩形）上的表达式为：其中待定系数可由四个顶点的取值以及的光滑性（即的连续性）以及边界条件唯一确定Matlab高维插值命令命令2 interp2功能： 二维数据内插值格式： （1） ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI

22、，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)Xi(i,j),yi(i,j)。用户可以输入行向量和列向量Xi与Yi，此时，输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（Not a Number）。格式（2）：ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第一种情形进行计算。格式（3）：Z

23、I = interp2(Z,n) %作n次递归计算，在Z的%每两个元素之间插入它们%的二维插值，这样，Z的%阶数将不断增加。%interp2(Z)等价于%interp2(z,1)。格式（4）：ZI=interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值：linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。例X,Y = meshgrid(-3:.25:3);Z = peaks(X,Y); %山峰函数XI,YI = meshgrid(-3:.125:3);ZZ = interp2(X,Y,

24、Z,XI,YI);surfl(X,Y,Z);hold on;surfl(XI,YI,ZZ+15)axis(-3 3 -3 3 -5 20);hold off例：years = 1950:10:1990; %从1950年到1990年，每间隔为十年service = 10:10:30; %服务年限 wage = 150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243; %工资水平w = interp2(service

25、,years,wage,15,1975)插值结果为：w =190.6288例：测得平板表面3*5网格点处的温度分别为： 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.再输入以下命令:xi=1:0.2:5;yi=1:0.2

26、:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi) %画出插值后的温度分布曲面图. 命令3 interp3功能 三维数据插值（查表）格式（1） VI = interp3(X,Y,Z,V,XI,YI,ZI) %找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度，不同方向（行或列）的向量，这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)

27、中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。格式（2）：VI = interp3(V,XI,YI,ZI) %缺省地，X=1:N，Y=1:M，Z=1:P，其中，M,N,P=size(V)，再按上面的情形计算。格式（3）：VI = interp3(V,n) %作n次递归计算，在V的每两个元素之间插入它们的三维插值。这样，V的阶数将不断增加。interp3(V)等价于interp3(V,1)。格式（4）：VI = interp3(,method) %用指定的算法method作插值计算： linear：线性插值（缺省算法）； cubic：三次插值； spline：三次样条插值； n

28、earest：最邻近插值。说明 在所有的算法中，都要求X,Y,Z是单调且有相同的格点形式。当X，Y,Z是等距且单调时，用算法linear，cubic，nearest,可得到快速插值 例 x,y,z,v = flow(20); xx,yy,zz = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);vv = interp3(x,y,z,v,xx,yy,zz); 命令4 interpn功能 n维数据插值格式（1） VI = interpn(X1,X2,Xn,V,Y1,Y2,Yn) %返回由参量X1,X2,Xn,V确定的n元函数V=V(X1,X2,Xn)在点（Y1,Y2

29、,Yn）处的插值。参量Y1,Y2,Yn是同型的矩阵或向量。若Y1,Y2,Yn是向量，则可以是不同长度，不同方向（行或列）的向量。它们将通过命令ndgrid生成同型的矩阵，再作计算。若点(Y1,Y2,Yn)中有位于点（X1,X2,Xn）之外的点，则相应地返回特殊变量NaN。格式(2): VI = interpn(V,Y1,Y2,Yn) %缺省地，X1=1:size(V,1)，X2=1:size(V,2)，Xn=1:size(V,n)，再按上面的情形计算。格式（3）: VI = interpn(V,ntimes) %作ntimes次递归计算，在V的每两个元素之间插入它们的n维插值。这样，V的阶数将

30、不断增加。interpn(V)等价于interpn(V,1)。格式（4）:VI = interpn(,method) %用指定的算法method计算： linear：线性插值（缺省算法）； cubic：三次插值； spline：三次样条插值法； nearest：最邻近插值算法。二维散乱数据插值散乱节点插值问题可简述为：已知求矩形区域R 上散乱分布着个节点以及相对应的竖坐标取值，求二元函数通过这些，并使用该函数计算其他的点的取值。常用的方法是反距离加权平均法，或称Shepard方法。基本思想是：在点处，定义其插值函数的函数值为节点处函数值按与节点距离的某种形式反比作为权重的加权平均。例如，记，则插值函数（曲面）可定义为其中这样定义的插值曲面是全局相关的，对曲面上任一