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

數(shù)據(jù)擬合與線性最小二乘法數(shù)據(jù)擬合與線性最小二乘法教學(xué)內(nèi)容數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題的求解思路線性最小二乘法原理算例MATLAB工具箱Curvefit演示教學(xué)內(nèi)容數(shù)據(jù)擬合問(wèn)題的提法

設(shè)

R=at+ba,b為待定系數(shù)求電阻R隨溫度t的變化規(guī)律。已知熱敏電阻數(shù)據(jù)：溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032引例1：熱敏電阻電阻值的變化規(guī)律設(shè)R=at+b求電阻R隨溫度t的變化規(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ī)律t(h)0.250.5數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題：已知一維（二維，…）數(shù)據(jù)，即平面上的n個(gè)點(diǎn)(xi，yi)，i=1,2,…,n，xi互不相同，尋求一個(gè)函數(shù)(曲線)y=f(x)，使f(x)在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合的最好，如下圖所示(圖中δi為(xi，yi)與y=f(x)的距離)。

Oxyδi(xi，yi)數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題：已知一維（二維，…）數(shù)據(jù)，即數(shù)據(jù)擬合問(wè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

的平方和最小，稱為最小二乘準(zhǔn)則。

擬合準(zhǔn)則還有如最小一乘準(zhǔn)則、極大極小準(zhǔn)則等。數(shù)據(jù)擬合問(wèn)題的求解思路線性最小二乘法是解決數(shù)據(jù)擬合最常用的方線性最小二乘法原理1.理論———基本理論之a(chǎn)k的確定

根據(jù)最小二乘準(zhǔn)則，記J(a1，a2，…，am)=為求a1，a2，…，am是J達(dá)到最小，只需要利用極值的必要條件

得到關(guān)于a1，…，am的線性方程組

線性最小二乘法原理1.理論———基本理論之a(chǎn)k的確定根線性最小二乘法原理記

，A=(a1，a2，…，am)T，y=(y1，…，yn)T，

方程組(3)可表為RTRA=RTy(4)(4)稱為法方程組，當(dāng){r1(x)，…，rm(x)}線性無(wú)關(guān)時(shí)，R列滿秩，RTR可逆，于是方程組(4)有唯一解A=(RTR)-1RTy(5)可以看出，只要f(x)關(guān)于待定系數(shù)a1，…，am線性，在最小二乘準(zhǔn)則（2）下得到的方程組(3)關(guān)于a1，a2，…，am也一定是線性的，故稱線性最小二乘法。

線性最小二乘法原理記，A=(a1，a2，…，am)T，y=線性最小二乘法原理2.理論______函數(shù)rk(x)的選取

對(duì)數(shù)據(jù)(xi，yi)用線性最小二乘法作擬合時(shí)，首要的、也是關(guān)鍵的一步是恰當(dāng)?shù)剡x取r1(x)，r2(x)，…，rm(x)。n

如果通過(guò)機(jī)理分析，能夠知道y與x之間應(yīng)該有什么樣的函數(shù)關(guān)系，則r1(x)，…，rm(x)容易確定。n

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

直線y=a1x+a2u

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

雙曲線（一支）y=a1/x+a2u

指數(shù)曲線：擬合前需作變量代換，化為線性函數(shù)。對(duì)已知數(shù)據(jù)，用什么樣的曲線擬合最好，可以在直觀判斷的基礎(chǔ)上，選擇幾種曲線分別作擬合，然后比較，看那條曲線的最小二乘指標(biāo)J最小。

線性最小二乘法原理2.理論______函數(shù)rk(x)的選取線性最小二乘法原理3.求解方法l

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

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

比較多種擬合結(jié)果，選擇其中較好的一種或者某幾種作為備選結(jié)果；

注·意：通常需要將非線性函數(shù)rk(x)的轉(zhuǎ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。

線性最小二乘法原理3.求解方法l

描出數(shù)據(jù)的圖示；注·意算例【例】給定數(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.135【解】：（1）、先描出數(shù)據(jù)的圖示

算例【例】給定數(shù)據(jù)(xi，f(xi))，i=0,1,2,3,算例（2）選定不同的數(shù)學(xué)函數(shù)（模型）或者rk(x)進(jìn)行擬合

l

直線模型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算例（2）選定不同的數(shù)學(xué)函數(shù)（模型）或者rk(x)進(jìn)行擬合求解法方程組

得到a=1.6380，b=3.3520，于是得到該模型下的最小二乘擬合曲線為g(x)=1.6380+3.3520x。

算例求解法方程組得到a=1.6380，b=3.3520，算例算例

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算例

l

多項(xiàng)式模型y=a0+a1x+a2x2i算例求解法方程組得到

a=3.6294，b=0.5406，c=0.9371，于是得到該模型下的最小二乘擬合曲線為g(x)=3.6294+0.5406x+0.9371x2。

算例求解法方程組得到算例

l

雙曲線模型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.1531390.1342280.118203r1(x)11111r2(x)1.001.251.501.752.00算例

l

雙曲線模型y=1/(a0+a1x)i0算例求解法方程組得到

a0=0.27139，a1=-0.07768，于是得到該模型下的最小二乘擬合曲線為g(x)=1/(0.27139-0.07768x)。

算例求解法方程組得到算例

l

指數(shù)曲線模型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.13534917r1(x)11111r2(x)1.001.251.501.752.00算例

l

指數(shù)曲線模型y=aebxi01234x算例求解法方程組得到

A=1.1225，b=0.5057，a=eA=3.072525924，于是得到此模型下的最小二乘擬合曲線

g(x)=3.072525924e0.5057x

。

算例求解法方程組得到算例(3)比較上面的結(jié)果如下:i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46直線模型4.995.836.677.508.340.49e-1多項(xiàng)式模型5.1075.7696.5497.4458.460.85e-3雙曲線模型5.1625.7386.4577.3838.6180.42e-1指數(shù)函數(shù)模型5.0955.7816.5607.4448.4480.12e-2經(jīng)過(guò)比較，使用多項(xiàng)式與指數(shù)函數(shù)擬合效果較好。

算例(3)比較上面的結(jié)果如下:i01234xi1.001用Matlab進(jìn)行數(shù)據(jù)擬合1.多項(xiàng)式曲線擬合: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.用Matlab進(jìn)行數(shù)據(jù)擬合1.多項(xiàng)式曲線擬合:poly例1已知觀測(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)式曲線擬合這些數(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程序如下:例1已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy0-0.4470.11.9t=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)t=0:0.1:1.2x=0:0.1:1例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(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例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,解:描出散點(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]解:描出散點(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解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]a例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(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擬合曲線為:y=-0.3012t+29.3804例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,例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ù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.2.89例3一個(gè)15.4cm×30.48cm的混凝土柱在加壓實(shí)驗(yàn)已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.解選取指數(shù)函數(shù)作擬合時(shí),在擬合前需作變量代換,化為k1,k2

的線性函數(shù).于是,令即已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(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-')在命令窗口輸入:x=[500*1.0e-61000*1.0已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.擬合曲線為:令則求得于是已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線有:直線多項(xiàng)式一般n=2,3,不宜過(guò)高.雙曲線(一支)指數(shù)曲線在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線有:直線多項(xiàng)式一般n=2,3,2.非線性曲線擬合: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.2.非線性曲線擬合:lsqcurvefit.功能:x=例4已知觀測(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的值,使得曲線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已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.例4已知觀測(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的值,使得曲線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)例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.編寫(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)行后顯示x=3.00224.03040.9404resnorm=0.0912編寫(xiě)下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;程序例4已知觀測(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的值,使得曲線f(x)=aex+bx2+cx3

與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.說(shuō)明:最小二乘意義上的最佳擬合函數(shù)為f(x)=3ex+4.03x2+0.94x3.此時(shí)的殘差是:0.0912.例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.f(x)=3ex+4.03x2+0.94x3.擬合函數(shù)為:f(x)=3ex+4.03x2+0.94x3.擬合作業(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.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]用命令lsqcurvefit(‘f’,a0,x,y)

作業(yè)：用函數(shù)f(x)=a1*exp(-a2*x)+a3*ex作業(yè)：考察某種纖維的強(qiáng)度與其拉伸倍數(shù)的關(guān)系,下表是實(shí)際測(cè)定的24個(gè)纖維樣品的強(qiáng)度與相應(yīng)的拉伸倍數(shù)的記錄:試建立纖維強(qiáng)度與拉伸倍數(shù)之間的關(guān)系。作業(yè)：考察某種纖維的強(qiáng)度與其拉伸倍數(shù)的關(guān)系,下表是實(shí)際測(cè)定的作業(yè)：煉鋼廠出鋼時(shí)所用盛鋼水的鋼包，由于鋼水對(duì)耐火材料的侵蝕，容積不斷增大，我們希望找出使用次數(shù)與增大容積之間的函數(shù)關(guān)系.實(shí)驗(yàn)數(shù)據(jù)如下：表4.2鋼包使用次數(shù)與增大容積使用次數(shù)23456789增大容積6.428.29.589.59.7109.939.99使用次數(shù)10111213141516增大容積10.4910.5910.610.810.610.910.76分別選擇函數(shù)擬合鋼包容積與使用次數(shù)的關(guān)系,在同一坐標(biāo)系內(nèi)作出函數(shù)圖形.作業(yè)：煉鋼廠出鋼時(shí)所用盛鋼水的鋼包，由于鋼水對(duì)耐火材料的侵蝕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ě)出具體函數(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ì)算：x1=[2:16];下面給出分式函數(shù)擬合程序：初始參數(shù)b0的上述方程組有兩種解法：手工，Matlab,下面介紹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)')注意：如果出現(xiàn)復(fù)數(shù)解，則只取實(shí)部上述方程組有兩種解法：手工，Matlab,下面介紹Matla作業(yè)給藥方案。

顯然,要設(shè)計(jì)給藥方案,必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律.為此,從實(shí)驗(yàn)和理論兩方面著手.在實(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作業(yè)給藥方案。顯然,要設(shè)計(jì)給藥方案,必須知道給例5給藥方案。

近似直線關(guān)系,即c(t)有按負(fù)指數(shù)規(guī)律減少的趨勢(shì).例5給藥方案。近似直線關(guān)系,即c(t)有按負(fù)指數(shù)例5給藥方案1.確定血藥濃度的變化規(guī)律假設(shè):a)藥物向體外排除的速率與中心室的血藥濃度成正比,比例系數(shù)為k(>0),稱排除速率.b)中心室血液容積為常數(shù)V,t=0瞬時(shí)注入藥物的劑量為d,血藥濃度立即為由假設(shè)a),中心室的血藥濃度c(t)應(yīng)滿足微分方程由假設(shè)b),方程的初始條件為:求解得:即血藥濃度c(t)按指數(shù)規(guī)律下降.例5給藥方案1.確定血藥濃度的變化規(guī)律假設(shè):a)藥物2.給藥方案設(shè)計(jì)簡(jiǎn)單實(shí)用的給藥方案是:每隔一定時(shí)間,重復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并保持在c1-c2

之間.x0yc1c22.給藥方案設(shè)計(jì)簡(jiǎn)單實(shí)用的給藥方案是:每隔一定時(shí)間2.給藥方案設(shè)計(jì)簡(jiǎn)單實(shí)用的給藥方案是:每隔一定時(shí)間,重復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并保持在c1-c2

之間.為此,初次劑量需加大到D0.由式得到:顯然,當(dāng)c1,c2

給定后,要確定給藥方案必須知道參數(shù)V

和k.2.給藥方案設(shè)計(jì)簡(jiǎn)單實(shí)用的給藥方案是:每隔一定時(shí)間2.由實(shí)驗(yàn)數(shù)據(jù)作曲線擬合以確定參數(shù)問(wèn)題化為由數(shù)據(jù)ti,yi

(i=1,…,8)擬合直線記為了用線性最小二乘法擬合的系數(shù)V

和k,先取對(duì)數(shù)得用Matlab作線性最小二乘法擬合,得到2.由實(shí)驗(yàn)數(shù)據(jù)作曲線擬合以確定參數(shù)問(wèn)題化為由數(shù)據(jù)ti,問(wèn)題化為由數(shù)據(jù)ti,yi

(i=1,…,8)擬合直線記為了用線性最小二乘法擬合的系數(shù)V

和k,先取對(duì)數(shù)得用Matlab作線性最小二乘法擬合,得到由實(shí)驗(yàn)數(shù)據(jù)d=300(mg)算出:問(wèn)題化為由數(shù)據(jù)ti,yi(i=1,…,8)擬擬合曲線為:擬合曲線為:3.結(jié)論將k,V

和給出的c1=10,c2=25代入得:D0=375.5,D=225.3,給藥方案不妨定為:D0=375mg,D=225mg,小時(shí).3.結(jié)論將k,V和給出的c1=10,c2=25范例：薄膜滲透率的測(cè)定

一、問(wèn)題：某種醫(yī)用薄膜，具有從高濃度的溶液向低濃度的溶液擴(kuò)散的功能，在試制時(shí)需測(cè)定薄膜被物質(zhì)分子穿透的能力。測(cè)定方法：用面積為S的薄膜將容器分成體積分別為的兩部份，在兩部分中分別注滿該物質(zhì)的兩種不同濃度的溶液。此時(shí)該物質(zhì)分子就會(huì)從高濃度溶液穿過(guò)薄膜向低濃度溶液中擴(kuò)散。平均每單位時(shí)間通過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K表征了薄膜被該物質(zhì)分子穿透的能力，稱為滲透率。定時(shí)測(cè)量容器中薄膜某一側(cè)的溶液濃度，以此確定K。VAVBS范例：薄膜滲透率的測(cè)定一、問(wèn)題：VAVBS二、問(wèn)題分析考察時(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è)的滲透量。二、問(wèn)題分析設(shè)，對(duì)

以容器A側(cè)為例，在時(shí)段[t，t+Δt]物質(zhì)質(zhì)量的增量為：分別表示在時(shí)刻t膜兩側(cè)溶液設(shè)的濃度，濃度單位:

由于平均每單位時(shí)間通過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)為K。

因此，在時(shí)段[t，t+Δt]，從B側(cè)滲透至A側(cè)的該物質(zhì)的質(zhì)量為：以容器A側(cè)為例，在時(shí)段[t，t+Δt]物質(zhì)質(zhì)量的增量為：于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：分別表示在初始時(shí)刻兩側(cè)溶液的濃度其中（1）2)注意到整個(gè)容器的溶液中含有該物質(zhì)的質(zhì)量不變,與初始時(shí)刻該物質(zhì)的含量相同，因此于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：分別從而：加上初值條件：代入式（1）得：便可得出CB(t)的變化規(guī)律，從而根據(jù)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合，估計(jì)出參數(shù)K,。從而：加上初值條件：代入式（1）得：便可得出CB(t)的變化三、數(shù)學(xué)模型假設(shè)：1）薄膜兩側(cè)的溶液始終是均勻的；2）平均每單位時(shí)間通過(guò)單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比。3）薄膜是雙向同性的即物質(zhì)從膜的任何一側(cè)向另一側(cè)滲透的性能是相同的?；诩僭O(shè)和前面的分析，B側(cè)的濃度CB(t)應(yīng)滿足如下微分方程和初始條件：三、數(shù)學(xué)模型四、求解方法：1.函數(shù)擬合法前面得到的模型是一個(gè)帶初值的一階線性微分方程，解之得：?jiǎn)栴}歸結(jié)為利用CB在時(shí)刻tj的測(cè)量數(shù)據(jù)Cj(j=1,2,...,N)來(lái)辨識(shí)K和。四、求解方法：1.函數(shù)擬合法前面得到的模型是一個(gè)帶初值的一引入從而

用函數(shù)CB(t)來(lái)擬合所給的實(shí)驗(yàn)數(shù)據(jù)，從而估計(jì)出其中的參數(shù)a,b,K。將代入上式有：引入從而用函數(shù)CB(t)來(lái)擬合所給的實(shí)驗(yàn)數(shù)據(jù)用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)在工作空間中執(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用MATLAB軟件進(jìn)行計(jì)算.nongdutest1進(jìn)一步求得：2.非線性規(guī)劃法

利用CB在時(shí)刻tj的測(cè)量數(shù)據(jù)Cj(j=1,2,...,N)來(lái)辨識(shí)K和。問(wèn)題可轉(zhuǎn)化為求函數(shù)即求函數(shù)的最小值點(diǎn)（K，a，b）。進(jìn)一步求得：2.非線性規(guī)劃法利用CB在時(shí)刻tj的測(cè)3.導(dǎo)函數(shù)擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1,2,...,N)若知道一組數(shù)據(jù)則可用最小二乘擬合的方法來(lái)求出函數(shù)中的未知參數(shù)K和h。3.導(dǎo)函數(shù)擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看即為求參數(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)即為求參數(shù)K,a使下列誤差函數(shù)達(dá)到最?。涸搯?wè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)用MATLAB軟件進(jìn)行計(jì)算.％求數(shù)據(jù)點(diǎn)(j=1,2,...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')3)輸出結(jié)果:x=0.10090.014％作函數(shù)4.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中：

如果得到了參數(shù)K的一個(gè)較好的近似值K*，則將CB(t)關(guān)于K在K*處展開(kāi)，略去K的二次及以上的項(xiàng)得CB(t)的一個(gè)近似式通過(guò)極小化4.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中確定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一個(gè)新的近似值，用同樣的方法再求新的修正值K。這個(gè)過(guò)程可以不斷重復(fù)，直到修正值足夠小為止。1）當(dāng)K的初值取為k=0.3時(shí)，出現(xiàn)奇異情況，迭代不收斂；2）當(dāng)K的初值取為k=0.2時(shí)，經(jīng)四次迭代，已經(jīng)收斂到一個(gè)很好的解。迭代結(jié)果如下表。確定a,b,d,再由K=d/0.02b得到K*的修正五、結(jié)果及誤差分析

幾種方法得出的結(jié)果及相應(yīng)的誤差總結(jié)于下表，誤差為計(jì)算數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)之差的平方和。注：導(dǎo)函數(shù)擬合法得出的參數(shù)值精度有限，線性化迭代法要求參數(shù)的初值比較接近精確值。因此可將導(dǎo)函數(shù)擬合法和線性化迭代法結(jié)合起來(lái)使用，把前者得到的參數(shù)K的值作為迭代法中K的初值，這樣可使迭代法收斂或收斂更快。3）取K的初值為k=0.1009,只一次迭代就得到2）中的最后結(jié)果。五、結(jié)果及誤差分析幾種方法得出的結(jié)果及相應(yīng)的數(shù)模培訓(xùn)_數(shù)據(jù)擬合方法函數(shù)擬合法的擬合效果函數(shù)擬合法的擬合效果求解參數(shù)辨識(shí)模型的方法：函數(shù)擬合；非線性規(guī)劃；導(dǎo)函數(shù)擬合；線性化迭代；其它方法。求解參數(shù)辨識(shí)模型的方法：布置“函數(shù)擬合”實(shí)驗(yàn)?zāi)康?.掌握用MATLAB計(jì)算函數(shù)擬合的方法內(nèi)容2.用函數(shù)擬合方法解決實(shí)際問(wèn)題。?給藥方案一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案。在快速靜脈注射下，所謂給藥方案是指，每次注射計(jì)量多大，間隔時(shí)間多長(zhǎng)。藥物進(jìn)入肌體后隨血液輸送到全身，在這過(guò)程中不斷被吸收、分解、代謝，最終排出體外。藥物向體外排出的速率與血藥濃度成正比。單位體積血液中的藥物含量，稱血藥濃度。臨床上，每種藥物有一個(gè)最小有效濃度c1和最大治療濃度c2。設(shè)計(jì)給藥方案時(shí)，要使血藥濃度保持在c1——c2之間，設(shè)本題研究的藥物的c1=10(g/ml),c2=25(g/ml),對(duì)某人用快速靜脈注射方式一次注入該藥物300mg后，在一定時(shí)刻t（小時(shí)）采集血樣，測(cè)得血藥濃度c(g/ml),如下頁(yè)表。試設(shè)計(jì)該藥的給藥方案。布置“函數(shù)擬合”實(shí)驗(yàn)?zāi)康?.掌握用MATLAB計(jì)算函實(shí)驗(yàn)報(bào)告要求1.目的。2.內(nèi)容。3.模型；算法；計(jì)算結(jié)果；結(jié)果分析；附程序。4.心得體會(huì)和建議。t0.250.511.523468c19.2118.1515.3614.112.899.327.455.243.01實(shí)驗(yàn)報(bào)告要求1.目的。2.內(nèi)容。3.模型；

一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案。在快速靜脈注射下，所謂給藥方案是指，每次注射計(jì)量多大，間隔時(shí)間多長(zhǎng)。藥物進(jìn)入肌體后隨血液輸送到全身，在這過(guò)程中不斷被吸收、分解、代謝，最終排出體外。藥物向體外排出的速率與血藥濃度成正比。單位體積血液中的藥物含量，稱血藥濃度。臨床上，每種藥物有一個(gè)最小有效濃度c1和最大治療濃度c2。設(shè)計(jì)給藥方案時(shí)，要使血藥濃度保持在c1——c2之間，設(shè)本題研究的藥物的c1=10(g/ml),c2=25(g/ml),對(duì)某人用快速靜脈注射方式一次注入該藥物300mg后，在一定時(shí)刻t（小時(shí)）采集血樣，測(cè)得血藥濃度c(g/ml),如下表：模擬訓(xùn)練題：給藥方案一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案。在快t0.250.511.523468c19.2118.1515.3614.112.899.327.455.243.01試設(shè)計(jì)該藥的給藥方案。t0.250.511.演講完畢，謝謝觀看！演講完畢，謝謝觀看！數(shù)據(jù)擬合與線性最小二乘法數(shù)據(jù)擬合與線性最小二乘法教學(xué)內(nèi)容數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題的求解思路線性最小二乘法原理算例MATLAB工具箱Curvefit演示教學(xué)內(nèi)容數(shù)據(jù)擬合問(wèn)題的提法

設(shè)

R=at+ba,b為待定系數(shù)求電阻R隨溫度t的變化規(guī)律。已知熱敏電阻數(shù)據(jù)：溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032引例1：熱敏電阻電阻值的變化規(guī)律設(shè)R=at+b求電阻R隨溫度t的變化規(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ī)律t(h)0.250.5數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題：已知一維（二維，…）數(shù)據(jù)，即平面上的n個(gè)點(diǎn)(xi，yi)，i=1,2,…,n，xi互不相同，尋求一個(gè)函數(shù)(曲線)y=f(x)，使f(x)在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合的最好，如下圖所示(圖中δi為(xi，yi)與y=f(x)的距離)。

Oxyδi(xi，yi)數(shù)據(jù)擬合問(wèn)題的提法數(shù)據(jù)擬合問(wèn)題：已知一維（二維，…）數(shù)據(jù)，即數(shù)據(jù)擬合問(wè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

的平方和最小，稱為最小二乘準(zhǔn)則。

擬合準(zhǔn)則還有如最小一乘準(zhǔn)則、極大極小準(zhǔn)則等。數(shù)據(jù)擬合問(wèn)題的求解思路線性最小二乘法是解決數(shù)據(jù)擬合最常用的方線性最小二乘法原理1.理論———基本理論之a(chǎn)k的確定

根據(jù)最小二乘準(zhǔn)則，記J(a1，a2，…，am)=為求a1，a2，…，am是J達(dá)到最小，只需要利用極值的必要條件

得到關(guān)于a1，…，am的線性方程組

線性最小二乘法原理1.理論———基本理論之a(chǎn)k的確定根線性最小二乘法原理記

，A=(a1，a2，…，am)T，y=(y1，…，yn)T，

方程組(3)可表為RTRA=RTy(4)(4)稱為法方程組，當(dāng){r1(x)，…，rm(x)}線性無(wú)關(guān)時(shí)，R列滿秩，RTR可逆，于是方程組(4)有唯一解A=(RTR)-1RTy(5)可以看出，只要f(x)關(guān)于待定系數(shù)a1，…，am線性，在最小二乘準(zhǔn)則（2）下得到的方程組(3)關(guān)于a1，a2，…，am也一定是線性的，故稱線性最小二乘法。

線性最小二乘法原理記，A=(a1，a2，…，am)T，y=線性最小二乘法原理2.理論______函數(shù)rk(x)的選取

對(duì)數(shù)據(jù)(xi，yi)用線性最小二乘法作擬合時(shí)，首要的、也是關(guān)鍵的一步是恰當(dāng)?shù)剡x取r1(x)，r2(x)，…，rm(x)。n

如果通過(guò)機(jī)理分析，能夠知道y與x之間應(yīng)該有什么樣的函數(shù)關(guān)系，則r1(x)，…，rm(x)容易確定。n

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

直線y=a1x+a2u

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

雙曲線（一支）y=a1/x+a2u

指數(shù)曲線：擬合前需作變量代換，化為線性函數(shù)。對(duì)已知數(shù)據(jù)，用什么樣的曲線擬合最好，可以在直觀判斷的基礎(chǔ)上，選擇幾種曲線分別作擬合，然后比較，看那條曲線的最小二乘指標(biāo)J最小。

線性最小二乘法原理2.理論______函數(shù)rk(x)的選取線性最小二乘法原理3.求解方法l

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

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

比較多種擬合結(jié)果，選擇其中較好的一種或者某幾種作為備選結(jié)果；

注·意：通常需要將非線性函數(shù)rk(x)的轉(zhuǎ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。

線性最小二乘法原理3.求解方法l

描出數(shù)據(jù)的圖示；注·意算例【例】給定數(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.135【解】：（1）、先描出數(shù)據(jù)的圖示

算例【例】給定數(shù)據(jù)(xi，f(xi))，i=0,1,2,3,算例（2）選定不同的數(shù)學(xué)函數(shù)（模型）或者rk(x)進(jìn)行擬合

l

直線模型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算例（2）選定不同的數(shù)學(xué)函數(shù)（模型）或者rk(x)進(jìn)行擬合求解法方程組

得到a=1.6380，b=3.3520，于是得到該模型下的最小二乘擬合曲線為g(x)=1.6380+3.3520x。

算例求解法方程組得到a=1.6380，b=3.3520，算例算例

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算例

l

多項(xiàng)式模型y=a0+a1x+a2x2i算例求解法方程組得到

a=3.6294，b=0.5406，c=0.9371，于是得到該模型下的最小二乘擬合曲線為g(x)=3.6294+0.5406x+0.9371x2。

算例求解法方程組得到算例

l

雙曲線模型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.1531390.1342280.118203r1(x)11111r2(x)1.001.251.501.752.00算例

l

雙曲線模型y=1/(a0+a1x)i0算例求解法方程組得到

a0=0.27139，a1=-0.07768，于是得到該模型下的最小二乘擬合曲線為g(x)=1/(0.27139-0.07768x)。

算例求解法方程組得到算例

l

指數(shù)曲線模型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.13534917r1(x)11111r2(x)1.001.251.501.752.00算例

l

指數(shù)曲線模型y=aebxi01234x算例求解法方程組得到

A=1.1225，b=0.5057，a=eA=3.072525924，于是得到此模型下的最小二乘擬合曲線

g(x)=3.072525924e0.5057x

。

算例求解法方程組得到算例(3)比較上面的結(jié)果如下:i01234xi1.001.251.501.752.00f(xi)5.105.796.537.458.46直線模型4.995.836.677.508.340.49e-1多項(xiàng)式模型5.1075.7696.5497.4458.460.85e-3雙曲線模型5.1625.7386.4577.3838.6180.42e-1指數(shù)函數(shù)模型5.0955.7816.5607.4448.4480.12e-2經(jīng)過(guò)比較，使用多項(xiàng)式與指數(shù)函數(shù)擬合效果較好。

算例(3)比較上面的結(jié)果如下:i01234xi1.001用Matlab進(jìn)行數(shù)據(jù)擬合1.多項(xiàng)式曲線擬合: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.用Matlab進(jìn)行數(shù)據(jù)擬合1.多項(xiàng)式曲線擬合:poly例1已知觀測(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)式曲線擬合這些數(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程序如下:例1已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy0-0.4470.11.9t=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)t=0:0.1:1.2x=0:0.1:1例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(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例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,解:描出散點(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]解:描出散點(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解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]a例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(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擬合曲線為:y=-0.3012t+29.3804例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,例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ù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.2.89例3一個(gè)15.4cm×30.48cm的混凝土柱在加壓實(shí)驗(yàn)已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.解選取指數(shù)函數(shù)作擬合時(shí),在擬合前需作變量代換,化為k1,k2

的線性函數(shù).于是,令即已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(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-')在命令窗口輸入:x=[500*1.0e-61000*1.0已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.擬合曲線為:令則求得于是已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來(lái)描述,即假設(shè)式中,在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線有:直線多項(xiàng)式一般n=2,3,不宜過(guò)高.雙曲線(一支)指數(shù)曲線在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線有:直線多項(xiàng)式一般n=2,3,2.非線性曲線擬合: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.2.非線性曲線擬合:lsqcurvefit.功能:x=例4已知觀測(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的值,使得曲線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已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.例4已知觀測(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的值,使得曲線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)例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.編寫(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)行后顯示x=3.00224.03040.9404resnorm=0.0912編寫(xiě)下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;程序例4已知觀測(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的值,使得曲線f(x)=aex+bx2+cx3

與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.說(shuō)明