數(shù)學(xué)建模插值與擬合第一頁，共五十一頁，2022年，8月28日我們經(jīng)常會遇到大量的數(shù)據(jù)需要處理，而處理數(shù)據(jù)的關(guān)鍵就在于這些算法，例如數(shù)據(jù)擬合、參數(shù)估計、插值等數(shù)據(jù)處理算法。此類問題在MATLAB中有很多現(xiàn)成的函數(shù)可以調(diào)用，熟悉MATLAB，這些方法都能游刃有余的用好。一、概述第二頁，共五十一頁，2022年，8月28日數(shù)據(jù)擬合在很多賽題中有應(yīng)用，與圖形處理有關(guān)的問題很多與插值和擬合有關(guān)系，例如98年美國賽A題，生物組織切片的三維插值處理，94年A題逢山開路，山體海拔高度的插值計算，2003年吵的沸沸揚揚的“非典”問題也要用到數(shù)據(jù)擬合算法，觀察數(shù)據(jù)的走向進行處理，2005年的雨量預(yù)報的評價的插值計算。2001年的公交車調(diào)度擬合問題，2003年的飲酒駕車擬合問題。第三頁，共五十一頁，2022年，8月28日插值問題——雨量預(yù)報的評價預(yù)測點和實測點的圖形插值后的圖形第四頁，共五十一頁，2022年，8月28日擬合問題——飲酒駕車喝兩瓶酒的擬合曲線喝1-5瓶酒的擬合曲線第五頁，共五十一頁，2022年，8月28日在實際中，常常要處理由實驗或測量所得到的一些離散數(shù)據(jù)。插值與擬合方法就是要通過這些數(shù)據(jù)去確定某一類已知函數(shù)的參數(shù)或?qū)で竽硞€近似函數(shù)，使所得到的近似函數(shù)與已知數(shù)據(jù)有較高的擬合精度。如果要求這個近似函數(shù)（曲線或曲面）經(jīng)過所已知的所有數(shù)據(jù)點，則稱此類問題為插值問題。（不需要函數(shù)表達式）二、基本概念第六頁，共五十一頁，2022年，8月28日如果不要求近似函數(shù)通過所有數(shù)據(jù)點，而是要求它能較好地反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法稱為數(shù)據(jù)擬合。（必須有函數(shù)表達式）近似函數(shù)不一定（曲線或曲面）通過所有的數(shù)據(jù)點。第七頁，共五十一頁，2022年，8月28日1、聯(lián)系都是根據(jù)實際中一組已知數(shù)據(jù)來構(gòu)造一個能夠反映數(shù)據(jù)變化規(guī)律的近似函數(shù)的方法。2、區(qū)別插值問題不一定得到近似函數(shù)的表達形式，僅通過插值方法找到未知點對應(yīng)的值。數(shù)據(jù)擬合要求得到一個具體的近似函數(shù)的表達式。

三、插值與擬合的區(qū)別和聯(lián)系第八頁，共五十一頁，2022年，8月28日四、插值的使用及求解當(dāng)數(shù)據(jù)量不夠，需要補充，且認定已有數(shù)據(jù)可信時,通常利用函數(shù)插值方法。實際問題當(dāng)中碰到的函數(shù)f(x)是各種各樣的，有的表達式很復(fù)雜，有的甚至給不出數(shù)學(xué)的式子，只提供了一些離散數(shù)據(jù)，警如，某些點上的函數(shù)值和導(dǎo)數(shù)值。4.1引言第九頁，共五十一頁，2022年，8月28日選用不同類型的插值函數(shù)，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值）（2）分段線性插值（3）Hermite（4）三次樣條插值。4.2插值方法第十頁，共五十一頁，2022年，8月28日Matlab實現(xiàn)：實現(xiàn)分段線性插值不需要編制函數(shù)程序，它自身提供了內(nèi)部的功能函數(shù)interp1(一維插值)interp2(二維)interp3(三維)intern(n維)4.3MATLAB實現(xiàn)插值第十一頁，共五十一頁，2022年，8月28日用MATLAB作插值計算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結(jié)果‘nearest’最鄰近插值；‘linear’

線性插值；‘spline’三次樣條插值；‘cubic’立方插值；缺省時分段線性插值．注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍．第十二頁，共五十一頁，2022年，8月28日例：從1點12點的11小時內(nèi)，每隔1小時測量一次溫度，測得的溫度的數(shù)值依次為：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,'k+',h,t,'b',hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius')第十三頁，共五十一頁，2022年，8月28日xy機翼下輪廓線例已知飛機下輪廓線上數(shù)據(jù)如下，求x每改變0.1時的y值．第十四頁，共五十一頁，2022年，8月28日x0=[035791112131415];y0=[01.21.72.02.12.01.81.21.01.6];x=0:0.1:15;y1=interp1(x0,y0,x,'cubic');y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');subplot(3,1,1)plot(x0,y0,'k+',x,y1,'r')gridtitle('cubic')subplot(3,1,2)plot(x0,y0,'k+',x,y2,'r')gridtitle('piecewiselinear')subplot(3,1,3)plot(x0,y0,'k+',x,y3,'r')gridtitle('spline')返回第十五頁，共五十一頁，2022年，8月28日

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

最鄰近插值；‘linear’

雙線性插值；‘cubic’

雙三次插值；缺省時雙線性插值.第十六頁，共五十一頁，2022年，8月28日例：測得平板表面3×5網(wǎng)格點處的溫度分別為：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個單位的地方進行插值.第十七頁，共五十一頁，2022年，8月28日再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.第十八頁，共五十一頁，2022年，8月28日

通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較．第十九頁，共五十一頁，2022年，8月28日x=1200:400:4000;y=1200:400:3600;z=[11301250128012301040900500700;...13201450142014001300700900850;...139015001500140090011001060950;...15001200110013501450120011501010;...15001200110015501600155013801070;...15001550160015501600160016001550;...1480150015501510143013001200980];

輸入原始數(shù)據(jù):第二十頁，共五十一頁，2022年，8月28日figure(1);meshz(x,y,z)xlabel('X'),ylabel('Y'),zlabel('Z')

xi=1200:50:4000;yi=1200:50:3600;

figure(2)z1i=interp2(x,y,z,xi,yi','nearest');surfc(xi,yi,z1i)xlabel('X'),ylabel('Y'),zlabel('Z')

figure(3)z2i=interp2(x,y,z,xi,yi');surfc(xi,yi,z2i)xlabel('X'),ylabel('Y'),zlabel('Z')

figure(4)z3i=interp2(x,y,z,xi,yi','cubic');surfc(xi,yi,z3i)xlabel('X'),ylabel('Y'),zlabel('Z')

figure(5)subplot(1,3,1),contour(xi,yi,z1i,10,'r');subplot(1,3,2),contour(xi,yi,z2i,10,'r');subplot(1,3,3),contour(xi,yi,z3i,10,'r');注:surfc函數(shù)功能在矩形區(qū)域內(nèi)顯示三維帶陰影曲面圖,且在曲面下面畫出等高線。meshc函數(shù)

生成帶等高線網(wǎng)線圖meshz函數(shù)

生成帶垂簾的網(wǎng)線圖數(shù)據(jù)插值與圖形繪制:第二十一頁，共五十一頁，2022年，8月28日

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

cz

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

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

雙線性插值‘cubic’

雙三次插值'v4'-MATLAB提供的插值方法缺省時,雙線性插值第二十二頁，共五十一頁，2022年，8月28日例在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）×（-50，150）里的哪些地方船要避免進入．第二十三頁，共五十一頁，2022年，8月28日4.作出水深小于5的海域范圍,即z=5的等高線.2.在矩形區(qū)域(75,200)×(-50,150)進行插值。1.輸入插值基點數(shù)據(jù)3.作海底曲面圖第二十四頁，共五十一頁，2022年，8月28日%程序一：插值并作海底曲面圖

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)第二十五頁，共五十一頁，2022年，8月28日海底曲面圖第二十六頁，共五十一頁，2022年，8月28日%程序二：插值并作出水深小于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)

第二十七頁，共五十一頁，2022年，8月28日水深小于5的海域范圍第二十八頁，共五十一頁，2022年，8月28日5.1引言

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

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

第三十頁，共五十一頁，2022年，8月28日直線擬合問題引例1溫度t(oC)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求60oC時的電阻R．

設(shè)

R=at+ba,b為待定系數(shù)第三十一頁，共五十一頁，2022年，8月28日曲線擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時間的變化規(guī)律c(t).在直角坐標(biāo)系下作圖如下(plot)第三十二頁，共五十一頁，2022年，8月28日曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個點（xi,yi)i=1,…,n,尋求一個函數(shù)（曲線）y=f(x),使f(x)在某種準(zhǔn)則下與所有數(shù)據(jù)點最為接近，即曲線擬合得最好．

+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離第三十三頁，共五十一頁，2022年，8月28日擬合與插值的關(guān)系

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似，由于近似的要求不同，二者在數(shù)學(xué)方法上是完全不同的．實例：下面數(shù)據(jù)是某次實驗所得，希望得到X和f之間的關(guān)系？MATLAB(cn)問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合．若要求所求曲線（面）通過所給所有數(shù)據(jù)點，就是插值問題；第三十四頁，共五十一頁，2022年，8月28日最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：第三十五頁，共五十一頁，2022年，8月28日曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(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個點（xi,yi)與曲線y=f(x)的距離i

的平方和最?。?/p>

問題歸結(jié)為，求

a1,a2,…,am

使

J(a1,a2,…,am)

最?。谌?，共五十一頁，2022年，8月28日用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)輸入同長度的數(shù)組x，y擬合多項式次數(shù)2.多項式在x處的值y可用以下命令計算：y=polyval（a，x）第三十七頁，共五十一頁，2022年，8月28日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ù)點和擬合曲線的圖形2）計算結(jié)果：Ａ=-9.810820.1293-0.0317用多項式擬合的命令第三十八頁，共五十一頁，2022年，8月28日如何預(yù)報人口的增長

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

擬合的精度:誤差平方和。可以算出：a=-283.2320b=0.1480模型：y=-283.2320

+0.1480x

第四十一頁，共五十一頁，2022年，8月28日則可看成是線性方程,用polyfit命令計算得：模型二：指數(shù)增長模型可變?yōu)閅A=+BXa=2.33,b=0.0179則所求模型為:第四十二頁，共五十一頁，2022年，8月28日程序如下：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('原曲線','模型一曲線','模型二曲線')第四十三頁，共五十一頁，2022年，8月28日結(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第四十四頁，共五十一頁，2022年，8月28日結(jié)果分析：（1）Q1=0.2915<0.7437=Q2.線性模型更適合中國人口的增長。（2）預(yù)報：1999年12.55億，13.43億（3）統(tǒng)計年鑒：2005年13.3億，2010年14億模型I2005年13.43億，2010年14.16億模型II2005年14.94億，2010年16.33億第四十五頁，共五十一頁，2022年，8月28日1.lsqcurvefit已知數(shù)據(jù)點：xdata=（xdata1,xdata2,…,xdatan）,ydata=（ydata1,ydata2,…,ydatan）

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

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

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

第四十六頁，共五十一頁，2022年，8月28日輸