數(shù)學(xué)建模教程擬合與插值第1頁(yè)在大量應(yīng)用領(lǐng)域中，人們經(jīng)常面臨這么問(wèn)題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求曲線或曲面。對(duì)這個(gè)問(wèn)題有兩種方法。一個(gè)是插值法，數(shù)據(jù)假定是正確，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生情況。另一個(gè)方法是曲線擬合或回歸。人們?cè)O(shè)法找出某條光滑曲線，它最正確地?cái)M合數(shù)據(jù)，但無(wú)須要經(jīng)過(guò)任何數(shù)據(jù)點(diǎn)。本專題主要目標(biāo)是：了解插值和擬合基本內(nèi)容；掌握用Matlab求解插值與擬合問(wèn)題基本命令。

函數(shù)插值與曲線擬合都是要依據(jù)一組數(shù)據(jù)結(jié)構(gòu)一個(gè)函數(shù)作為近似，因?yàn)榻埔蟛灰粯?，二者?shù)學(xué)方法上是完全不一樣。第2頁(yè)內(nèi)容提要1.擬合問(wèn)題引例及基本理論2.Matlab求解擬合問(wèn)題3.應(yīng)用實(shí)例4.插值問(wèn)題引例及基本理論5.Maltab求解插值問(wèn)題6.應(yīng)用實(shí)例第3頁(yè)擬合問(wèn)題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時(shí)電阻R。

設(shè)

R=at+ba,b為待定系數(shù)一、擬合問(wèn)題第4頁(yè)擬合問(wèn)題引例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).作半對(duì)數(shù)坐標(biāo)系(semilogy)下圖形第5頁(yè)曲線擬合問(wèn)題提法已知一組（二維）數(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)

i

i為點(diǎn)（xi,yi)與曲線y=f(x)距離第6頁(yè)線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數(shù){r1(x),…rm(x)}選取

1.經(jīng)過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，經(jīng)過(guò)直觀判斷確定f(x)：第7頁(yè)曲線擬合問(wèn)題最慣用解法——線性最小二乘法基本思緒第一步:先選定一組函數(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平方和最小。記

問(wèn)題歸結(jié)為，求

a1,a2,…am

使J(a1,a2,…am)

最小。第8頁(yè)線性最小二乘法求解：預(yù)備知識(shí)超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)方程組即Ra=y其中超定方程普通是不存在解矛盾方程組。假如有向量a使得到達(dá)最小，則稱a為上述超定方程最小二乘解。第9頁(yè)線性最小二乘法求解定理：當(dāng)RTR可逆時(shí)，超定方程組（3）存在最小二乘解，

且即為方程組RTRa=RTy------正則（正規(guī)）方程組解：a=(RTR)-1RTy所以，曲線擬合最小二乘法要處理問(wèn)題，實(shí)際上就是求以下超定方程組最小二乘解問(wèn)題。其中Ra=y（3）第10頁(yè)用MATLAB解擬合問(wèn)題1、線性最小二乘擬合2、非線性最小二乘擬合第11頁(yè)用MATLAB作線性最小二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已經(jīng)有命令:a=polyfit(x,y,m)3.對(duì)超定方程組可得最小二乘意義下解。，用2.多項(xiàng)式在x處值y計(jì)算命令：y=polyval（a，x）輸出擬合多項(xiàng)式系數(shù)a=[a1,…,am,am+1]’（數(shù)組）輸入同長(zhǎng)度數(shù)組X，Y擬合多項(xiàng)式

次數(shù)第12頁(yè)即要求出二次多項(xiàng)式:中使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合第13頁(yè)1）輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];R=[(x.^2)',x',ones(11,1)]；

A=R\y'MATLAB(zxec1)解法1．解超定方程方法2）計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]第14頁(yè)第15頁(yè)2）計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]解法2．用多項(xiàng)式擬合命令MATLAB(zxec2)1）輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和擬合曲線圖形第16頁(yè)1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非線性最小二乘擬合兩個(gè)求非線性最小二乘擬合函數(shù)：lsqcurvefit、lsqnonlin。相同點(diǎn)和不一樣點(diǎn)：兩個(gè)命令都要先建立M-文件fun.m，定義函數(shù)f(x)，但定義f(x)方式不一樣。

lsqcurvefit用以求含參量x（向量）向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T使得第17頁(yè)輸入格式:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub);

(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,lb,ub,options);(4)[x,resnorm]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,resnorm,residual]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);

fun是一個(gè)事先建立定義函數(shù)F(x,xdata)M-文件,自變量為x和xdata說(shuō)明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化第18頁(yè)lsqnonlin用以求含參量x（向量）向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

，使得

最小。其中fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai2.lsqnonlin已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）第19頁(yè)輸入格式：1）x=lsqnonlin（‘fun’，x0）；2）x=lsqnonlin（‘fun’，x0，lb，ub）；3）x=lsqnonlin（‘fun’，x0，，lb，ub，options）；4）[x，resnorm]=lsqnonlin（‘fun’，x0，…）；5）[x，resnorm

，residual]=lsqnonlin（‘fun’，x0，…）；說(shuō)明：x=lsqnonlin（‘fun’，x0，options）；fun是一個(gè)事先建立定義函數(shù)f(x)M-文件，自變量為x迭代初值選項(xiàng)見(jiàn)無(wú)約束優(yōu)化第20頁(yè)例2用下面一組數(shù)據(jù)擬合中參數(shù)a，b，k該問(wèn)題即解最優(yōu)化問(wèn)題：第21頁(yè)

1）編寫(xiě)M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];

x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令lsqcurvefit第22頁(yè)3）運(yùn)算結(jié)果：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結(jié)論：a=0.0063,b=-0.0034,k=0.2542第23頁(yè)1）編寫(xiě)M-文件curvefun2.m

functionf=curvefun2(x)tdata=100:100:1000;cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函數(shù)curvefun2自變量是x，cdata和tdata是已知參數(shù)，故應(yīng)將cdatatdata值寫(xiě)在curvefun2.m中解法2用命令lsqnonlinx=（a，b，k）第24頁(yè)3）運(yùn)算結(jié)果為

f=1.0e-003*（0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792）x=0.0063-0.00340.2542能夠看出，兩個(gè)命令計(jì)算結(jié)果是相同。4）結(jié)論：即擬合得a=0.0063b=-0.0034k=0.2542第25頁(yè)插值問(wèn)題第26頁(yè)擬合與插值關(guān)系說(shuō)明：函數(shù)插值與曲線擬合都是要依據(jù)一組數(shù)據(jù)結(jié)構(gòu)一個(gè)函數(shù)作為近似，因?yàn)榻埔蟛灰粯?，二者?shù)學(xué)方法上完全不一樣。實(shí)例：下面數(shù)據(jù)是某次試驗(yàn)所得，希望得到x和f之間關(guān)系？MATLAB(cn)問(wèn)題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求曲線或曲面處理方案：若不要求曲線（面）經(jīng)過(guò)全部數(shù)據(jù)點(diǎn)，而是要求它反應(yīng)對(duì)象整體改變趨勢(shì)，就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）經(jīng)過(guò)所給全部數(shù)據(jù)點(diǎn)，就是插值問(wèn)題；第27頁(yè)最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：第28頁(yè)拉格朗日插值分段線性插值一維插值一、插值定義二、插值方法三、用Matlab解插值問(wèn)題第29頁(yè)返回二維插值一、二維插值定義二、網(wǎng)格節(jié)點(diǎn)插值法三、用Matlab解插值問(wèn)題最鄰近插值分片線性插值雙線性插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)插值散點(diǎn)數(shù)據(jù)插值第30頁(yè)一維插值定義已知n+1個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)求任一插值點(diǎn)處插值

第31頁(yè)結(jié)構(gòu)一個(gè)(相對(duì)簡(jiǎn)單)函數(shù)經(jīng)過(guò)全部節(jié)點(diǎn),即再用計(jì)算插值，即

返回第32頁(yè)拉格朗日插值問(wèn)題函數(shù)值

已知

y=f(x)在n+1個(gè)點(diǎn)結(jié)構(gòu)n次多項(xiàng)式

pn(x),使得從而得到

f(x)

近似計(jì)算式

第33頁(yè)一、線性插值（n=1，一次插值）求解L1(x)=a1

x+a0已知使得L1(xi)

=yi.(i=0,1)假如令則稱l0(x),l1(x)為x0,x1上線性插值基函數(shù)。這時(shí)依據(jù)點(diǎn)斜式得到并稱其為一次Lagrange插值多項(xiàng)式。f(x)

≈L1(x)=y0l0(x)+y1l1(x)第34頁(yè)二、拋物線插值（n=2，二次插值）求解L2(x)=a2x2+a1

x+a0使得L2(xi)=yi,i=0,1,2.關(guān)于二次多項(xiàng)式結(jié)構(gòu)采取以下方法：令已知并由插值條件得到L2(x)=A(x-x1)(x-x2)+B(x-x0)(x-x2)+C(x-x0)(x-x1)L2(x0)=y0,L2(x1)=y1

,L2(x)=y2第35頁(yè)于是得到這時(shí)：f(x)

≈L2(x)=y0l0(x)+y1l1(x)+y2l2(x)假如令則有并稱其為二次Lagrange插值多項(xiàng)式?！瘮?shù)表示x0,x1,x2上二次插值基函數(shù)第36頁(yè)稱為拉格朗日插值基函數(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.處理此問(wèn)題拉格朗日插值多項(xiàng)式公式以下其中Li(x)為n次多項(xiàng)式：拉格朗日(Lagrange)插值第37頁(yè)拉格朗日(Lagrange)插值尤其地:兩點(diǎn)一次(線性)插值多項(xiàng)式:三點(diǎn)二次(拋物)插值多項(xiàng)式:第38頁(yè)拉格朗日多項(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é)果圖形.例第39頁(yè)第40頁(yè)分段線性插值計(jì)算量與n無(wú)關(guān);n越大，誤差越小.

xjxj-1xj+1x0xnxoy第41頁(yè)例用分段線性插值法求插值,并觀察插值誤差.在[-6,6]中平均選取41個(gè)點(diǎn)作插值,結(jié)果如圖示第42頁(yè)比分段線性插值更光滑。

xyxi-1xiab在數(shù)學(xué)上，光滑程度定量描述是：函數(shù)(曲線)k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線含有k階光滑性。光滑性階次越高，則越光滑。是否存在較低次分段多項(xiàng)式到達(dá)較高階光滑性方法？三次樣條插值就是一個(gè)很好例子。三次樣條插值第43頁(yè)g(x)為被插值函數(shù)。三次樣條插值第44頁(yè)例用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值第45頁(yè)用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范圍。第46頁(yè)例：在1-1211小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試預(yù)計(jì)每隔1/10小時(shí)溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數(shù)據(jù)將是很多)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)第47頁(yè)第48頁(yè)

xy機(jī)翼下輪廓線例已知飛機(jī)下輪廓線上數(shù)據(jù)以下，求x每改變0.1時(shí)y值。第49頁(yè)第50頁(yè)二維插值定義

xyO第一個(gè)（網(wǎng)格節(jié)點(diǎn)）：第51頁(yè)

已知m

n個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)結(jié)構(gòu)一個(gè)二元函數(shù)經(jīng)過(guò)全部已知節(jié)點(diǎn),即再用計(jì)算插值，即第52頁(yè)第二種（散亂節(jié)點(diǎn)）：

yx0第53頁(yè)已知n個(gè)節(jié)點(diǎn)其中互不相同，結(jié)構(gòu)一個(gè)二元函數(shù)經(jīng)過(guò)全部已知節(jié)點(diǎn),即再用計(jì)算插值，即第54頁(yè)注意：最鄰近插值普通不連續(xù)。含有連續(xù)性最簡(jiǎn)單插值是分片線性插值。最鄰近插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O二維或高維情形最鄰近插值，與被插值點(diǎn)最鄰近節(jié)點(diǎn)函數(shù)值即為所求。第55頁(yè)將四個(gè)插值點(diǎn)（矩形四個(gè)頂點(diǎn)）處函數(shù)值依次簡(jiǎn)記為：

分片線性插值xy

(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第56頁(yè)插值函數(shù)為：第二片(上三角形區(qū)域)：(x,y)滿足插值函數(shù)為：注意：(x,y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)；分兩片函數(shù)表示式以下：第一片(下三角形區(qū)域)：(x,y)滿足第57頁(yè)雙線性插值是一片一片空間二次曲面組成。雙線性插值函數(shù)形式以下：其中有四個(gè)待定系數(shù)，利用該函數(shù)在矩形四個(gè)頂點(diǎn)（插值節(jié)點(diǎn)）函數(shù)值，得到四個(gè)代數(shù)方程，恰好確定四個(gè)系數(shù)。雙線性插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O第58頁(yè)

要求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í),雙線性插值第59頁(yè)例：測(cè)得平板表面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)畫(huà)出原始數(shù)據(jù)，畫(huà)出粗糙溫度分布曲圖.第60頁(yè)第61頁(yè)再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫(huà)出插值后溫度分布曲面圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位地方進(jìn)行插值.第62頁(yè)第63頁(yè)

經(jīng)過(guò)此例對(duì)最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法插值效果進(jìn)行比較。第64頁(yè)原始數(shù)據(jù)圖第65頁(yè)最鄰近插值第66頁(yè)雙線性插值第67頁(yè)雙三次插值第68頁(yè)等高線圖第69頁(yè)

插值函數(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í),雙線性插值第70頁(yè)一室模型：將整個(gè)機(jī)體看作一個(gè)房室，稱中心室，室內(nèi)血藥濃度是均勻?？焖凫o脈注射后，濃度馬上上升；然后快速下降。當(dāng)濃度太低時(shí)，達(dá)不到預(yù)期治療效果；當(dāng)濃度太高，又可能造成藥品中毒或副作用太強(qiáng)。臨床上，每種藥品有一個(gè)最小有效濃度c1和一個(gè)最大有效濃度c2。設(shè)計(jì)給藥方案時(shí)，要使血藥濃度保持在c1~c2之間。本題設(shè)c1=10，c2=25(ug/ml).擬合問(wèn)題實(shí)例給藥方案——一個(gè)新藥用于臨床之前，必須設(shè)計(jì)給藥方案.藥品進(jìn)入機(jī)體后血液輸送到全身，在這個(gè)過(guò)程中不停地被吸收、分布、代謝，最終排出體外，藥品在血液中濃度，即單位體積血液中藥品含量，稱為血藥濃度。第71頁(yè)

在試驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注入該藥品300mg后,在一定時(shí)刻t(小時(shí))采集血藥,測(cè)得血藥濃度c(ug/ml)以下表:

t(h)0.250.511.52346