曲线拟合和插值.docx
《曲线拟合和插值.docx》由会员分享,可在线阅读,更多相关《曲线拟合和插值.docx(30页珍藏版)》请在冰豆网上搜索。
曲线拟合和插值
曲线拟合与插值
在大量的应用领域中,人们常常面临用一个解析函数描述数据(一般是测量值)的任务。
对那个问题有两种方式。
在插值法里,数据假定是正确的,要求以某种方式描述数据点之间所发生的情形。
这里讨论的方式是曲线拟合或回归。
人们设法找出某条滑腻曲线,它最佳地拟合数据,但没必要要通过任何数据点。
图说明了这两种方式。
标有'o'的是数据点;连接数据点的实线描画了线性内插,虚线是数据的最佳拟合。
1曲线拟合
曲线拟合涉及回答两个大体问题:
最佳拟合意味着什么?
应该用什么样的曲线?
可用许多不同的方式概念最佳拟合,并存在无穷数量的曲线。
所以,从这里开始,咱们走向何方?
正如它证明的那样,当最佳拟合被解释为在数据点的最小误差平方和,且所用的曲线限定为多项式时,那么曲线拟合是相当简捷的。
数学上,称为多项式的最小二乘曲线拟合。
如图。
虚线和标志的数据点之间的垂直距离是在该点的误差。
对各数据点距离求平方,并把平方距离全加起来,就是误差平方和。
这条虚线是使误差平方和尽可能小的曲线,即是最佳拟合。
最小二乘那个术语仅仅是使误差平方和最小的省略说法。
图2阶曲线拟合
在MATLAB中,函数polyfit求解最小二乘曲线拟合问题。
为了论述那个函数的用法,让咱们以上面图中的数据开始。
»x=[0.1.2.3.4.5.6.7.8.91];
»y=[];
为了用polyfit,咱们必需给函数给予上面的数据和咱们希望最佳拟合数据的多项式的阶次或度。
若是咱们选择n=1作为阶次,取得最简单的线性近似。
通常称为线性回归。
相反,若是咱们选择n=2作为阶次,取得一个2阶多项式。
此刻,咱们选择一个2阶多项式。
»n=2;%polynomialorder
»p=polyfit(x,y,n)
p=
polyfit的输出是一个多项式系数的行向量。
其解是y=-+-。
为了将曲线拟合解与数据点比较,让咱们把二者都绘成图。
»xi=linspace(0,1,100);%x-axisdataforplotting
»z=polyval(p,xi);
为了计算在xi数据点的多项式值,挪用MATLAB的函数polyval。
»plot(x,y,'o',x,y,xi,z,':
')
画出了原始数据x和y,用'o'标出该数据点,在数据点之间,再用直线重画原始数据,并用点':
'线,画出多项式数据xi和z。
»xlabel('x'),ylabel('y=f(x)'),title('SecondOrderCurveFitting')
将图作标志。
这些步骤的结果表示于前面的图中。
多项式阶次的选择是有点任意的。
两点决定一直线或一阶多项式。
三点决定一个平方或2阶多项式。
按此进行,n+1数据点唯一地肯定n阶多项式。
于是,在上面的情形下,有11个数据点,咱们可选一个高达10阶的多项式。
但是,高阶多项式给出很差的数值特性,人们不该选择比所需的阶次高的多项式。
另外,随着多项式阶次的提高,近似变得不够滑腻,因为较高阶次多项式在变零前,可多次求导。
例如,选一个10阶多项式
»pp=polyfit(x,y,10);
»formatshorte%changedisplayformat
»pp.'%displaypolynomialcoefficientsasacolumn
ans=
+005
+006
+006
+006
+006
+006
+005
+005
+004
+002
要注意在此刻情形下,多项式系数的规模与前面的2阶拟合的比较。
还要注意在最小和最大+006)系数之间有7个数量级的幅度差。
将那个解作图,并把此图与原始数据及2阶曲线拟合相较较,结果如何呢?
»zz=polyval(pp,xi);%evaluate10thorderpolynomial
»plot(x,y,'o',xi,z,':
',xi,zz)%plotdata
»xlabel('x'),ylabel('y=f(x)'),title('2ndand10thOrdercurveFitting')
在下面的图中,原始数据标以'o',2阶曲线拟合是虚线,10阶拟合是实线。
注意,在10阶拟合中,在左侧和右边的极值处,数据点之间出现大的纹波。
当企图进行高阶曲线拟合时,这种纹波现象常常发生。
按照图,显然,‘越多就越好’的观念在这里不适用。
图2阶和10阶曲线拟合
一维插值
正如在前一节对曲线拟合所描述的那样,插值概念为对数据点之间函数的估值方式,这些数据点是由某些集合给定。
当人们不能专门快地求出所需中间点的函数值时,插值是一个有价值的工具。
例如,当数据点是某些实验测量的结果或是太长的计算进程时,就有这种情形。
或许最简单插值的例子是MATLAB的作图。
按缺省,MATLAB用直线连接所用的数据点以作图。
那个线性插值猜想中间值落在数据点之间的直线上。
固然,当数据点个数的增加和它们之间距离的减小时,线性插值就更精准。
例如,
»x1=linspace(0,2*pi,60);
»x2=linspace(0,2*pi,6);
»plot(x1,sin(x1),x2,sin(x2),'-')
»xlabel('x'),ylabel('sin(x)'),title('LinearInterpolation')
图线性插值
图是sine函数的两个图,一个在数据点之间用60个点,它比另一个只用6个点更滑腻和更精准。
如曲线拟合一样,插值要作决策。
按照所作的假设,有多种插值。
而且,能够在一维以上空间中进行插值。
即若是有反映两个变量函数的插值,z=f(x,y),那么就可在x之间和在y之间,找出z的中间值进行插值。
MATLAB在一维函数interp1和在二维函数interp2中,提供了许多的插值选择。
其中的每一个函数将在下面论述。
为了说明一维插值,考虑下列问题,12小时内,一小时测量一次室外温度。
数据存储在两个MATLAB变量中。
»hours=1:
12;%indexforhourdatawasrecorded
»temps=[589152529313022252724];%recordedtemperatures
»plot(hours,temps,hours,temps,'+')%viewtemperatures
»title('Temperature')
»xlabel('Hour'),ylabel('DegreesCelsius')
图在线性插值下室外温度曲线
正如图看到的,MATLAB画出了数据点线性插值的直线。
为了计算在任意给按时刻的温度,人们可试着对可视的图作解释。
另外一种方式,可用函数interp1。
»t=interp1(hours,temps,%estimatetemperatureathour=
t=
»t=interp1(hours,temps,%estimatetemperatureathour=
t=
22
»t=interp1(hours,temps,[])%findtempatmanypoints!
t=
interp1的缺省用法是由interp1(x,y,xo)来描述,这里x是独立变量(横坐标),y是应变量(纵坐标),xo是进行插值的一个数值数组。
另外,该缺省的利用假定为线性插值。
若不采用直线连接数据点,咱们可采用某些更滑腻的曲线来拟合数据点。
最常常利用的方式是用一个3阶多项式,即3次多项式,来对接踵数据点之间的各段建模,每一个3次多项式的头两个导数与该数据点相一致。
这种类型的插值被称为3次样条或简称为样条。
函数interp1也能执行3次样条插值。
»t=interp1(hours,temps,,'spline')%estimatetemperatureathour=
t=
»t=interp1(hours,temps,,'spline')%estimatetemperatureathour=
t=
»t=interp1(hours,temps,[],'spline')
t=
注意,样条插值取得的结果,与上面所示的线性插值的结果不同。
因为插值是一个估量或猜想的进程,其意义在于,应用不同的估量规则致使不同的结果。
一个最常常利用的样条插值是对数据光滑。
也就是,给定一组数据,利用样条插值在更细的距离求值。
例如,
»h=1:
:
12;%estimatetemperatureevery1/10hour
»t=interp1(hours,temps,h,'spline');
»plot(hours,temps,'-',hours,temps,'+',h,t)%plotcomparativeresults
»title('SpringfieldTemperature')
»xlabel('Hour'),ylabel('DegreesCelsius')
在图中,虚线是线性插值,实线是光滑的样条插值,标有'+'的是原始数据。
如要求在时刻轴上有更细的分辨率,并利用样条插值,咱们有一个更光滑、但不必然更精准地对温度的估量。
尤其应注意,在数据点,样条解的斜率不突然改变。
作为那个光滑插值的回报,3次样条插值要求更大量的计算,因为必需找到3次多项式以描述给定数据之间的特征。
图在不同插值下室外温度曲线
在讨论二维插值之前,了解interp1所强制的二个强约束是很重要的。
第一,人们不能要求有独立变量范围之外的结果,例如,interp1(hours,temps,致使一个错误,因为hours在1到12之间转变。
第二,独立变量必需是单调的。
即独立变量在值上必需老是增加的或老是减小的。
在咱们的例子里,hours是单调的。
但是,若是咱们已经概念独立变量为一天的实际时刻,
»time_of_day=[7:
121:
6]%startat7AM,endat6PM
time_of_day=
789101112123456
则独立变量将不是单调的,因为time_of_day增加到12,然后跌到1,再然后增加。
若是用time_of_day代替interp1中的hours,将会返回一个错误。
一样的理由,人们不能对temps插值来找出产生某温度的时刻(小时),因为temps不是单调的。
二维插值
二维插值是基于与一维插值一样的大体思想。
但是,正如名字所隐含的,二维插值是对两变量的函数z=f(x,y)进行插值。
为了说明那个附加的维数,考虑一个问题。
设人们对平板上的温度散布估量感兴趣,给定的温度值取自平板表面均匀散布的格栅。
收集了下列的数据:
»width=1:
5;%indexforwidthofplate.,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-维相联系(见图。
为了估量在中间点的温度,咱们必需对它们进行辨识。
»wi=1:
:
5;%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
»xlabel('WidthofPlate'),ylabel('DegreesCelsius')
»title(['TemperatureatDepth='num2str(d)])
另一种方式,咱们能够在两个方向插值。
先在三维坐标画出原始数据,看一下该数据的粗糙程度(见图。
»mesh(width,depth,temps)%usemeshplot
»xlabel('WidthofPlate'),ylabel('DepthofPlate')
»zlabel('DegreesCelsius'),axis('ij'),grid
图在深度d=2处的平板温度
图平板温度
然后在两个方向上插值,以光滑数据。
»di=1:
:
3;%choosehigherresolutionfordepth
»wi=1:
:
5;%choosehigherresolutionforwidth
»zcubic=interp2(width,depth,temps,wi,di,'cubic');%cubic
»mesh(wi,di,zcubic)
»xlabel('WidthofPlate'),ylabel('DepthofPlate')
»zlabel('DegreesCelsius'),axis('ij'),grid
上面的例子清楚地证明了,二维插值更为复杂,只是因为有更多的量要维持跟踪。
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-轴插值的一个数值数组。
图二维插值后的平板温度
可选的参数method能够是'linear','cubic'或'nearest'。
在这种情形下,cubic不意味着3次样条,而是利用3次多项式的另一种算法。
linear方式是线性插值,仅用作连接图上数据点。
nearest方式只选择最接近各估量点的粗略数据点。
在所有的情形下,假定独立变量x和y是线性距离和单调的。
M文件举例
虽然对于许多应用,函数interp1和interp2是很有效的,但它们限制为对单调向量进行插值。
在某些情形,那个限制太严格。
例如,考虑下面的插值:
»x=linspace(0,5);
»y=1-exp(-x).*sin(2*pi*x);
»plot(x,y)
图函数1-exp(-x).*sin(2*pi*x)的曲线
函数interp1可用来在任何值或x的值上估量y值。
»yi=interp1(x,y,
yi=
但是,interp1不能找出对应于某些y值的x值。
例如,如在图上所示,考虑寻觅y=处的x值:
图给y值在函数曲线上求x的值
»plot(x,y,[0,5],[])
从图上,咱们看到有四个交点。
利用interp1,咱们取得:
»xi=interp1(y,x,
Errorusing==>table1
Firstcolumnofthetablemustbemonotonic.
那个函数interp1失败,由于y不是单调的。
在本章精通MATLAB工具箱所说明的M文件例子,消除单调性的要求。
»table=[x;y].';%createcolumnorientedtablefromdata
»xi=mminterp(table,2,
xi=
这里利用了线性插值,函数mminterp估量了y=处的四个点。
由于函数mminterp的一般性质,要插值的数据是由面向列矩阵给出,在上面的例子中称作为表(table)。
第二个输入参量是被搜索矩阵table的列,第三个参量是要找的值。
那个精通MATLAB工具箱函数的主体由下面给出:
functiony=mminterp(tab,col,val)
%MMINTERP1-DTableSearchbyLinearInterpolation.
%Y=MMINTERP(TAB,COL,VAL)linearlyinterpolatesthetable
%TABsearchingforthescalarvalueVALinthecolumnCOL.
%AllcrossingsarefoundandTAB(:
COL)neednotbemonotonic.
%EachcrossingisreturnedasaseparaterowinYandYhasas
%manycolumnsas,thecolumnCOLofYcontains
%thevalueVAL.IfVALisnotfoundinthetable,Y=[].
%Copyright(c)1996byPrentice-Hall,Inc.
[rt,ct]=size(tab);
iflength(val)>1,error('VALmustbeascalar.'),end
ifcol>ct|col<1,error('Chosencolumnoutsidetablewidth.'),end
ifrt<2,error('Tabletoosmallornotorientedincolumns.'),end
above=tab(:
col)>val;%Truewhere>VAL
below=tab(:
col)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
»table=[x;y;z].';
»t=mminterp(table,2,%sameinterpolationasearlier
t=
»t=mminterp(table,3,%secondthirdcolumnnow
t=
这些最后的结果估量了x和y在z=处的值。
虽然逐条地对函数mminterp解释如何工作是很有帮忙的,但如此做要求有更多的篇幅和时刻。
解释mminterp如何工作最容易的方式是创建一个小表格,然后,在重要的语句末尾删除分号以后,挪用函数。
如此,中间值将帮忙用户理解函数是如何找到与所需值相符的数据值和如何执行插值。
前面已论述了interp1的用法。
当用于线性插值时,只要所要求的插值点的个数少,interp1工作专门好。
在要求许多插值点情形下,由于所用的算法,interp1工作较慢。
为了克服那个问题,精通MATLAB工具箱包括了函数mmtable,
MMTABLE由线性插值实现一维单调表搜索
YI=MMTABLE(TAB,COL,VALS)线性地对表TAB进行插值,在列COL中搜索值为VALS
TAB(:
COL)必需是单调的,但没必要等价地生成空间。
YI