數(shù)模培訓(xùn)-數(shù)據(jù)擬合方法

數(shù)據(jù)擬合與線(xiàn)性最小二乘法教學(xué)內(nèi)容數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題的求解思緒線(xiàn)性最小二乘法原理算例MATLAB工具箱Curvefit演示設(shè)R=at+ba,b為待定系數(shù)求電阻R隨溫度t的變化規(guī)律。知熱敏電阻數(shù)據(jù)：溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032引例1：熱敏電阻電阻值的變化規(guī)律t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01對(duì)某人用快速靜脈注射方式一次性注射某種藥物300mg后，經(jīng)過(guò)時(shí)間t采集血樣，測(cè)得血藥濃度c如下表：求血藥濃度隨時(shí)間的變化規(guī)律c(t).半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形Log10c(t)=at+b引例2：血藥濃度的變化規(guī)律數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題：知一維〔二維，…〕數(shù)據(jù)，即平面上的n個(gè)點(diǎn)(xi，yi)，i=1,2,…,n，xi互不一樣，尋求一個(gè)函數(shù)(曲線(xiàn))y=f(x)，使f(x)在某種準(zhǔn)那么下與一切數(shù)據(jù)點(diǎn)最為接近，即曲線(xiàn)擬合的最好，如以下圖所示(圖中δi為(xi，yi)與y=f(x)的間隔)。Oxyδi(xi，yi)數(shù)據(jù)擬合問(wèn)題的求解思緒線(xiàn)性最小二乘法是處理數(shù)據(jù)擬合最常用的方法。根本思緒：令f(x)=a1r1(x)+a2r2(x)+…+amrm(x)(1)其中rk(x)是事先選定的一組函數(shù)，ak是待定系數(shù)(k=1,2,…,m，m<n)。擬合準(zhǔn)那么是使n個(gè)點(diǎn)(xi，yi)，i=1,2,…,n，與y=f(xi)的間隔δi的平方和最小，稱(chēng)為最小二乘準(zhǔn)那么。擬合準(zhǔn)那么還有如最小一乘準(zhǔn)那么、極大極小準(zhǔn)那么等。線(xiàn)性最小二乘法原理1.實(shí)際———根本實(shí)際之a(chǎn)k確實(shí)定根據(jù)最小二乘準(zhǔn)那么，記J(a1，a2，…，am)=為求a1，a2，…，am是J到達(dá)最小，只需求利用極值的必要條件得到關(guān)于a1，…，am的線(xiàn)性方程組線(xiàn)性最小二乘法原理記，A=(a1，a2，…，am)T，y=(y1，…，yn)T，方程組(3)可表為RTRA=RTy(4)(4)稱(chēng)為法方程組，當(dāng){r1(x)，…，rm(x)}線(xiàn)性無(wú)關(guān)時(shí)，R列滿(mǎn)秩，RTR可逆，于是方程組(4)有獨(dú)一解A=(RTR)-1RTy(5)可以看出，只需f(x)關(guān)于待定系數(shù)a1，…，am線(xiàn)性，在最小二乘準(zhǔn)那么〔2〕下得到的方程組(3)關(guān)于a1，a2，…，am也一定是線(xiàn)性的，故稱(chēng)線(xiàn)性最小二乘法。線(xiàn)性最小二乘法原理2.實(shí)際______函數(shù)rk(x)的選取對(duì)數(shù)據(jù)(xi，yi)用線(xiàn)性最小二乘法作擬合時(shí)，首要的、也是關(guān)鍵的一步是恰當(dāng)?shù)剡x取r1(x)，r2(x)，…，rm(x)。n

假設(shè)經(jīng)過(guò)機(jī)理分析，可以知道y與x之間應(yīng)該有什么樣的函數(shù)關(guān)系，那么r1(x)，…，rm(x)容易確定。n

假設(shè)無(wú)法知道y與x之間的關(guān)系，可以將數(shù)據(jù)(xi，yi)，i=1,2,…,n作圖，直觀(guān)地判別應(yīng)該用什么樣的曲線(xiàn)去作擬合。常用的曲線(xiàn)有u

直線(xiàn)y=a1x+a2u

多項(xiàng)式y(tǒng)=a1xm+…+amx+am+1(普通m=2,3,不宜過(guò)高)u

雙曲線(xiàn)〔一支〕y=a1/x+a2u

指數(shù)曲線(xiàn)：擬合前需作變量代換，化為線(xiàn)性函數(shù)。對(duì)知數(shù)據(jù)，用什么樣的曲線(xiàn)擬合最好，可以在直觀(guān)判別的根底上，選擇幾種曲線(xiàn)分別作擬合，然后比較，看那條曲線(xiàn)的最小二乘目的J最小。線(xiàn)性最小二乘法原理3.求解方法l

描出數(shù)據(jù)的圖示；l

察看并選擇不同的數(shù)學(xué)函數(shù)進(jìn)展擬合；l

比較多種擬合結(jié)果，選擇其中較好的一種或者某幾種作為備選結(jié)果；注·意：通常需求將非線(xiàn)性函數(shù)rk(x)的轉(zhuǎn)化成線(xiàn)性的函數(shù)Rk(x)，然后再用Rk(x)進(jìn)展擬合，計(jì)算中通常需求列下表：i01…nxix0x1…xnyi=f(xi)y0y1…ynR1(x)R1(x0)R1(x1)…R1(xn)………………Rm(x)Rm(x0)Rm(x1)…Rm(xn)這樣就容易確定出法方程組RTRA=RTy。上表中后面的m行即為RT。算例【例】給定數(shù)據(jù)(xi，f(xi))，i=0,1,2,3,4，見(jiàn)下表，使選擇適當(dāng)?shù)哪Ｐ?，求最小二乘擬合函數(shù)g(x)。i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46Yi=lnf(xi)1.6291.7561.8762.0082.【解】：〔1〕、先描出數(shù)據(jù)的圖示算例〔2〕選定不同的數(shù)學(xué)函數(shù)〔模型〕或者rk(x)進(jìn)展擬合

l

直線(xiàn)模型y=a+bx選取y=a+bx，此時(shí)，r1(x)=1，r2(x)=x。要求y=a+bx與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，yi=f(xi)。列表計(jì)算如下：i01234xi1.001.251.501.752.00yi=f(xi)5.105.796.537.458.46r1(x)11111r2(x)1.001.251.501.752.00求解法方程組得到a=1.6380，b=3.3520，于是得到該模型下的最小二乘擬合曲線(xiàn)為g(x)=1.6380+3.3520x。算例算例

l

多項(xiàng)式模型y=a0+a1x+a2x2選取Y=a+bx+cx2，此時(shí)，r1(x)=1，r2(x)=x，r3(x)=x2。要求y=a+bx+cx2與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，yi=f(xi)。列表計(jì)算如下：i01234xi1.001.251.501.752.00Yi=f(xi)5.105.796.537.458.46r1(x)=111111r2(x)=x1.001.251.501.752.00r3(x)=x21.001.56252.253.06254.00算例求解法方程組得到a=3.6294，b=0.5406，c=0.9371，于是得到該模型下的最小二乘擬合曲線(xiàn)為g(x)=3.6294+0.5406x+0.9371x2。算例

l

雙曲線(xiàn)模型y=1/(a0+a1x)選取y=1/(a0+a1x)，令Y=1/y=a0+a1x，此時(shí)，r1(x)=1，r2(x)=x。要求Y=a0+a1x與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，Yi=1/f(xi)。列表計(jì)算如下：i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46Yi=1/f(xi)0.1960780.1727120.1530.1342280.118203r1(x)11111r2(x)1.001.251.501.752.00算例求解法方程組得到a0=0.27，a1=-0.07768，于是得到該模型下的最小二乘擬合曲線(xiàn)為g(x)=1/(0.27-0.07768x)。算例

l

指數(shù)曲線(xiàn)模型y=aebx選擇y=aebx，取對(duì)數(shù)lny=lna+bx，令Y=lny，A=lna，取r1(x)=1，r2(x)=x，要求Y=A+bx與(xi，Yi)，i=0,1,2,3,4，做最小二乘擬合，Yi=lnf(xi)。計(jì)算結(jié)果如下：i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46Yi=lnf(xi)1.6292405401.7561322921.876406942.008214032.34917r1(x)11111r2(x)1.001.251.501.752.00算例求解法方程組得到A=1.1225，b=0.5057，a=eA=3.072525924，于是得到此模型下的最小二乘擬合曲線(xiàn)g(x)=3.072525924e0.5057x。算例(3)比較上面的結(jié)果如下:i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46直線(xiàn)模型4.995.836.677.508.340.49e-1多項(xiàng)式模型5.1075.7696.5497.4458.460.85e-3雙曲線(xiàn)模型5.1625.7386.4577.3838.6180.42e-1指數(shù)函數(shù)模型5.0955.7816.5607.4448.4480.12e-2經(jīng)過(guò)比較，運(yùn)用多項(xiàng)式與指數(shù)函數(shù)擬合效果較好。用Matlab進(jìn)展數(shù)據(jù)擬合1.多項(xiàng)式曲線(xiàn)擬合:polyfit.y0=polyval(p,x0)p=polyfit(x,y,m)其中,x,y為知數(shù)據(jù)點(diǎn)向量,分別表示橫,縱坐標(biāo),m為擬合多項(xiàng)式的次數(shù),結(jié)果前往m次擬合多項(xiàng)式系數(shù),從高次到低次存放在向量p中.可求得多項(xiàng)式在x0處的值y0.例1知觀(guān)測(cè)數(shù)據(jù)點(diǎn)如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分別用3次和6次多項(xiàng)式曲線(xiàn)擬合這些數(shù)據(jù)點(diǎn).x=0:0.1:1y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)編寫(xiě)Matlab程序如下:t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)holdonplot(t,s,'r-','linewidth',2)plot(t,s,'b--','linewidth',2)gridx=0:0.1:1y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)例2用切削機(jī)床進(jìn)展金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需求測(cè)定刀具的磨損速度.在一定的時(shí)間丈量刀具的厚度,得數(shù)據(jù)如表所示:切削時(shí)間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時(shí)間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')a=-0.301229.3804holdonplot(t,y1),holdoffa=polyfit(t,y,1)y1=-0.3012*t+29.3804例2用切削機(jī)床進(jìn)展金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需求測(cè)定刀具的磨損速度.在一定的時(shí)間丈量刀具的厚度,得數(shù)據(jù)如表所示:切削時(shí)間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時(shí)間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm擬合曲線(xiàn)為:y=-0.3012t+29.3804例3一個(gè)15.4cm×30.48cm的混凝土柱在加壓實(shí)驗(yàn)中的應(yīng)力-應(yīng)變關(guān)系測(cè)試點(diǎn)的數(shù)據(jù)如表所示1.552.472.933.03知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描畫(huà),即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.2.89知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描畫(huà),即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.解選取指數(shù)函數(shù)作擬合時(shí),在擬合前需作變量代換,化為k1,k2的線(xiàn)性函數(shù).于是,令即在命令窗口輸入:x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6]y=[3.103*1.0e+32.465*1.0e+31.953*1.0e+31.517*1.0e+31.219*1.0e+3]z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=[1.552.472.933.032.89]plot(x,w,'*')y1=exp(8.3009)*x.*exp(-494.5209*x)plot(x,w,'*',x,y1,'r-')知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描畫(huà),即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.擬合曲線(xiàn)為:令那么求得于是在實(shí)踐運(yùn)用中常見(jiàn)的擬合曲線(xiàn)有:直線(xiàn)多項(xiàng)式普通n=2,3,不宜過(guò)高.雙曲線(xiàn)(一支)指數(shù)曲線(xiàn)2.非線(xiàn)性曲線(xiàn)擬合:lsqcurvefit.功能:x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)根據(jù)給定的數(shù)據(jù)xdata,ydata(對(duì)應(yīng)點(diǎn)的橫,縱坐標(biāo)),按函數(shù)文件fun給定的函數(shù),以x0為初值作最小二乘擬合,前往函數(shù)fun中的系數(shù)向量x和殘差的平方和resnorm.例4知觀(guān)測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線(xiàn)f(x)=aex+bx2+cx3與知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.首先編寫(xiě)存儲(chǔ)擬合函數(shù)的函數(shù)文件.functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.^2+x(3)*xdata.^3保管為文件nihehanshu.m例4知觀(guān)測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線(xiàn)f(x)=aex+bx2+cx3與知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.編寫(xiě)下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;ydata=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];x0=[0,0,0];[x,resnorm]=lsqcurvefit(@nihehanshu,x0,xdata,ydata)編寫(xiě)下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;ydata=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];x0=[0,0,0];[x,resnorm]=lsqcurvefit(@nihehanshu,x0,xdata,ydata)程序運(yùn)轉(zhuǎn)后顯示x=3.00224.03040.9404resnorm=0.0912例4知觀(guān)測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線(xiàn)f(x)=aex+bx2+cx3與知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.闡明:最小二乘意義上的最正確擬合函數(shù)為f(x)=3ex+4.03x2+0.94x3.此時(shí)的殘差是:0.0912.f(x)=3ex+4.03x2+0.94x3.擬合函數(shù)為:作業(yè)：用函數(shù)f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合以下數(shù)據(jù)點(diǎn)：xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.91.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.00.22110.17040.2636]用命令lsqcurvefit(‘f’,a0,x,y)作業(yè)：調(diào)查某種纖維的強(qiáng)度與其拉伸倍數(shù)的關(guān)系,下表是實(shí)踐測(cè)定的24個(gè)纖維樣品的強(qiáng)度與相應(yīng)的拉伸倍數(shù)的記錄:試建立纖維強(qiáng)度與拉伸倍數(shù)之間的關(guān)系。作業(yè)：煉鋼廠(chǎng)出鋼時(shí)所用盛鋼水的鋼包，由于鋼水對(duì)耐火資料的侵蝕，容積不斷增大，我們希望找出運(yùn)用次數(shù)與增大容積之間的函數(shù)關(guān)系.實(shí)驗(yàn)數(shù)據(jù)如下：表4.2鋼包運(yùn)用次數(shù)與增大容積使用次數(shù)23456789增大容積6.428.29.589.59.7109.939.99使用次數(shù)10111213141516增大容積10.4910.5910.610.810.610.910.76分別選擇函數(shù)擬合鋼包容積與運(yùn)用次數(shù)的關(guān)系,在同一坐標(biāo)系內(nèi)作出函數(shù)圖形.x1=[2:16];y1=[6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76];b01=[0.1435,0.084];%初始參數(shù)值fun1=inline('x./(b(1)+b(2)*x)','b','x');%定義函數(shù)[b1,r1,j1]=nlinfit(x1,y1,fun1,b01);y=x1./(0.1152+0.0845*x1);%根據(jù)b1寫(xiě)出詳細(xì)函數(shù)plot(x1,y1,'*',x1,y,'-or');下面給出分式函數(shù)擬合程序：初始參數(shù)b0的計(jì)算，由于確定兩個(gè)參數(shù)值，因此我們選擇知數(shù)據(jù)中的兩點(diǎn)〔2,6.42〕和〔16,10.76〕代入方程，得到方程組：可決系數(shù)計(jì)算：上述方程組有兩種解法：手工，Matlab,下面引見(jiàn)Matlab解方程組的方法[x,y]=solve('6.42*(2*a+b)=2','10.76*(16*a+b)=16')取點(diǎn)：(2,6.42),(8,9.93),(10,10.49)代入上述方程[a,b,c]=solve('log(b)+c*2=log(6.42/a-1)','log(b)+c*10=log(10.49/a-1)','log(b)+c*8=log(9.93/a-1)')留意：假設(shè)出現(xiàn)復(fù)數(shù)解，那么只取實(shí)部作業(yè)給藥方案。顯然,要設(shè)計(jì)給藥方案,必需知道給藥后血藥濃度隨時(shí)間變化的規(guī)律.為此,從實(shí)驗(yàn)和實(shí)際兩方面著手.在實(shí)驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時(shí)辰t(小時(shí))采集血樣,測(cè)得血藥濃度c如表:血藥濃度c(t)的測(cè)試數(shù)據(jù)t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01例5給藥方案。近似直線(xiàn)關(guān)系,即c(t)有按負(fù)指數(shù)規(guī)律減少的趨勢(shì).例5給藥方案1.確定血藥濃度的變化規(guī)律假設(shè):a)藥物向體外排除的速率與中心室的血藥濃度成正比,比例系數(shù)為k(>0),稱(chēng)排除速率.b)中心室血液容積為常數(shù)V,t=0瞬時(shí)注入藥物的劑量為d,血藥濃度立刻為由假設(shè)a),中心室的血藥濃度c(t)應(yīng)滿(mǎn)足微分方程由假設(shè)b),方程的初始條件為:求解得:即血藥濃度c(t)按指數(shù)規(guī)律下降.2.給藥方案設(shè)計(jì)簡(jiǎn)單適用的給藥方案是:每隔一定時(shí)間,反復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并堅(jiān)持在c1-c2之間.x0yc1c22.給藥方案設(shè)計(jì)簡(jiǎn)單適用的給藥方案是:每隔一定時(shí)間,反復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并堅(jiān)持在c1-c2之間.為此,初次劑量需加大到D0.由式得到:顯然,當(dāng)c1,c2給定后,要確定給藥方案必需知道參數(shù)V和k.2.由實(shí)驗(yàn)數(shù)據(jù)作曲線(xiàn)擬合以確定參數(shù)問(wèn)題化為由數(shù)據(jù)ti,yi(i=1,…,8)擬合直線(xiàn)記為了用線(xiàn)性最小二乘法擬合的系數(shù)V和k,先取對(duì)數(shù)得用Matlab作線(xiàn)性最小二乘法擬合,得到問(wèn)題化為由數(shù)據(jù)ti,yi(i=1,…,8)擬合直線(xiàn)記為了用線(xiàn)性最小二乘法擬合的系數(shù)V和k,先取對(duì)數(shù)得用Matlab作線(xiàn)性最小二乘法擬合,得到由實(shí)驗(yàn)數(shù)據(jù)d=300(mg)算出:擬合曲線(xiàn)為:3.結(jié)論將k,V和給出的c1=10,c2=25代入得:D0=375.5,D=225.3,給藥方案無(wú)妨定為:D0=375mg,D=225mg,小時(shí).范例：薄膜浸透率的測(cè)定一、問(wèn)題：某種醫(yī)用薄膜，具有從高濃度的溶液向低濃度的溶液分散的功能，在試制時(shí)需測(cè)定薄膜被物質(zhì)分子穿透的才干。測(cè)定方法：用面積為S的薄膜將容器分成體積分別為的兩部份，在兩部分中分別注滿(mǎn)該物質(zhì)的兩種不同濃度的溶液。此時(shí)該物質(zhì)分子就會(huì)從高濃度溶液穿過(guò)薄膜向低濃度溶液中分散。平均每單位時(shí)間經(jīng)過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K表征了薄膜被該物質(zhì)分子穿透的才干，稱(chēng)為浸透率。定時(shí)丈量容器中薄膜某一側(cè)的溶液濃度，以此確定K。VAVBS二、問(wèn)題分析調(diào)查時(shí)段[t，t+Δt]薄膜兩側(cè)容器中該物質(zhì)質(zhì)量的變化。設(shè)，對(duì)容器的B部分溶液濃度的測(cè)試結(jié)果如下表：〔濃度單位〕1〕在容器的一側(cè)，物質(zhì)質(zhì)量的添加是由于另一側(cè)的物質(zhì)向該側(cè)浸透的結(jié)果，因此物質(zhì)質(zhì)量的增量應(yīng)等于另一側(cè)的該物質(zhì)向這側(cè)的浸透量。以容器A側(cè)為例，在時(shí)段[t，t+Δt]物質(zhì)質(zhì)量的增量為：分別表示在時(shí)辰t膜兩側(cè)溶液設(shè)的濃度，濃度單位:由于平均每單位時(shí)間經(jīng)過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)為K。因此，在時(shí)段[t，t+Δt]，從B側(cè)浸透至A側(cè)的該物質(zhì)的質(zhì)量為：于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：分別表示在初始時(shí)辰兩側(cè)溶液的濃度其中〔1〕2)留意到整個(gè)容器的溶液中含有該物質(zhì)的質(zhì)量不變,與初始時(shí)辰該物質(zhì)的含量一樣，因此從而：加上初值條件：代入式〔1〕得：便可得出CB(t)的變化規(guī)律，從而根據(jù)實(shí)驗(yàn)數(shù)據(jù)進(jìn)展擬合，估計(jì)出參數(shù)K,。三、數(shù)學(xué)模型假設(shè)：1〕薄膜兩側(cè)的溶液一直是均勻的；2〕平均每單位時(shí)間經(jīng)過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比。3〕薄膜是雙向同性的即物質(zhì)從膜的任何一側(cè)向另一側(cè)浸透的性能是一樣的。基于假設(shè)和前面的分析，B側(cè)的濃度CB(t)應(yīng)滿(mǎn)足如下微分方程和初始條件：四、求解方法：1.函數(shù)擬合法前面得到的模型是一個(gè)帶初值的一階線(xiàn)性微分方程，解之得：?jiǎn)栴}歸結(jié)為利用CB在時(shí)辰tj的丈量數(shù)據(jù)Cj(j=1,2,...,N)來(lái)辨識(shí)K和。引入從而用函數(shù)CB(t)來(lái)擬合所給的實(shí)驗(yàn)數(shù)據(jù)，從而估計(jì)出其中的參數(shù)a,b,K。將代入上式有：用MATLAB軟件進(jìn)展計(jì)算.1〕編寫(xiě)函數(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)在任務(wù)空間中執(zhí)行以下命令(test1.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)3)輸出結(jié)果:x=0.007-0.0030.1012即k=0.1012,a=0.007,b=-0.003,nongdutest1進(jìn)一步求得：2.非線(xiàn)性規(guī)劃法利用CB在時(shí)辰tj的丈量數(shù)據(jù)Cj(j=1,2,...,N)來(lái)辨識(shí)K和。問(wèn)題可轉(zhuǎn)化為求函數(shù)即求函數(shù)的最小值點(diǎn)〔K，a，b〕。3.導(dǎo)函數(shù)擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1,2,...,N)假設(shè)知道一組數(shù)據(jù)那么可用最小二乘擬合的方法來(lái)求出函數(shù)中的未知參數(shù)K和h。即為求參數(shù)K,a使以下誤差函數(shù)到達(dá)最?。涸搯?wèn)題等價(jià)于用函數(shù)f(K,a,CB)=K(0.01a-0.02CB)來(lái)擬合數(shù)據(jù)(j=1,2,...,N)用MATLAB軟件進(jìn)展計(jì)算.％求數(shù)據(jù)點(diǎn)(j=1,2,...,N)tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...610626639650659];[d,ifail]=e01bef(tdata,cdata);[cj,dcj]=e01bgf(tdata,cdata,d,tdata);1〕編寫(xiě)函數(shù)M-文件baomof.mfunctionf=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中x(1)=K;x(2)=h2)編寫(xiě)命令M文件(baomo21.m)3)輸出結(jié)果:x=0.10090.014即k=0.1009,h=0.014％作函數(shù)擬合x(chóng)0=[0.2,0.1];x=lsqcurvefit('baomof',x0,cdata,dcj')4.線(xiàn)性化迭代法前面帶初始條件的一階線(xiàn)性微分方程的解為其中：假設(shè)得到了參數(shù)K的一個(gè)較好的近似值K*，那么將CB(t)關(guān)于K在K*處展開(kāi)，略去K的二次及以上的項(xiàng)得CB(t)的一個(gè)近似式經(jīng)過(guò)極小化確定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一個(gè)新的近似值，用同樣的方法再求新的修正值K。這個(gè)過(guò)程可以不斷反復(fù)，直到修正值足夠小為止