实验三插值.docx

实验三插值.docx

《实验三插值.docx》由会员分享，可在线阅读，更多相关《实验三插值.docx（27页珍藏版）》请在冰豆网上搜索。

实验三插值

试验三Matlab插值与拟合

一.实验目的

了解插值的基本内容和原理；掌握用matlab求解插值问题,包括一维插值和二维插值的各种常用方法；

二.实验原理和方法

插值法是实用的数值方法，是函数逼近的重要方法。

在生产和科学实验中，自变量x与因变量y的函数y=f（x）的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。

当要求知道观测点之外的函数值时，需要估计函数值在该点的值。

如何根据观测点的值，构造一个比较简单的函数y=φ（x），使函数在观测点的值等于已知的数值或导数值，进而用简单函数y=φ（x）在点x处的值来估计未知函数y=f（x）在x点的值。

寻找这样的函数φ（x），办法是很多的。

φ（x）可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ（x）可以是任意光滑（任意阶导数连续）的函数或是分段函数；函数类的不同，自然地有不同的逼近效果。

一维插值

已知

个数据节点:

构造一个（相对简单）函数

（称为插值函数）,通过全部结点即

（j＝0,1,…n）再用

计算插值，即

数学上插值方法非常多，只介绍几种常用方法：

（1）拉格朗日（Lagrange）插值

已知函数f（x）（称为被插值函数）在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn。

求一n次多项式函数Pn（x）（称为插值函数），使其满足：

Pn（xi）=yi,i=0,1,…,n..

解决此问题的拉格朗日插值多项式公式如下：

其中

为

次多项式：

称为拉格朗日插值基函数,可以验证该多项式通过所给数据点。

特别地:

两点一次（线性）插值多项式为：

三点二次（抛物）插值多项式为:

matlab没有Lagrange插值的函数，需要自己编写。

Lagrange插值插值的最大优点是函数具有很好的解析性质,无穷次可微。

缺点是：

（1）可能出现严重的振荡现象，在头尾附近会出现一些不好的振荡现象（龙格现象）

（2）多项式的系数对数据非常敏感。

导致Lagrange插值如上缺点的原因在于：

多项式次数过高不一定是好事。

而克服这一缺点的办法之一就是采用分段低次插值，其中最简单的是分段线性插值。

（2）分段线性插值

分段线性插值就是在每个子区间上做通过两个端点的直线。

该方法具有良好的收敛性，即

，其中

为被插值函数，可以通过增加插值节点的个数来控制插值误差，是一种较实用的插值方法，Matlab有分段线性插值函数，但总体光滑程度不够。

（3）分段三次埃尔米特（Hermite）插值

对插值函数的要求可以根据需要进行改变，由此产生了不同的插值方法，分段三次埃尔米特插值就是其一。

在插值问题中，如果除了指定在插值节点的函数值以外，同时指定插值节点处的导数值，这时插值问题变为：

给定一组观察数据

及

，要求一个分段多项式函数

，满足：

这就相当于所求分段多项式在每个小段

上满足四个条件：

可以确定四个待定参数，三次多项式恰好有四个系数，所以可以考虑用三次多项式作为插值函数，这就是分段三次Hermite插值。

（4）分段三次样条插值

在数学上，光滑程度的定量描述是：

函数（曲线）的k阶导数存在且连续，则称该曲线具有k阶光滑性。

光滑性的阶次越高，则曲线光滑程度越好。

分段线性插值具有零阶光滑性，也就是不光滑；分段三次埃尔米特插值具有一阶光滑性，但是这样的光滑程度在工程设计和机械加工等实际应用中是远远不够的。

是否存在较低次的分段多项式达到较高阶光滑性的方法？

三次样条插值就是一个很好的例子。

什么是样条？

所谓样条，就是一个细的、可弯曲的木制或塑料条，在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线（放样）所用的工具。

1946年,Schoenberg将样条引入数学,即所谓的样条函数。

三次样条函数本质上是一段一段的三次多项式拼合而成的曲线，在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的。

满足以上条件的分段函数

称为三次样条插值函数。

可以看出，三次样条插值是分段三次多项式插值，但是不仅要求过所给数据节点，而且要求在所给节点处的二阶导数也是连续的，从而构造出来的插值函数比较光滑。

在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个方程，还差两个。

为此常用的方法是在两个边界节点

和

附加要求，这就是所谓的边界条件。

根据实际问题的不同，三次样条插值常用到下列三类边界条件。

（1）

边界条件：

。

即指定两个边界节点的一阶导数值为给定值

（2）

边界条件：

。

即指定两个边界节点的二阶导数值为给定值

特别地，当

和

都为零时，称为自然边界条件。

（3）周期性边界条件：

总而言之，理论上三次样条插值是确定的。

matlab插值命令

命令1interp1

功能一维数据插值。

该命令对数据点之间计算内插值。

它找出一元函数f（x）在中间点的数值。

其中函数f（x）由所给数据决定。

各个参量之间的关系示意图如下图

图数据点与插值点关系示意图

格式：

（1）yi=interp1（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的默认的边界条件。

即要求第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第二个点一样。

‘pchip’：

分段三次Hermite插值。

对于该方法，命令interp1调用函数pchip，用于对向量x与y执行分段三次内插值。

该方法保留单调性与数据的外形；

‘cubic’：

与’pchip’操作相同；

‘v5cubic’：

在MATLAB5.0中的三次插值。

对于超出x范围的xi的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。

对其他的方法，interp1将对超出的分量执行外插值算法。

（4）yi=interp1（x,Y,xi,method,'extrap'）%对于超出x范围的xi中的分量将执行特殊的外插值法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处的值，可使用yi=ppval（pp,xi）来计算。

ppval（pp,xi）等价于interp1（x,Y,xi,method,'extrap'）.

例：

已知标准正态分布函数值

。

利用分段线性插值计算

x=[2.342.35];

y=[0.990360.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=interp1（x,y,xx）;

plot（x,y,'kd',xx,yy）

例

year=1900:

10:

2010;%年份从1900到2010，每间隔为10年

product=[75.99591.972105.711123.203131.669150.697...

179.323203.212226.505249.633256.344267.893];%产量

p1995=interp1（year,product,1995）%插值计算1995年产量，采用

%默认的线性差值方法。

x=1900:

1:

2010;%1900到2010年，每间隔1年

%采用分段三次Hermite插值

y=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'）,holdon,plot（xi,yi,'r:

'）,holdoff

例：

在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：

5，8，9，15，25，29，31，30，22，25，27，24。

试估计每隔1/10小时的温度值。

解：

hours=1:

12;

temps=[589152529313022252724];

h=1:

0.1:

12;

t=interp1（hours,temps,h,'spline'）;%（直接输出数据将是很多的）

plot（hours,temps,'+',h,t,'r:

'）%作图

xlabel（'Hour'）,ylabel（'DegreesCelsius'）

命令2spline

功能：

三次样条数据插值

格式

（1）：

pp=spline（x,y）

说明：

返回分段多项式形式的三次样条插值函数供函数ppval使用，

Pp为x必须是一个向量.

1.若y也是一个向量，而且length（x）=length（y）则所求的样条函数S（x（i））=y（i）,此时用非扭结（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的要求同格式

（1）。

例：

对离散地分布在y=exp（x）sin（x）函数曲线上的数据点进行样条插值计算：

x=[024581212.817.219.920];

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.151.122.362.361.46.49.060];

cs=spline（x,[0y0]）;%数组[0y0]所含分两个数比x多两个，%使用指定一阶导数的边界条件

xx=linspace（-4,4,101）;

plot（x,y,'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,:

）,'-'）;holdon;

plot（x,Y（2,:

）,'o',xx,YY（2,:

）,':

'）;holdoff;

命令3csape

功能：

使用各种边界条件的三次样条插值，返回PP形式的样条函数。

基本使用格式：

（1）pp=csape（x,y）；%使用默认的边界条件，即not-a-knot条件。

（2）pp=csape（x,y,conds,valconds）

说明：

conds指定边界类型，有如下几种：

'complete':

边界条件,即给定边界一阶导数.一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。

也即

'not-a-knot'：

非扭结条件，就是第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第一个点一样。

'periodic'：

周期性边界条件

'second'：

边界条件，指定定边界二阶导数，在参数valconds给出。

二阶导数的缺省值为[0,0]

'variational'：

自然样条，边界二阶导数为0，如果选择该方法，则忽略给定的二阶导数

对于一些特殊的边界条件，可以通过conds的一个1x2矩阵来表示，conds元素的取值可为0，1，2。

conds（i）=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。

详细情况请使用帮助helpcsape。

例：

考虑数据

x

1

2

4

5

y

1

3

4

2

边界条件S''

（1）=2.5,S''（5）=-3,

程序如下：

x=[1245];y=[1342];

pp=csape（x,y,'second',[2.5,-3]）;

pp.coefs

xi=1:

0.1:

5;

yi=ppval（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）当然应该是在插值节点所形成的矩形区域内；分片线性插值函数是连续的；

（3）双线性插值

双线性插值是一片一片的空间二次曲面构成，双线性插值函数的形式如下：

其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数，但分片双线性插值曲面在所给区域上不连续。

（4）分片双三次样条插值

分片双三次样条插值函数在每一片（即小矩形）上的表达式为：

其中待定系数

可由四个顶点的取值以及

的光滑性（即

的连续性）以及边界条件唯一确定

Matlab高维插值命令

命令2interp2

功能：

二维数据内插值

格式：

（1）ZI=interp2（X,Y,Z,XI,YI）%返回矩阵ZI，其元素包含对应于参量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（NotaNumber）。

格式

（2）：

ZI=interp2（Z,XI,YI）%缺省地，X=1:

n、Y=1:

m，其中[m,n]=size（Z）。

再按第一种情形进行计算。

格式（3）：

ZI=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,Z,XI,YI）;

surfl（X,Y,Z）;holdon;

surfl（XI,YI,ZZ+15）

axis（[-33-33-520]）;

holdoff

例：

years=1950:

10:

1990;%从1950年到1990年，每间隔为十年

service=10:

10:

30;%服务年限

wage=[150.697199.592187.625

179.323195.072250.287

203.212179.092322.767

226.505153.706426.730

249.633120.281598.243];%工资水平

w=interp2（service,years,wage,15,1975）

插值结果为：

w=

190.6288

例：

测得平板表面3*5网格点处的温度分别为：

8281808284

7963616581

8484828586

试作出平板表面的温度分布曲面z=f（x,y）的图形。

1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.

输入以下命令：

x=1:

5;

y=1:

3;

temps=[8281808284;7963616581;8484828586];

mesh（x,y,temps）

2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.

再输入以下命令:

xi=1:

0.2:

5;

yi=1:

0.2:

3;

zi=interp2（x,y,temps,xi,yi’,'cubic'）;

mesh（xi,yi,zi）%画出插值后的温度分布曲面图.

命令3interp3

功能三维数据插值（查表）

格式

（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）中有位于点（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’：

三次样条插值；

‘nearest’：

最邻近插值。

说明在所有的算法中，都要求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）;

命令4interpn

功能n维数据插值

格式

（1）VI=interpn（X1,X2,,…,Xn,V,Y1,Y2,…,Yn）%返回由参量X1,X2,…,Xn,V确定的n元函数V=V（X1,X2,…,Xn）在点（Y1,Y2,…,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的阶数将不断增加。

interpn（V）等价于interpn（V,1）。

格式（4）:

VI=interpn（…,method）%用指定的算法method计算：

‘linear’：

线性插值（缺省算法）；

‘cubic’：

三次插值；

‘spline’：

三次样条插值法；

‘nearest’：

最邻近插值算法。

二维散乱数据插值

散乱节点插值问题可简述为：

已知求矩形区域R上散乱分布着

个

节点

以及相对应的竖坐标取值

，求二元函数

通过这些

，并使用该函数计算其他的点

的取值。

常用的方法是反距离加权平均法，或称Shepard方法。

基本思想是：

在点

处，定义其插值函数的函数值为节点处函数值按

与节点距离的某种形式反比作为权重的加权平均。

例如，记

，则插值函数（曲面）可定义为

其中

这样定义的插值曲面是全局相关的，对曲面上任一