第九講 數(shù)據(jù)插值與擬合

1、第九講第九講 曲線(xiàn)擬合與插值曲線(xiàn)擬合與插值在工程實(shí)踐和科學(xué)實(shí)驗(yàn)中，常常需要從一組實(shí)驗(yàn)觀測(cè)數(shù)據(jù) niyxii, 1 , 0),(揭示自變量x與因變量y之間的關(guān)系， 一般可以用一個(gè)近似的函數(shù)關(guān)系式y(tǒng)f（x）來(lái)表示 通常可以采用兩種方法：曲線(xiàn)擬合和插值 擬合主要是考慮到觀測(cè)數(shù)據(jù)受隨機(jī)誤差的影響，尋求整體誤差最小、較好反映觀測(cè)數(shù)據(jù)的近似函數(shù)，并不保證所得到的函數(shù)一定滿(mǎn)足 )(iixfy 曲線(xiàn)擬合的目的是根據(jù)實(shí)驗(yàn)獲得的數(shù)據(jù)去建立因變量與自變量之間有效的經(jīng)驗(yàn)函數(shù)關(guān)系，為進(jìn)一步的深入研究提供線(xiàn)索 插值函數(shù)一般是已知函數(shù)的線(xiàn)性組合或者稱(chēng)為加權(quán)平均插值在工程實(shí)踐和科學(xué)實(shí)驗(yàn)中有著非常廣泛而又十分重要的應(yīng)用，例如，

2、信息技術(shù)中的圖像重建、圖像放大中為避免圖像的扭曲失真的插值補(bǔ)點(diǎn)、建筑工程的外觀設(shè)計(jì)?；瘜W(xué)工程實(shí)驗(yàn)數(shù)據(jù)與模型的分析、天文觀測(cè)數(shù)據(jù)、地理信息數(shù)據(jù)的處理如(天氣預(yù)報(bào)）以及社會(huì)經(jīng)濟(jì)現(xiàn)象的統(tǒng)計(jì)分析等等插值則要求函數(shù)在每個(gè)觀測(cè)點(diǎn)處一定要滿(mǎn)足 )(iixfy 1、船在該海域會(huì)擱淺嗎船在該海域會(huì)擱淺嗎 在某海域測(cè)得一些點(diǎn)(x，y)處的水深z（單位：英尺）由下表給出，水深數(shù)據(jù)是在低潮時(shí)測(cè)得的船的吃水深度為5英尺，問(wèn)在矩形區(qū)域(75,200)*(-50，150）里的哪些地方船要避免進(jìn)人 一、實(shí)例及其模型分析分析由于測(cè)量點(diǎn)是散亂分布的，先在平面上作出測(cè)量點(diǎn)的分布圖，再利用二維插值方法補(bǔ)充一些點(diǎn)的水深，然后作出海底曲

3、面圖和等高線(xiàn)圖，并求出水深小于5的海域范圍在化學(xué)反應(yīng)中，為研究某化合物的濃度隨時(shí)間的變化規(guī)律，測(cè)得一組數(shù)據(jù)如表2、濃度的變化規(guī)律、濃度的變化規(guī)律表中的數(shù)據(jù)反映了濃度隨時(shí)間變化的函數(shù)關(guān)系，它是一種離散關(guān)系若需要推斷20，40分鐘時(shí)的濃度值，能否用一個(gè)顯函數(shù)y=f(t)來(lái)擬合表中的離散數(shù)據(jù)，然后再計(jì)算濃度值f(20), f（40）?問(wèn)題分析問(wèn)題分析（1）首先將這些離散數(shù)據(jù)分布在直角坐標(biāo)系下，由此可發(fā)現(xiàn)濃度與時(shí)間之間呈現(xiàn)什么規(guī)律這種數(shù)據(jù)分布在直角坐標(biāo)系下的圖形被稱(chēng)為散點(diǎn)圖；（2）根據(jù)散點(diǎn)圖，判段它接近于哪類(lèi)函數(shù)曲線(xiàn)， 即確定函數(shù)形式（3）函數(shù)形式確定以后，關(guān)鍵是要確定函數(shù)中含有的 待定參數(shù)。 最常用

4、的確定待定系數(shù)的方法是，曲線(xiàn)擬合的最小二乘法 二、二、 插值與擬合插值與擬合1、插值方法、插值方法（1）分段線(xiàn)性插值）分段線(xiàn)性插值分段線(xiàn)性插值的提法如下： （2）分段三次埃爾米特插值分段三次埃爾米特插值 在插值問(wèn)題中，如果除了插值節(jié)點(diǎn)的函數(shù)值給定外，還要求在節(jié)點(diǎn)的導(dǎo)數(shù)值為給定值，即插值問(wèn)題變?yōu)橄喈?dāng)于在每一小段上應(yīng)滿(mǎn)足四個(gè)條件（方程），可以確定四個(gè)待定參數(shù)三次多項(xiàng)式正好有四個(gè)系數(shù)，所以可以考慮用三次多項(xiàng)式函數(shù)作為插值函數(shù)，這就是分段三分段三次埃爾米特插值次埃爾米特插值，它與分段線(xiàn)性插值一起都稱(chēng)為分段多項(xiàng)式插值 （3）三次樣條插值）三次樣條插值 2、曲線(xiàn)擬合的最小二乘法、曲線(xiàn)擬合的最小二乘法給定平

5、面上的點(diǎn), 2 , 1),(niyxi進(jìn)行曲線(xiàn)擬合有多種方法，最小二乘法是解決曲線(xiàn)擬合最常用的一種方法 最小二乘法的原理是求f(x),使 niiiniiyxf1212)(達(dá)到最小簡(jiǎn)單地說(shuō)，最小二乘法準(zhǔn)則就是使所有散點(diǎn)到曲線(xiàn)的距離平方和最小 線(xiàn)性最小二乘法線(xiàn)性最小二乘法擬合函數(shù)可由一些簡(jiǎn)單的“基函數(shù)”（例如冪函數(shù)，三角函數(shù)等等） )(,),(),(10 xxxn來(lái)線(xiàn)性表示)()()()(1100 xcxcxcxfmm現(xiàn)在要確定系數(shù) ,10mccc使達(dá)到極小為此 三、插值的matlab實(shí)現(xiàn)MATLAB中的插值函數(shù)為interp1，其調(diào)用格式為 ) ,( 1intmethodxiyxerpyi 其中

6、x，y為插值點(diǎn)，yi為在被插值點(diǎn)xi處的插值結(jié)果，x, y為向量。 注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過(guò)x的范圍。 linearsplinecubicnearestmethod) ,( 1intmethodxiyxerpyi MATLAB提供的插值方法有幾種 表示采用的插值方法：分段線(xiàn)性插值pchip：三次Hermite插值（立方插值）：三次分段樣條插值：最近點(diǎn)等值方式缺省時(shí)表示線(xiàn)性插值 例例1 在一 天24小時(shí)內(nèi)，從零點(diǎn)開(kāi)始每間隔2小時(shí)測(cè)得的環(huán)境溫度數(shù)據(jù)分別為 12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，推測(cè)中午（即13點(diǎn)）時(shí)的溫度x=0

7、：2：24；y=12 9 9 10 18 24 28 27 25 20 18 15 13；x113 ；y1interp1（x，y，x1，spline）若要得到一天24小時(shí)的溫度曲線(xiàn) x=0：2：24；y=12 9 9 10 18 24 28 27 25 20 18 15 13xi0：13600：24；yi=interp1(x,y,xi,spline );plot(x, y, o, xi, yi)2、高維插值、高維插值 N維插值函數(shù)interpN（） 其中N可以為2，3，如N2為二維插值，調(diào)用格式為) ,(2intmethodyixizyxerpzi 其中 x，y，z為插值節(jié)點(diǎn)，zi為被插值點(diǎn)(

8、xi，yi)處的插值結(jié)果 且， xi， yi為被插值節(jié)點(diǎn)構(gòu)成的新的網(wǎng)格數(shù)據(jù)methods代表的意思和可選擇的插值方法和前面一樣注意：注意：所有的插值方法都要求x和y是單調(diào)的網(wǎng)格，x和 y可以是等距的也可以是不等距的 (1) 網(wǎng)格數(shù)據(jù)插值問(wèn)題例例2 氣旋變化情況的可視化下表是氣象學(xué)家測(cè)量得到的氣象資料，它們分別表示在南半球地時(shí)按不同緯度。不同月份的平均氣旋數(shù)字根據(jù)這些數(shù)據(jù)，繪制出氣旋分布曲面圖形y=5:10:85;x=1:12;x,y=meshgrid(x,y);plot(x,y,*);pausez=2.4,1.6,2.4,3.2,1.0,0.5,0.4,0.2,0.5,0.8,2.4,3.6;

9、 18.7 21.4 16.2 9.2 2.8 1.7 1.4 2.4 5.8 9.2 10.3 16; 20.8 18.5 18.2 16.6 12.9 10.1 8.3 11.2 12.5 21.1 23.9 25.5; 22.1 20.1 20.5 25.1 29.2 32.6 33.0 31.0 28.6 32.0 28.1 25.6; 37.3 28.8 27.8 37.2 40.3 41.7 46.2 39.9 35.9 40.3 38.2 43.4; 48.2 36.6 35.5 40 37.6 35.4 35 34.7 35.7 39.5 40 41.9; 25.6 24.2

10、25.5 24.6 21.1 22.2 20.2 21.2 22.6 28.5 25.3 24.3; 5.3 5.3 5.4 4.9 4.9 7.1 5.3 7.3 7 8.6 6.3 6.6; 0.3,0,0,0.3,0,0,0.1,0.2,0.3,0,0.1,0.3;figuresurf(x,y,z)pausexi,yi=meshgrid(1:12,5:1:85);zi=interp2(x,y,z,xi,yi,spline );figuremesh(xi,yi,zi)xlabel(月份),ylabel(緯度),zlabel(氣旋),axis(0 12 0 90 0 50)title(南半球

11、氣旋可視化圖形)（2）、一般二維分布的數(shù)據(jù)插值）、一般二維分布的數(shù)據(jù)插值 在實(shí)際應(yīng)用問(wèn)題中，大部分的數(shù)據(jù)以實(shí)測(cè)的多組（xi,yi,zi）給出，所以不能直接使用interp2（）函數(shù)。 Matlab中提供了另一個(gè)函數(shù)griddata( )，用來(lái)專(zhuān)門(mén)解決這類(lèi)問(wèn)題。其調(diào)用格式如下Z=griddata(x,y,z,x0,y0,method)x,y,z是已知樣本點(diǎn)的坐標(biāo)，可以是任意分布的。X0,y0是期望的插值位置，即被插值節(jié)點(diǎn), 可以是單點(diǎn)，向量或者網(wǎng)格型矩陣插值方法，除了上面的 方法外，還有一個(gè)是4.0版本提供的一個(gè)插值方法，選項(xiàng)為v4四、曲線(xiàn)擬合的四、曲線(xiàn)擬合的matlab實(shí)現(xiàn)實(shí)現(xiàn)1、已知函數(shù)原型

12、的、已知函數(shù)原型的（1）多項(xiàng)式擬合）多項(xiàng)式擬合假設(shè)已知函數(shù)原型為 11nnnaxaxayMatlab提供的擬合函數(shù)為 a=polyfit(xdata,ydata,n)其中n表示多項(xiàng)式的最高階數(shù)，xdata，ydata為將要擬合的數(shù)據(jù)，它是用數(shù)組的方式輸入 輸出參數(shù)a為擬合多項(xiàng)式 的系數(shù),11nnaaaa 注：注：多項(xiàng)式在x處的值y可用下面程序計(jì)算 y=polyval(a,x)T=19.1 25.0 30.1 36.0 40.0 45.1 50.0;R=76.30 77.80 79.25 80.80 82.35 83.90 85.10;PR=polyfit(T,R,1);t=10:60;r=po

13、lyval(PR,t);plot(T,R,*,t,r)解：Matlab程序（2）一般函數(shù)線(xiàn)性組合的曲線(xiàn)擬合）一般函數(shù)線(xiàn)性組合的曲線(xiàn)擬合假設(shè)已知函數(shù)原型為 )()()()(1100 xcxcxcxfmm通過(guò)求解線(xiàn)性方程可得待定系數(shù)，一般方法：X= %已知數(shù)據(jù)x的列向量Y= %已知數(shù)據(jù)y的列向量A=f1(X),f2(X),fm(X) %系數(shù)矩陣，fm（）為基函數(shù)c=Ay解：matlab程序X=0 0.1 0.2 0.3 0.4 0.5 0.6Y=2 2.20254 2.40715 2.61592 2.83096 3.05448 3.28876A=ones(size(X),exp(X),exp(-X

14、);c1=AY;C=c1x=0:0.05:1;y=C(1)+C(2)*exp(x)+C(3)*exp(-x);plot(X,Y,*,x,y)（3）一般的曲線(xiàn)擬合）一般的曲線(xiàn)擬合假設(shè)已知函數(shù)原型是一般的函數(shù)，可以是多項(xiàng)式，可以是線(xiàn)性，也可以是非線(xiàn)性的，一般情況下用這個(gè)來(lái)求解非線(xiàn)性情況 Matlab在優(yōu)化工具箱中提供的求解一般的曲線(xiàn)擬合函數(shù)lsqcurvefit(),其調(diào)用格式如下其中Fun表示函數(shù)Fun(p,data)的M函數(shù)文件，p0表示函數(shù)的初值.。p=lsqcurvefit(Fun,p0,xdata,ydata)注：若要求解點(diǎn)x處的函數(shù)值可用程序f=Fun(p,x)計(jì)算 (1)函數(shù)原型m文

15、件function y=fname(a,t)y=a(1)*exp(-a(2)*t);注：注：因?yàn)楹竺婵赡芤玫接?jì)算函數(shù)在一些點(diǎn)上的值，因此寫(xiě)函數(shù)原型表達(dá)式時(shí)，記得用點(diǎn)運(yùn)算點(diǎn)運(yùn)算tk=0.2 0.3 0.4 0.5 0.6 0.7 0.8;Ik=3.16 2.38 1.75 1.34 1.00 0.74 0.56;a=lsqcurvefit(fname,1,1,tk,Ik)x=0:0.05:1;y=fname(a,x);plot(tk,Ik,*,x,y)運(yùn)行結(jié)果：a = 5.6361 2.8906(2)擬合程序2、函數(shù)原型未知、函數(shù)原型未知已知一組數(shù)據(jù)，用什么樣的曲線(xiàn)擬合最好呢？可以根據(jù)散點(diǎn)圖進(jìn)行直觀判斷，在此基礎(chǔ)上，選擇幾種曲線(xiàn)分別擬合，然后觀察哪條曲線(xiàn)的最小二乘指標(biāo)最小。圖(a)，數(shù)據(jù)接近于直線(xiàn)，故宜采用線(xiàn)性函數(shù)y=a+bx擬合； 解解：先將表中的數(shù)據(jù)用曲線(xiàn)表示X=1.1052 1.2214 1.3499 1.4918 1.6487 1.8221 2.0138 2.2255 2.4596 2.7183 3.6693;Y=0.6795 0.6006 0.5309 0.4693 0.4148 0.3666 0.3241 0.2865 0.2532 0.22