第六講插值與擬合

第六講插值與擬合第一頁，共六十一頁，編輯于2023年，星期一插值與擬合一、插值的基本原理二、擬合的基本原理三、插值與擬合的關(guān)系四、插值的MATLAB實(shí)現(xiàn)五、擬合的Matlab實(shí)現(xiàn)第二頁，共六十一頁，編輯于2023年，星期一

我們經(jīng)常會遇到大量的數(shù)據(jù)需要處理，而處理數(shù)據(jù)的關(guān)鍵就在于這些算法，例如數(shù)據(jù)擬合、參數(shù)估計(jì)、插值等數(shù)據(jù)處理算法。此類問題在MATLAB中有很多現(xiàn)成的函數(shù)可以調(diào)用，熟悉MATLAB，這些方法都能游刃有余的用好。一、概述第三頁，共六十一頁，編輯于2023年，星期一

數(shù)據(jù)擬合在很多賽題中有應(yīng)用，與圖形處理有關(guān)的問題很多與插值和擬合有關(guān)系，例如98年美國賽A題，生物組織切片的三維插值處理，94年A題逢山開路，山體海拔高度的插值計(jì)算，2003年吵的沸沸揚(yáng)揚(yáng)的“非典”問題也要用到數(shù)據(jù)擬合算法，觀察數(shù)據(jù)的走向進(jìn)行處理，2005年的雨量預(yù)報(bào)的評價(jià)的插值計(jì)算。2001年的公交車調(diào)度擬合問題，2003年的飲酒駕車擬合問題。第四頁，共六十一頁，編輯于2023年，星期一喝兩瓶酒的擬合曲線喝1-5瓶酒的擬合曲線第五頁，共六十一頁，編輯于2023年，星期一

在實(shí)際中，常常要處理由實(shí)驗(yàn)或測量所得到的一些離散數(shù)據(jù)。插值與擬合方法就是要通過這些數(shù)據(jù)去確定某一類已知函數(shù)的參數(shù)或?qū)で竽硞€(gè)近似函數(shù)，使所得到的近似函數(shù)與已知數(shù)據(jù)有較高的擬合精度。如果要求這個(gè)近似函數(shù)（曲線或曲面）經(jīng)過所已知的所有數(shù)據(jù)點(diǎn)，則稱此類問題為插值問題。（不需要函數(shù)表達(dá)式）二、基本概念第六頁，共六十一頁，編輯于2023年，星期一

如果不要求近似函數(shù)通過所有數(shù)據(jù)點(diǎn)，而是要求它能較好地反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法稱為數(shù)據(jù)擬合。（必須有函數(shù)表達(dá)式）近似函數(shù)不一定（曲線或曲面）通過所有的數(shù)據(jù)點(diǎn)。第七頁，共六十一頁，編輯于2023年，星期一1、聯(lián)系都是根據(jù)實(shí)際中一組已知數(shù)據(jù)來構(gòu)造一個(gè)能夠反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法。2、區(qū)別插值問題不一定得到近似函數(shù)的表達(dá)形式，僅通過插值方法找到未知點(diǎn)對應(yīng)的值。數(shù)據(jù)擬合要求得到一個(gè)具體的近似函數(shù)的表達(dá)式。

三、插值與擬合的區(qū)別和聯(lián)系第八頁，共六十一頁，編輯于2023年，星期一四、插值的使用及求解

當(dāng)數(shù)據(jù)量不夠，需要補(bǔ)充，且認(rèn)定已有數(shù)據(jù)可信時(shí),通常利用函數(shù)插值方法。實(shí)際問題當(dāng)中碰到的函數(shù)f(x)

是各種各樣的，有的表達(dá)式很復(fù)雜，有的甚至給不出數(shù)學(xué)的式子，只提供了一些離散數(shù)據(jù)，警如，某些點(diǎn)上的函數(shù)值和導(dǎo)數(shù)值。4.1引言第九頁，共六十一頁，編輯于2023年，星期一

選用不同類型的插值函數(shù)，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值）（2）分段線性插值（3）Hermite（4）三次樣條插值。4.2插值方法第十頁，共六十一頁，編輯于2023年，星期一一維插值的定義已知n+1個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)求任一插值點(diǎn)處的插值節(jié)點(diǎn)可視為由產(chǎn)生,，表達(dá)式復(fù)雜,，或無顯式形式,，或未知.。第十一頁，共六十一頁，編輯于2023年，星期一

構(gòu)造一個(gè)(相對簡單的)函數(shù)通過全部節(jié)點(diǎn),即再用計(jì)算插值，即返回第十二頁，共六十一頁，編輯于2023年，星期一

稱為拉格朗日插值基函數(shù)。

已知函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿足：

Pn(xi)=yi,i=0,1,…,n.

解決此問題的拉格朗日插值多項(xiàng)式公式如下其中Li(x)為n次多項(xiàng)式：拉格朗日(Lagrange)插值第十三頁，共六十一頁，編輯于2023年，星期一拉格朗日(Lagrange)插值特別地:兩點(diǎn)一次(線性)插值多項(xiàng)式:三點(diǎn)二次(拋物)插值多項(xiàng)式:第十四頁，共六十一頁，編輯于2023年，星期一

拉格朗日多項(xiàng)式插值的這種振蕩現(xiàn)象叫Runge現(xiàn)象

采用拉格朗日多項(xiàng)式插值：選取不同插值節(jié)點(diǎn)個(gè)數(shù)n+1，其中n為插值多項(xiàng)式的次數(shù)，當(dāng)n分別取2,4,6,8,10時(shí)，繪出插值結(jié)果圖形.例ToMatlablch(larg1)第十五頁，共六十一頁，編輯于2023年，星期一分段線性插值計(jì)算量與n無關(guān);n越大，誤差越小.xjxj-1xj+1x0xnxoy第十六頁，共六十一頁，編輯于2023年，星期一ToMATLABxch11，xch12，xch13，xch14例用分段線性插值法求插值,并觀察插值誤差.1.在[-6,6]中平均選取5個(gè)點(diǎn)作插值(xch11)4.在[-6,6]中平均選取41個(gè)點(diǎn)作插值(xch14)2.在[-6,6]中平均選取11個(gè)點(diǎn)作插值(xch12)3.在[-6,6]中平均選取21個(gè)點(diǎn)作插值(xch13)第十七頁，共六十一頁，編輯于2023年，星期一比分段線性插值更光滑。xyxi-1xiab

在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個(gè)很好的例子。三次樣條插值第十八頁，共六十一頁，編輯于2023年，星期一

三次樣條插值g(x)為被插值函數(shù)。第十九頁，共六十一頁，編輯于2023年，星期一例用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值(ych)ToMATLABych(larg1)第二十頁，共六十一頁，編輯于2023年，星期一Matlab實(shí)現(xiàn)：實(shí)現(xiàn)分段線性插值不需要編制函數(shù)程序，它自身提供了內(nèi)部的功能函數(shù)interp1(一維插值)intep2(二維)interp3(三維)intern(n維)4.3MATLAB實(shí)現(xiàn)插值第二十一頁，共六十一頁，編輯于2023年，星期一用MATLAB作插值計(jì)算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)xi處的插值結(jié)果‘nearest’

最鄰近插值；‘linear’

線性插值；‘spline’

三次樣條插值；‘cubic’

立方插值；缺省時(shí)分段線性插值．

注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍．第二十二頁，共六十一頁，編輯于2023年，星期一例：從1點(diǎn)12點(diǎn)的11小時(shí)內(nèi)，每隔1小時(shí)測量一次溫度，測得的溫度的數(shù)值依次為：5，8，9，15，25，29，31，30，22，25，27，24．試估計(jì)每隔1/10小時(shí)的溫度值．ToMATLAB(temp)hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)第二十三頁，共六十一頁，編輯于2023年，星期一xy機(jī)翼下輪廓線例已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變0.1時(shí)的y值．ToMATLAB(plane)返回第二十四頁，共六十一頁，編輯于2023年，星期一

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍．z=interp2(x0,y0,z0,x,y,’method’)被插值點(diǎn)插值方法用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’

最鄰近插值；‘linear’

雙線性插值；‘cubic’

雙三次插值；缺省時(shí)雙線性插值.第二十五頁，共六十一頁，編輯于2023年，星期一例：測得平板表面3×5網(wǎng)格點(diǎn)處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形．輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲線圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.第二十六頁，共六十一頁，編輯于2023年，星期一再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.ToMATLAB(wendu)第二十七頁，共六十一頁，編輯于2023年，星期一

通過此例對最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較．ToMATLAB(moutain)返回第二十八頁，共六十一頁，編輯于2023年，星期一

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

要求cx取行向量，cy取為列向量．被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值'v4'-MATLAB提供的插值方法缺省時(shí),雙線性插值第二十九頁，共六十一頁，編輯于2023年，星期一

例在某海域測得一些點(diǎn)(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）×（-50，150）里的哪些地方船要避免進(jìn)入．第三十頁，共六十一頁，編輯于2023年，星期一ToMATLABhd1返回4.作出水深小于5的海域范圍,即z=5的等高線.2.在矩形區(qū)域(75,200)×(-50,150)進(jìn)行插值。1.輸入插值基點(diǎn)數(shù)據(jù)3.作海底曲面圖第三十一頁，共六十一頁，編輯于2023年，星期一%程序一：插值并作海底曲面圖

x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];z=[48686889988949];x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');meshc(x1,y1,z1)第三十二頁，共六十一頁，編輯于2023年，星期一海底曲面圖第三十三頁，共六十一頁，編輯于2023年，星期一%程序二：插值并作出水深小于5的海域范圍。x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,'v4');

%插值z1(z1>=5)=nan;

%將水深大于5的置為nan，這樣繪圖就不會顯示出來meshc(x1,y1,z1)

第三十四頁，共六十一頁，編輯于2023年，星期一水深小于5的海域范圍第三十五頁，共六十一頁，編輯于2023年，星期一實(shí)驗(yàn)作業(yè)1

山區(qū)地貌：在某山區(qū)測得一些地點(diǎn)的高程如下表：(平面區(qū)域1200≤x

≤4000,1200≤y

≤3600)，試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進(jìn)行比較．第三十六頁，共六十一頁，編輯于2023年，星期一5.1引言

對于情況較復(fù)雜的實(shí)際問題（因素不易化簡，作用機(jī)理不詳）可直接使用數(shù)據(jù)組建模，尋找簡單的因果變量之間的數(shù)量關(guān)系，從而對未知的情形作預(yù)報(bào)。這樣組建的模型為擬合模型。擬合模型的組建主要是處理好觀測數(shù)據(jù)的誤差，使用數(shù)學(xué)表達(dá)式從數(shù)量上近似因果變量之間的關(guān)系。擬合模型的組建是通過對有關(guān)變量的觀測數(shù)據(jù)的觀察、分析和選擇恰當(dāng)?shù)臄?shù)學(xué)表達(dá)方式得到的。五、擬合的使用及求解第三十七頁，共六十一頁，編輯于2023年，星期一5.2擬合模型的分類

5.2.1直線擬合5.2.2曲線擬合5.2.3觀察數(shù)據(jù)修勻?qū)τ谝呀o一批實(shí)測數(shù)據(jù)，由于實(shí)測方法、實(shí)驗(yàn)環(huán)境等一些外界因素的影響，不可避免地會產(chǎn)生隨機(jī)干擾和誤差。我們自然希望根據(jù)數(shù)據(jù)分布的總趨勢去剔除觀察數(shù)據(jù)中的偶然誤差，這就是所謂的數(shù)據(jù)修勻（或稱數(shù)據(jù)平滑）問題。

第三十八頁，共六十一頁，編輯于2023年，星期一直線擬合問題引例1溫度t(oC)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求60oC時(shí)的電阻R．

設(shè)

R=at+ba,b為待定系數(shù)第三十九頁，共六十一頁，編輯于2023年，星期一曲線擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時(shí)間的變化規(guī)律c(t).在直角坐標(biāo)系下作圖如下(plot)MATLAB(aa1)第四十頁，共六十一頁，編輯于2023年，星期一曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個(gè)點(diǎn)（xi,yi)i=1,…,n,

尋求一個(gè)函數(shù)（曲線）y=f(x),

使f(x)

在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好．

+++++++++xyy=f(x)(xi,yi)ii為點(diǎn)（xi,yi)與曲線y=f(x)的距離第四十一頁，共六十一頁，編輯于2023年，星期一曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(shù)

r1(x),r2(x),…,rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…,am

為待定系數(shù)．

第二步:確定a1,a2,…,am

的準(zhǔn)則（最小二乘準(zhǔn)則）：使n個(gè)點(diǎn)（xi,yi)與曲線y=f(x)的距離i

的平方和最小

．記

問題歸結(jié)為，求

a1,a2,…,am

使

J(a1,a2,…,am)

最?。谒氖摚擦豁?，編輯于2023年，星期一線性最小二乘法的求解：預(yù)備知識超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。

如果有向量a使得達(dá)到最小，則稱a為上述超定方程的最小二乘解。第四十三頁，共六十一頁，編輯于2023年，星期一線性最小二乘法的求解

定理：當(dāng)RTR可逆時(shí)，超定方程組（3）存在最小二乘解，且即為方程組

RTRa=RTy的解：a=(RTR)-1RTy

所以，曲線擬合的最小二乘法要解決的問題，實(shí)際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y（3）第四十四頁，共六十一頁，編輯于2023年，星期一線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數(shù){r1(x),…rm(x)}的選取

1.通過機(jī)理分析建立數(shù)學(xué)模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過直觀判斷確定f(x)：第四十五頁，共六十一頁，編輯于2023年，星期一用MATLAB作線性最小二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)輸入同長度的數(shù)組x，y擬合多項(xiàng)式次數(shù)2.多項(xiàng)式在x處的值y可用以下命令計(jì)算：

y=polyval（a，x）第四十六頁，共六十一頁，編輯于2023年，星期一1）輸入以下命令：

x=0:0.1:1;y=[-0.4471.9783.286.167.087.34…7.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2）計(jì)算結(jié)果：Ａ=-9.810820.1293-0.0317用多項(xiàng)式擬合的命令MATLAB(zxec)第四十七頁，共六十一頁，編輯于2023年，星期一如何預(yù)報(bào)人口的增長

人口的增長是當(dāng)前世界上引起普遍關(guān)注的問題，并且我們會發(fā)現(xiàn)在不同的刊物預(yù)報(bào)同一時(shí)間的人口數(shù)字不相同，這顯然是由于用了不同的人口模型計(jì)算的結(jié)果。我國是世界第一人口大國，基本上地球每九個(gè)人中就有一個(gè)中國人。有效地控制我國人口的增長是使我過全面進(jìn)入小康社會、到21世紀(jì)中葉建成富強(qiáng)民主文明的社會主義國家的需要。而有效控制人口增長的前提是要認(rèn)識人口數(shù)量的變化規(guī)律，建立人口模型，作出較準(zhǔn)確的預(yù)報(bào)。例:如何預(yù)報(bào)人口的增長第四十八頁，共六十一頁，編輯于2023年，星期一例如：1949年—1994年我國人口數(shù)據(jù)資料如下：年份xi1949195419591964196919741979198419891994人口數(shù)yi5.46.06.77.08.19.19.810.311.311.8建模分析我國人口增長的規(guī)律,預(yù)報(bào)1999年我國人口數(shù)。第四十九頁，共六十一頁，編輯于2023年，星期一模型一：假設(shè)人口隨時(shí)間線性地增加模型：參數(shù)估計(jì)觀測值的模型：

擬合的精度:誤差平方和?？梢运愠觯篴=-283.2320b=0.1480模型：y=–1.93+0.146x

第五十頁，共六十一頁，編輯于2023年，星期一則可看成是線性方程,用polyfit命令計(jì)算得：模型二：指數(shù)增長模型可變?yōu)閅A=+BXa=2.33,b=0.0179則所求模型為:第五十一頁，共六十一頁，編輯于2023年，星期一程序如下：x=[1949195419591964196919741979198419891994];y=[5.46.06.77.08.19.19.810.311.311.8];a=polyfit(x,y,1);x1=[1949:10:1994];y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=exp(b(2))*exp(b(1)*x1);plot(x,y,'*')holdonplot(x1,y1,'--r')holdonplot(x1,y2,'-k')legend('原曲線','模型一曲線','模型二曲線')第五十二頁，共六十一頁，編輯于2023年，星期一結(jié)論的比較如下表：年份xi1949195419591964196919741979198419891994人口數(shù)yi5.46.06.77.08.19.19.810.311.311.8模型一值5.245.976.707.438.168.909.6210.3611.0911.82誤差0.160.030.00-0.43-0.060.200.18-0.060.01-0.02模型二值5.556.066.627.237.908.649.4410.3111.2612.31誤差-0.15-0.060.08-0.230.200.460.36-0.01-0.13-0.51第五十三頁，共六十一頁，編輯于2023年，星期一結(jié)果分析：（1）Q1=0.2915<0.7437=Q2.線性模型更適合中國人口的增長。（2）預(yù)報(bào)：1999年12.55億，13.43億（3）統(tǒng)計(jì)年鑒：

2005年13.3億，2010年14億模型I2005年13.43億，2010年14.16億模型II2005年14.94億，2010年16.33億第五十四頁，共六十一頁，編輯于2023年，星期一1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1,xdata2,…,xdatan）,ydata=（ydata1,ydata2,…,ydatan）

用MATLAB作非線性最小二乘擬合

MATLAB提供了求非線性最小二乘擬合的函數(shù)：lsqcurvefit．這個(gè)命令都要先建立M文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x,xdata1）,…,F（x，xdatan））T中的參變量x(向量),使得

第五十五頁，共六十一頁，編輯于2023年，星期一

輸入格式為:

(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);

(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個(gè)事先建立的定義函數(shù)F(x,xdata)的M文件,自變量為x和xdata說明:x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見無約束優(yōu)化第五十六頁，共六十一頁，編輯于2023年，星期一