數(shù)學實驗第七次講稿

數(shù)學實驗之七數(shù)據(jù)擬合2023/2/21實驗目的

[1]

了解最小二乘擬合的原理，掌握用MATLAB作線性最小二乘擬合的方法。[2]

通過實例學習如何用擬合方法解決實際問題；通過實例理解參數(shù)辨識問題的幾種方法。2023/2/22主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/23求電阻R隨溫度t的變化規(guī)律。已知熱敏電阻數(shù)據(jù)：溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032引例1：熱敏電阻電阻值的變化規(guī)律2023/2/24

設(shè)

R=at+ba,b為待定系數(shù)引例1：熱敏電阻電阻值的變化規(guī)律2023/2/25引例2：血藥濃度的變化規(guī)律

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).2023/2/26作圖觀察半對數(shù)坐標系(semilogy)下的圖形Log10c(t)=at+b2023/2/27主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/28曲線擬合的基本原理——問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個點（xi,yi)i=1,…n,

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

使f(x)

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

+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離2023/2/29問題的數(shù)學模型步驟：1)選定一類函數(shù)

f(x，a1，a2，

…，am)（1）其中

a1,a2,…am

為待定常數(shù)。2)確定參數(shù)a1,a2,…am，準則(最小二乘準則）：使n個點（xi,yi)

與曲線y=f(x，a1，a2，

…，am)的距離

i的平方和最小

。2023/2/210記問題歸結(jié)為:求

a1,a2,…am

使

J(a1,a2,…am)

最小。這樣的擬合稱為最小二乘擬合。問題的數(shù)學模型2023/2/211

特別，若選定一組函數(shù)

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

f(x,a1,a2,…am

)=a1r1(x)+a2r2(x)+…+amrm(x)其中

a1,a2,…am

為待定系數(shù)。問題的解法——線性最小二乘法

當f具有上述形式，且按最小二乘準則確定待定系數(shù)時，稱這樣的擬合為線性最小二乘擬合，稱對應(yīng)的函數(shù)f(x,a1,a2,…am

)

為已知數(shù)據(jù)點的線性最小二乘擬合函數(shù)。2023/2/2121)除了最小二乘準則（即各點誤差的平方和最?。阏J為還可以用怎樣的擬合準則？

2)比較起來，最小二乘準則有什么優(yōu)點？2023/2/213主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/214最小二乘擬合函數(shù)的選取+++++++++++++++f=a1+a2x將數(shù)據(jù)(xi,yi)i=1,…n作圖，通

過直觀判斷確定f：f=a1+a2x+a3x2f=a1+a2x+a3x22023/2/215最小二乘擬合函數(shù)的選取2.通過機理分析建立數(shù)學模型來確定f。+++++++++++++++f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)2023/2/216主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/217多項式擬合：作多項式f(x)=a1xm+…+amx+am+1函數(shù)擬合,可利用已有程序polyfit,其調(diào)用格式為:a=polyfit(x,y,m)用MATLAB作最小二乘擬合數(shù)據(jù)點擬合多項式次數(shù)系數(shù)2023/2/218例1.由數(shù)據(jù)溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032擬合R=f(t)用MATLAB作多項式擬合2023/2/2192.用命令polyfit(x,y,m)作最小二乘擬合3.編寫MATLAB程序dianzu1.m，并運行得到：a1=3.3940,a2=702.49181.選取擬合函數(shù)R=

a1t+a2用MATLAB作多項式擬合2023/2/220t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')用MATLAB作多項式擬合2023/2/2212.曲線擬合：

作一般的最小二乘曲線擬合，可利用已有程序lsqcurvefit,其調(diào)用格式為：

[a,resnorm,residual]=lsqcurvefit(‘f’,a0,x,y)

用MATLAB作曲線擬合數(shù)據(jù)點待定常數(shù)a的初值函數(shù)M文件擬合函數(shù)y=f(a,x)的函數(shù)M—文件誤差平方和各數(shù)據(jù)點處的誤差2023/2/222用MATLAB作最小二乘曲線擬合xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]例2：用函數(shù)

y(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合下列數(shù)據(jù)點：1.用命令lsqcurvefit(‘fun’,a0,x,y)2023/2/223a0=[1,1,1,0];xdata=[0:.1:2];%ydata省略;[a,resnorm,residual,flag,output]=lsqcurvefit('fitf',a0,xdata,ydata)xi=linspace(0,2,200);yi=fitf(a,xi);plot(xdata,ydata,'ro',xi,yi)xlabel('x'),ylabel('y=f(x)'),title('nonlinearcurvefitting')用MATLAB作最小二乘曲線擬合2023/2/224functiony=fitf(x,xdata)y=x(1)*exp(-x(2)*xdata)+x(3)*exp(-x(4)*xdata);主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/226范例1：靜脈注射的給藥方案問題背景

一種新藥用于臨床之前，必須設(shè)計給藥方案。在快速靜脈注射下，所謂給藥方案是指，每次注射量多大，間隔時間多長。藥物進入肌體后隨血液輸送到全身，在這過程中不斷被吸收、分解、代謝，最終排出體外。藥物向體外排出的速率與血藥濃度成正比。單位體積血液中的藥物含量，稱血藥濃度。臨床上，每種藥物有一個最小有效濃度c1和最大治療濃度c2。2023/2/227范例1：靜脈注射的給藥方案1.藥物排除速率與血藥濃度成正比，比例系數(shù)

k(>0);2.血液容積v,t=0瞬時注射劑量d,血藥濃度立即為d/v.模型假設(shè)：由假設(shè)1，3.快速靜脈注射下一室模型的血藥濃度：c(t)ctc00由假設(shè)2，2023/2/228若c1=10,c2=25(g/ml),

給藥方案設(shè)計歸結(jié)為根據(jù)數(shù)據(jù)(ti,ci)i=1,…n(d給定）擬合曲線c(t),以確定系數(shù)k,v.cc2c10t給藥方案設(shè)計：D0:初次劑量；:注射時間間隔；D:重復注射劑量。范例1：靜脈注射的給藥方案2023/2/229給藥方案血藥濃度數(shù)據(jù)(t=0注射300mg)

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01記則思考:取對數(shù)化為線性最小二乘,對結(jié)果有什么影響？2023/2/230主要內(nèi)容范例2：薄膜滲透率的測定布置實驗引例1,引例2擬合的基本原理擬合函數(shù)的選取用MATLAB作擬合計算范例1：靜脈注射的給藥方案2023/2/231范例2：薄膜滲透率的測定[1]問題背景；[2]假設(shè)；[3]問題的分析[4]數(shù)學模型；[5]參數(shù)辨識方法：

函數(shù)擬合法

非線性規(guī)劃法；

導函數(shù)擬合法；線性化迭代法。[6]結(jié)果分析VAVBS2023/2/232

某種醫(yī)用薄膜在試制時需測定其被物質(zhì)分子穿透的能力。測定方法：用面積為S的薄膜將容器分成兩部份，在兩部分中分別注滿該物質(zhì)的兩種不同濃度的溶液。此時該物質(zhì)分子就會從一側(cè)向另一側(cè)擴散。平均每單位時間通過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K稱為滲透率。定時測量容器中薄膜某一側(cè)的溶液濃度，以此確定K。VAVBS問題背景薄膜滲透率的測定2023/2/2333）薄膜是雙向同性的即物質(zhì)從膜的任何一側(cè)向另一側(cè)滲透的性能是相同的。VAVBS假設(shè)：1）薄膜兩側(cè)的溶液始終是均勻的；2）平均每單位時間通過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比。2023/2/234VAVBS

第一步：通過機理分析確定容器一側(cè)的濃度隨時間的變化規(guī)律：

C=CB(t)

第二步：利用實驗數(shù)據(jù)(tj,CB(tj))，和函數(shù)擬合的方法求出其中的未知參數(shù)，包括滲透率K.解決問題的思路：2023/2/235

考察時段[t，t+Δt]薄膜的一側(cè)容器中該物質(zhì)質(zhì)量的變化。

1）以容器A側(cè)為例，在時段[t，t+Δt]物質(zhì)質(zhì)量的增量為：

問題分析：

由假設(shè)2,在時段[t，t+Δt]，從B側(cè)滲透至A側(cè)的該物質(zhì)的質(zhì)量為：VAVBS2023/2/236于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：（1）問題分析：2023/2/237分別表示在初始時刻兩側(cè)溶其中2)注意到任意時刻，整個溶液中含有該物質(zhì)的質(zhì)量,與初始時刻該物質(zhì)的含量相同，因此問題分析：液的濃度2023/2/238從而：加上初值條件：代入式（1）得：便可得出CB(t)的變化規(guī)律，從而根據(jù)實驗數(shù)據(jù)進行擬合，估計出參數(shù)K,。問題分析：2023/2/239基于假設(shè)和前面的分析，B側(cè)的濃度CB(t)應(yīng)滿足如下微分方程和初始條件：數(shù)學模型2023/2/240模型求解方法

前面得到的模型是一個帶初值的一階線性微分方程，解之得：

問題歸結(jié)為利用CB在時刻tj的測量數(shù)據(jù)Cj(j=1,2,...,N)來辨識K和。1.函數(shù)擬合法2023/2/241引入從而

用函數(shù)CB(t)來擬合所給的實驗數(shù)據(jù)，從而估計出其中的參數(shù)a,b,K。若將其代入上式有：擬合函數(shù)化簡:2023/2/242用MATLAB軟件進行計算.1）編寫函數(shù)M-文件nongdu.mfunctionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中x(1)=a;x(2)=b;x(3)=k;2)編寫M文件(baomo.m)tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.10...6.266.396.506.59];x0=[0.2,0.05,0.05];x=lsqcurvefit(‘nongdu’,x0,tdata,cdata)編寫MATLAB程序:2023/2/2433)輸出結(jié)果:x=0.007-0.0030.1012

即k=0.1012,a=0.007,b=-0.003,結(jié)果:進一步求得：ToMATLAB(nongdu,baomo1)2023/2/2442.非線性規(guī)劃法

利用CB在時刻tj的測量數(shù)據(jù)Cj(j=1,2,...,N)來辨識K和。問題可轉(zhuǎn)化為求函數(shù)即求函數(shù)的最小值點（K，a，b）。2023/2/2453.導函數(shù)擬合法令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1,2,...,N)若知道一組數(shù)據(jù)則可用最小二乘擬合的方法來求出函數(shù)中的未知參數(shù)K和h。2023/2/246即為求參數(shù)K,a使下列誤差函數(shù)達到最?。?/p>

該問題等價于用函數(shù)f(K,a,CB)=K(0.01a-0.02CB)來擬合數(shù)據(jù)(j=1,2,...,N)用上述方法求出參數(shù)K,a.?3.導函數(shù)擬合法2023/2/2474.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中：

如果得到了參數(shù)K的一個較好的近似值K*，則將CB(t)關(guān)于K在K*處展開，略去K的二次及以上的項得CB(t)的一個近似式2023/2/2484.線性化迭代法通過極小化確定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一個新的近似值，用同樣的方法再求新的修正值K。這個過程可以不斷重復，直到修正值足夠小為止。2023/2/2491）當K的初值取為k=0.3時，出現(xiàn)奇異情況，迭代不收斂；2）當K的初值取為k=0.2時，經(jīng)四次迭代，已經(jīng)收斂到一個很好的解。迭代結(jié)果如下表。4.線性化迭代法2023/2/2503）取K的初值為0.1009,只一次迭代就得到2）中的最后結(jié)果。4.線性化迭代法導函數(shù)擬合法得出的參數(shù)值精度有限，線性化迭代法要求參數(shù)的初值比較接近精確值。可將導函數(shù)擬合法和線性化迭代法結(jié)合起來，把前者得到的參數(shù)K的值作為迭代法中K的初值，可使其收斂或收斂更快。2023/2/251結(jié)果及誤差分析

幾種方法得出的結(jié)果及相應(yīng)的誤差總結(jié)于下表，誤差為計算數(shù)據(jù)與實驗數(shù)據(jù)之差的平方和。2023/2/252％作出擬合函數(shù)曲線和數(shù)據(jù)點的圖形tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...

610626639650659];x0=[0.2,0.05,0.05];[x,resnorm,residual]=lsqcurvefit('nongdu',x0,tdata,cdata)t=linspace(100,1000,100);c=nongdu(x,t);plot(tdata,cdata,'o',t,c)％求擬合在各節(jié)點的誤差平方和c1=nongdu(x,tdata);e=c1-cdata;e1=sum(e.*e)函數(shù)擬合法的誤差分析程序2023/2/253函數(shù)擬合法的擬合效果返回2023/2/254人口問題馬爾薩斯人口模型假設(shè):人口出生率和死亡率是常數(shù)，因此人口的凈增長率為常數(shù).設(shè)時刻t的人口數(shù)量為p(t)，人口出生率為b，死亡率為d,(k=b-d)因此馬爾薩斯人口模型如下：該微分方程初值問題的解析解：馬爾薩斯人口模型為了驗證，采用某國家的人口歷史數(shù)據(jù),見下頁表.年份17901800181018201830184018501860187018801890人口/103392953087240963812866170692319231433385585015662948年份19001910192019301940195019601970198019902000人口/1037599591972105711122755131669150697179323203212226505248710281416首先確定參數(shù)k：因為所以利用1790年和1840年的數(shù)據(jù)計算得出：=0.0294.馬爾薩斯人口模型由此預測，1850年的人口數(shù)量為22,898,000，誤差為1%，1900年的人口數(shù)量為99,476,000，誤差為31%，2000年的人口數(shù)量為1,877,463,000，誤差為567%，2050年將達到80多個億？問題：（1）能不能用盡可能多的數(shù)據(jù)擬合微分方程的參數(shù)？（2）初始數(shù)據(jù)能不能也作為參數(shù)處理？這樣會得到更好的結(jié)果嗎？馬爾薩斯人口模型2023/2/259處理方法方法1：直接使用原始數(shù)據(jù)（t和p），使用非線性擬合求解方程的系數(shù)，并和原始數(shù)據(jù)作比較。方法2：因為，也就是lnp(t)=k(t-t0)+lnp0。記q(t)=lnp(t),q0=lnp0。上面的問題變成一個線性回歸問題，可以得到相應(yīng)的系數(shù)，和上面的方法已經(jīng)原始數(shù)據(jù)比較。方法3：還是使用方法1，僅僅對時間做一個變換，比如1790—0,1800—1，等等。重新使用方法1進行計算。1837年荷蘭生物數(shù)學家Verhulst提出改進.人口增長率不應(yīng)該是常數(shù)，假設(shè)增長率k是隨著人口數(shù)量接近最大數(shù)量M而線性遞減.邏輯斯諦人口模型

從而得到改進后的人口模型(邏輯斯諦增長模型)邏輯斯諦人口模型可以求得，該人口模型的解為：在這個方程中，將M，r，p0作為參數(shù)，使用前面類似的方法計算。并比較擬合數(shù)據(jù)和原始數(shù)據(jù)。掌握用MATLAB作最小二乘多項式擬合和曲線擬合的方法；2.通過實例學習如何用擬合方法解決實際問題，注意與插值方法的區(qū)別。3.鼓勵不囿于固定的模式或秩序，靈活調(diào)整思路，突破思維的呆板性，找到打破常規(guī)的解法。并在文獻檢索、動手和動腦等方面得到鍛煉，樹立創(chuàng)新意識。?布置實驗實驗目的2023/2/262實驗內(nèi)容?1.下表給出了近兩個世紀某國人口的統(tǒng)計數(shù)據(jù)，試預測2010年該國的人口。人口預測問題2023/2/2631)可先用以上數(shù)據(jù)擬合Multhus人口指數(shù)增長模型，根據(jù)檢驗結(jié)果進一步討論馬爾薩斯人口模型的改進。2)Malthus模型的基本假設(shè)是：人口的增長率為常數(shù)，記為r。記時刻t的人口為x(t），且初始時刻的人口為x0，于是得到如下微分方程：2023/2/2643)Malthus模型導致人口總數(shù)的劇烈上升，為此對該模型做一些改進。其中一種改進是認為增長率應(yīng)該是x的函數(shù)，并且隨著x增大出現(xiàn)資源的競爭導致增長率減小。因此一般可以假設(shè)r(x)=r-x/k。得到如下微分方程：2023/2/265實驗內(nèi)容?2.某年美國舊車價格的調(diào)查資料如下表其中xi表示轎車的使用年數(shù)，yi表示相應(yīng)的平均價格。試分析用什么形式的曲線來擬合上述的數(shù)據(jù)，并預測使用4.5年后轎車的平均價格大致為多少？舊車價格預測xi12345678910yi26151943149410877655384842902262042023/2/266?

經(jīng)濟增長模型3.增加生產(chǎn)、發(fā)展經(jīng)濟的主要因素有增加投資、勞動力以及技術(shù)革新等，在研究國民經(jīng)濟產(chǎn)值與這些因素的數(shù)量關(guān)系時，由于技術(shù)水平不像資金、勞動力那樣容易定量化，作為初步的模型，可認為技術(shù)水平不變，只討論產(chǎn)值和資金、勞動力之間的關(guān)系。用Q，K，L分別表示產(chǎn)值、資金、勞動力，要尋求Q(K,L)。經(jīng)過簡化與分析，在經(jīng)濟學中，推導出一個著名的Cobb-Douglas生產(chǎn)函數(shù)：Q(K,L)=aKαLβ，0<α,β<1（*）式中α,β，a要由經(jīng)濟統(tǒng)計數(shù)據(jù)確定。根據(jù)下表所給的統(tǒng)計數(shù)據(jù)，估計α,β，a的值。

2023/2/267?

經(jīng)濟增長模型2023/2/268

tQK