第X章 MATLAB在拟合与插值中的应用Word下载.docx
《第X章 MATLAB在拟合与插值中的应用Word下载.docx》由会员分享,可在线阅读,更多相关《第X章 MATLAB在拟合与插值中的应用Word下载.docx(13页珍藏版)》请在冰豆网上搜索。
plot(x,y,'
o'
x,y,xi,z,'
:
'
)
画出了原始数据x和y,用'
o'
标出该数据点,在数据点之间,再用直线重画原始数据,并用点'
线,画出多项式数据xi和z。
xlabel('
x'
),ylabel('
y=f(x)'
),title('
SecondOrderCurveFitting'
)
将图作标志。
这些步骤的结果表示于前面的图1中。
多项式阶次的选择是有点任意的。
两点决定一直线或一阶多项式。
三点决定一个平方或2阶多项式。
按此进行,n+1数据点唯一地确定n阶多项式。
于是,在上面的情况下,有11个数据点,我们可选一个高达10阶的多项式。
然而,高阶多项式给出很差的数值特性,我们不应选择比所需的阶次高的多项式。
此外,随着多项式阶次的提高,近似变得不够光滑,因为较高阶次多项式在变零前,可多次求导。
不妨选一个10阶多项式
pp=polyfit(x,y,10);
formatshorte
pp.'
则ans=
-4.6436e+005
2.2965e+006
-4.8773e+006
5.8233e+006
-4.2948e+006
2.0211e+006
-6.0322e+005
1.0896e+005
-1.0626e+004
4.3599e+002
-4.4700e-001
要注意在现在情况下,多项式系数的规模与前面的2阶拟合的比较。
还要注意在最小(-4.4700e-001)和最大(5.8233e+006)系数之间有7个数量级的幅度差。
将这个解作图,并把此图与原始数据及2阶曲线拟合相比较。
zz=polyval(pp,xi);
xi,z,'
xi,zz)
2ndand10thOrdercurveFitting'
在下面的图11.2中,原始数据标以'
,2阶曲线拟合是虚线,10阶拟合是实线。
注意,在10阶拟合中,在左边和右边的极值处,数据点之间出现大的纹波。
当企图进行高阶曲线拟合时,这种纹波现象经常发生。
根据图2,显然,‘越多就越好’在这里不适用。
图22阶和10阶曲线拟合
一维插值
正如曲线拟合所描述的那样,插值定义为对数据点之间函数的估值方法,这些数据点是由某些集合给定。
当我们不能很快地求出所需中间点的函数值时,插值是一个有价值的工具。
例如,当数据点是某些实验测量的结果或是过长的计算过程时,就有这种情况。
举例一维插值,考虑下列问题,(由于手头没有我们专业相关的实验数据,故采用老师经常提到的测温的例子,只是数据名称不一样)12小时内,一小时测量一次室外温度。
数据存储在两个MATLAB变量中。
hours=1:
12;
%indexforhourdatawasrecorded
temps=[589152529313022252724];
%recordedtemperatures
plot(hours,temps,hours,temps,'
+'
)%viewtemperatures
title('
Temperature'
Hour'
DegreesCelsius'
图3在线性插值下室外温度曲线
正如图3看到的,MATLAB画出了数据点线性插值的直线。
为了计算在任意给定时间的温度,人们可试着对可视的图作解释。
另外一种方法,可用函数interp1。
t=interp1(hours,temps,9.3)%estimatetemperatureathour=9.3
t=
22.9000
t=interp1(hours,temps,4.7)%estimatetemperatureathour=4.7
22
t=interp1(hours,temps,[3.26.57.111.7])%findtempatmanypoints!
10.2000
30.0000
30.9000
24.9000
若不采用直线连接数据点,我们可采用某些更光滑的曲线来拟合数据点。
最常用的方法是用一个3阶多项式,即3次多项式,来对相继数据点之间的各段建模,每个3次多项式的头两个导数与该数据点相一致。
这种类型的插值被称为3次样条或简称为样条。
函数interp1也能执行3次样条插值。
(这在工程中经常用到)
t=interp1(hours,temps,9.3,'
spline'
)%estimatetemperatureathour=9.3
21.8577
t=interp1(hours,temps,4.7,'
)%estimatetemperatureathour=4.7
22.3143
t=interp1(hours,temps,[3.26.57.111.7],'
)
9.6734
30.0427
31.1755
25.3820
样条插值得到的结果,与上面所示的线性插值的结果不同。
因为插值是一个估计或猜测的过程,其意义在于,应用不同的估计规则导致不同的结果。
一个最常用的样条插值是对数据平滑。
也就是,给定一组数据,使用样条插值在更细的间隔求值。
例如,
h=1:
0.1:
%estimatetemperatureevery1/10hour
t=interp1(hours,temps,h,'
);
plot(hours,temps,'
-'
hours,temps,'
h,t)%plotcomparativeresults
SpringfieldTemperature'
在图4中,虚线是线性插值,实线是平滑的样条插值,标有'
的是原始数据。
如要求在时间轴上有更细的分辨率,并使用样条插值,我们有一个更平滑、但不一定更精确地对温度的估计。
尤其应注意,在数据点,样条解的斜率不突然改变。
作为这个平滑插值的回报,3次样条插值要求更大量的计算,因为必须找到3次多项式以描述给定数据之间的特征。
图4在不同插值下室外温度曲线
二维插值
二维插值是基于与一维插值同样的基本思想。
然而,正如名字所隐含的,二维插值是对两变量的函数z=f(x,y)进行插值(比如钢筋混凝土实验中的正应力和剪应力都对挠度产生影响)。
这里依然考虑温度问题。
(数据由课件中改动而成)设人们对平板上的温度分布估计感兴趣,给定的温度值取自平板表面均匀分布的格栅。
采集了下列的数据:
width=1:
5;
%indexforwidthofplate(i.e.,thex-dimension)
depth=1:
3;
%indexfordepthofplate(i,e,,they-dimension)
temps=[8281808284;
7963616581;
8484828586]%temperaturedata
temps=
8281808284
7963616581
8484828586
如同在标引点上测量一样,矩阵temps表示整个平板的温度分布。
temps的列与下标depth或y-维相联系,行与下标width或x-维相联系(见图5)。
为了估计在中间点的温度,我们必须对它们进行辨识。
wi=1:
0.2:
%estimateacrosswidthofplate
d=2;
%atadepthof2
zlinear=interp2(width,depth,temps,wi,d);
%linearinterpolation
zcubic=interp2(width,depth,temps,wi,d,'
cubic'
%cubicinterpolation
plot(wi,zlinear,'
wi,zcubic)%plotresults
WidthofPlate'
title(['
TemperatureatDepth='
num2str(d)])
另一种方法,我们可以在两个方向插值。
先在三维坐标画出原始数据,看一下该数据的粗糙程度(见图6)。
mesh(width,depth,temps)%usemeshplot
DepthofPlate'
zlabel('
),axis('
ij'
),grid
图5在深度d=2处的平板温度
图6平板温度
然后在两个方向上插值,以平滑数据。
di=1:
%choosehigherresolutionfordepth
%choosehigherresolutionforwidth
zcubic=interp2(width,depth,temps,wi,di,'
%cubic
mesh(wi,di,zcubic)
该例子清楚地证明了,二维插值更为复杂,只是因为有更多的量要保持跟踪。
interp2的基本形式是interp2(x,y,z,xi,yi,method)。
这里x和y是两个独立变量,z是一个应变量矩阵。
x和y对z的关系是
z(i,:
)=f(x,y(i))和z(:
j)=f(x(j),y).
也就是,当x变化时,z的第i行与y的第i个元素y(i)相关,当y变化时,z的第j列与x的第j个元素x(j)相关,。
xi是沿x-轴插值的一个数值数组;
yi是沿y-轴插值的一个数值数组。
图7二维插值后的平板温度
虽然对于许多应用,函数interp1和interp2是很有用的,但它们限制为对单调向量进行插值。
在某些情况,这个限制太严格。
例如,考虑下面的插值:
x=linspace(0,5);
y=1-exp(-x).*sin(2*pi*x);
plot(x,y)
图8函数1-exp(-x).*sin(2*pi*x)的曲线
函数interp1可用来在任何值或x的值上估计y值。
yi=interp1(x,y,1.8)
yi=
1.1556
然而,interp1不能找出对应于某些y值的x值。
例如,如在图8上所示,考虑寻找y=1.1处的x值:
图8给y值在函数曲线上求x的值
plot(x,y,[0,5],[1.11.1])
从图8上,我们看到有四个交点。
使用interp1,我们得到:
xi=interp1(y,x,1.1)
?
?
Errorusing==>
table1
Firstcolumnofthetablemustbemonotonic.
这个函数interp1失败,由于y不是单调的。
如何消除了单调性的要求(我尝试搜索了一些资料,对此问题可由如下解答)
table=[x;
y].'
;
%createcolumnorientedtablefromdata
xi=mminterp(table,2,1.1)
xi=
0.52811.1000
0.95801.1000
1.58251.1000
1.88471.1000
这里使用了线性插值,函数mminterp估计了y=1.1处的四个点。
由于函数mminterp的一般性质,要插值的数据是由面向列矩阵给出,在上面的例子中称作为表(table)。
第二个输入参量是被搜索矩阵table的列,第三个参量是要找的值。
函数的主体由下面给出:
functiony=mminterp(tab,col,val)
%MMINTERP1-DTableSearchbyLinearInterpolation.
%Y=MMINTERP(TAB,COL,VAL)linearlyinterpolatesthetable
%TABsearchingforthescalarvalueVALinthecolumnCOL.
%AllcrossingsarefoundandTAB(:
COL)neednotbemonotonic.
%EachcrossingisreturnedasaseparaterowinYandYhasas
%manycolumnsasTAB.Naturally,thecolumnCOLofYcontains
%thevalueVAL.IfVALisnotfoundinthetable,Y=[].
[rt,ct]=size(tab);
iflength(val)>
1,error('
VALmustbeascalar.'
),end
ifcol>
ct|col<
Chosencolumnoutsidetablewidth.'
ifrt<
2,error('
Tabletoosmallornotorientedincolumns.'
above=tab(:
col)>
val;
%Truewhere>
VAL
below=tab(:
col)<
%Truewhere<
equal=tab(:
col)==val;
%Truewhere=VAL
ifall(above==0)|all(below==0),%handlesimplestcase
y=tab(find(equal),:
);
return
end
pslope=find(below(1:
rt-1)&
above(2:
rt));
%indiceswhereslopeis+
nslope=find(below(2:
rt)&
above(1:
rt-1));
%indiceswhereslopeis-
ib=sort([pslope;
nslope+1]);
%putindicesbelowinorder
ia=sort([nslope;
pslope+1]);
%putindicesaboveinorder
ie=find(equal);
%indiceswhereequaltoval
[tmp,ix]=sort([ib,ie]);
%findwhereequalsfitinresult
ieq=ix>
length(ib);
%Truewhereequalsvaluesfit
ry=length(tmp);
%#ofrowsinresulty
y=zeros(ry,ct);
%pokedataintoazeromatrix
alpha=(val-tab(ib,col))./(tab(ia,col)-tab(ib,col));
alpha=alpha(:
ones(1,ct));
%duplicateforallcolumns
y(~ieq,:
)=alpha.*tab(ia,:
)+(1-alpha).*tab(ib,:
%interpolatedvalues
y(ieq,:
)=tab(ie,:
%equalvalues
y(:
col)=val*ones(ry,1);
%removeroundofferror
正如所见的,mminterp利用了find和sort函数、逻辑数组和数组操作技术。
没有For循环和While循环。
不论用其中哪一种技术来实现将使运行变慢,尤其对大的表。
mminterp与含有大于或等于2的任意数列的表一起工作,如同函数interp1一样。
而且,在这种情况下,插值变量可以是任意的列。
z=sin(pi*x);
%addmoredatatotable
y;
z].'
t=mminterp(table,2,1.1)%sameinterpolationasearlier
t=
0.52811.10000.9930
0.95801.10000.1314
1.58251.1000-0.9639
1.88471.1000-0.3533
t=mminterp(table,3,-.5)%secondthirdcolumnnow
1.16690.7316-0.5000
1.83291.1377-0.5000
3.16710.9639-0.5000
3.83311.0187-0.5000
这些最后的结果估计了x和y在z=-0.5处的值。
小结
曲线的插值和拟合是一个很复杂的工作,但在MATLAB中能由几句轻松的命令来实现,为工程技术人员和科研工作者带来极大的方便,让人不禁感叹它的强大,实为工科学生必备之工具!
!
下面的表总结了在MATLAB中所具有的曲线拟合和插值函数。
可供同学们参考。
。
曲线拟合和插值函数
polyfit(x,y,n)
对描述n阶多项式y=f(x)的数据
进行最小二乘曲线拟合
interp1(x,y,xo)
1维线性插值
interp1(x,y,xo,'
1维3次样条插值
1维3次插值
interp2(x,y,Z,xi,yi)
2维线性插值
interp2(x,y,Z,xi,yi,'
2维3次插值
nearest'
2维最近邻插值
参考文献《精通Matlab综合辅导与指南》