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

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

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

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

直線y=a1x+a2u

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

雙曲線〔一支〕y=a1/x+a2u

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

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

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

比較多種擬合結(jié)果，選擇其中較好的一種或者某幾種作為備選結(jié)果；注·意：通常需求將非線性函數(shù)rk(x)的轉(zhuǎn)化成線性的函數(shù)Rk(x)，然后再用Rk(x)進展擬合，計算中通常需求列下表：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，見下表，使選擇適當(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)進展擬合

l

直線模型y=a+bx選取y=a+bx，此時，r1(x)=1，r2(x)=x。要求y=a+bx與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，yi=f(xi)。列表計算如下：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，于是得到該模型下的最小二乘擬合曲線為g(x)=1.6380+3.3520x。算例算例

l

多項式模型y=a0+a1x+a2x2選取Y=a+bx+cx2，此時，r1(x)=1，r2(x)=x，r3(x)=x2。要求y=a+bx+cx2與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，yi=f(xi)。列表計算如下：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，于是得到該模型下的最小二乘擬合曲線為g(x)=3.6294+0.5406x+0.9371x2。算例

l

雙曲線模型y=1/(a0+a1x)選取y=1/(a0+a1x)，令Y=1/y=a0+a1x，此時，r1(x)=1，r2(x)=x。要求Y=a0+a1x與(xi，yi)，i=0,1,2,3,4，做最小二乘擬合，Yi=1/f(xi)。列表計算如下：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，于是得到該模型下的最小二乘擬合曲線為g(x)=1/(0.27-0.07768x)。算例

l

指數(shù)曲線模型y=aebx選擇y=aebx，取對數(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)。計算結(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，于是得到此模型下的最小二乘擬合曲線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多項式模型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)過比較，運用多項式與指數(shù)函數(shù)擬合效果較好。用Matlab進展數(shù)據(jù)擬合1.多項式曲線擬合:polyfit.y0=polyval(p,x0)p=polyfit(x,y,m)其中,x,y為知數(shù)據(jù)點向量,分別表示橫,縱坐標,m為擬合多項式的次數(shù),結(jié)果前往m次擬合多項式系數(shù),從高次到低次存放在向量p中.可求得多項式在x0處的值y0.例1知觀測數(shù)據(jù)點如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分別用3次和6次多項式曲線擬合這些數(shù)據(jù)點.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)編寫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用切削機床進展金屬品加工時,為了適當(dāng)?shù)卣{(diào)整機床,需求測定刀具的磨損速度.在一定的時間丈量刀具的厚度,得數(shù)據(jù)如表所示:切削時間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm解:描出散點圖,在命令窗口輸入: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,'*')解:描出散點圖,在命令窗口輸入: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用切削機床進展金屬品加工時,為了適當(dāng)?shù)卣{(diào)整機床,需求測定刀具的磨損速度.在一定的時間丈量刀具的厚度,得數(shù)據(jù)如表所示:切削時間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm擬合曲線為:y=-0.3012t+29.3804例3一個15.4cm×30.48cm的混凝土柱在加壓實驗中的應(yīng)力-應(yīng)變關(guān)系測試點的數(shù)據(jù)如表所示1.552.472.933.03知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描畫,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.2.89知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描畫,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.解選取指數(shù)函數(shù)作擬合時,在擬合前需作變量代換,化為k1,k2的線性函數(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ù)曲線來描畫,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.擬合曲線為:令那么求得于是在實踐運用中常見的擬合曲線有:直線多項式普通n=2,3,不宜過高.雙曲線(一支)指數(shù)曲線2.非線性曲線擬合:lsqcurvefit.功能:x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)根據(jù)給定的數(shù)據(jù)xdata,ydata(對應(yīng)點的橫,縱坐標),按函數(shù)文件fun給定的函數(shù),以x0為初值作最小二乘擬合,前往函數(shù)fun中的系數(shù)向量x和殘差的平方和resnorm.例4知觀測數(shù)據(jù)點如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與知數(shù)據(jù)點在最小二乘意義上充分接近.首先編寫存儲擬合函數(shù)的函數(shù)文件.functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.^2+x(3)*xdata.^3保管為文件nihehanshu.m例4知觀測數(shù)據(jù)點如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與知數(shù)據(jù)點在最小二乘意義上充分接近.編寫下面的程序調(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)編寫下面的程序調(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)程序運轉(zhuǎn)后顯示x=3.00224.03040.9404resnorm=0.0912例4知觀測數(shù)據(jù)點如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與知數(shù)據(jù)點在最小二乘意義上充分接近.闡明:最小二乘意義上的最正確擬合函數(shù)為f(x)=3ex+4.03x2+0.94x3.此時的殘差是:0.0912.f(x)=3ex+4.03x2+0.94x3.擬合函數(shù)為:作業(yè)：用函數(shù)f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合以下數(shù)據(jù)點：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)查某種纖維的強度與其拉伸倍數(shù)的關(guān)系,下表是實踐測定的24個纖維樣品的強度與相應(yīng)的拉伸倍數(shù)的記錄:試建立纖維強度與拉伸倍數(shù)之間的關(guān)系。作業(yè)：煉鋼廠出鋼時所用盛鋼水的鋼包，由于鋼水對耐火資料的侵蝕，容積不斷增大，我們希望找出運用次數(shù)與增大容積之間的函數(shù)關(guā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)系,在同一坐標系內(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寫出詳細函數(shù)plot(x1,y1,'*',x1,y,'-or');下面給出分式函數(shù)擬合程序：初始參數(shù)b0的計算，由于確定兩個參數(shù)值，因此我們選擇知數(shù)據(jù)中的兩點〔2,6.42〕和〔16,10.76〕代入方程，得到方程組：可決系數(shù)計算：上述方程組有兩種解法：手工，Matlab,下面引見Matlab解方程組的方法[x,y]=solve('6.42*(2*a+b)=2','10.76*(16*a+b)=16')取點：(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ù)解，那么只取實部作業(yè)給藥方案。顯然,要設(shè)計給藥方案,必需知道給藥后血藥濃度隨時間變化的規(guī)律.為此,從實驗和實際兩方面著手.在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時辰t(小時)采集血樣,測得血藥濃度c如表:血藥濃度c(t)的測試數(shù)據(jù)t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01例5給藥方案。近似直線關(guān)系,即c(t)有按負指數(shù)規(guī)律減少的趨勢.例5給藥方案1.確定血藥濃度的變化規(guī)律假設(shè):a)藥物向體外排除的速率與中心室的血藥濃度成正比,比例系數(shù)為k(>0),稱排除速率.b)中心室血液容積為常數(shù)V,t=0瞬時注入藥物的劑量為d,血藥濃度立刻為由假設(shè)a),中心室的血藥濃度c(t)應(yīng)滿足微分方程由假設(shè)b),方程的初始條件為:求解得:即血藥濃度c(t)按指數(shù)規(guī)律下降.2.給藥方案設(shè)計簡單適用的給藥方案是:每隔一定時間,反復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并堅持在c1-c2之間.x0yc1c22.給藥方案設(shè)計簡單適用的給藥方案是:每隔一定時間,反復(fù)注入固定劑量D,使血藥濃度c(t)呈周期性變化,并堅持在c1-c2之間.為此,初次劑量需加大到D0.由式得到:顯然,當(dāng)c1,c2給定后,要確定給藥方案必需知道參數(shù)V和k.2.由實驗數(shù)據(jù)作曲線擬合以確定參數(shù)問題化為由數(shù)據(jù)ti,yi(i=1,…,8)擬合直線記為了用線性最小二乘法擬合的系數(shù)V和k,先取對數(shù)得用Matlab作線性最小二乘法擬合,得到問題化為由數(shù)據(jù)ti,yi(i=1,…,8)擬合直線記為了用線性最小二乘法擬合的系數(shù)V和k,先取對數(shù)得用Matlab作線性最小二乘法擬合,得到由實驗數(shù)據(jù)d=300(mg)算出:擬合曲線為:3.結(jié)論將k,V和給出的c1=10,c2=25代入得:D0=375.5,D=225.3,給藥方案無妨定為:D0=375mg,D=225mg,小時.范例：薄膜浸透率的測定一、問題：某種醫(yī)用薄膜，具有從高濃度的溶液向低濃度的溶液分散的功能，在試制時需測定薄膜被物質(zhì)分子穿透的才干。測定方法：用面積為S的薄膜將容器分成體積分別為的兩部份，在兩部分中分別注滿該物質(zhì)的兩種不同濃度的溶液。此時該物質(zhì)分子就會從高濃度溶液穿過薄膜向低濃度溶液中分散。平均每單位時間經(jīng)過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K表征了薄膜被該物質(zhì)分子穿透的才干，稱為浸透率。定時丈量容器中薄膜某一側(cè)的溶液濃度，以此確定K。VAVBS二、問題分析調(diào)查時段[t，t+Δt]薄膜兩側(cè)容器中該物質(zhì)質(zhì)量的變化。設(shè)，對容器的B部分溶液濃度的測試結(jié)果如下表：〔濃度單位〕1〕在容器的一側(cè)，物質(zhì)質(zhì)量的添加是由于另一側(cè)的物質(zhì)向該側(cè)浸透的結(jié)果，因此物質(zhì)質(zhì)量的增量應(yīng)等于另一側(cè)的該物質(zhì)向這側(cè)的浸透量。以容器A側(cè)為例，在時段[t，t+Δt]物質(zhì)質(zhì)量的增量為：分別表示在時辰t膜兩側(cè)溶液設(shè)的濃度，濃度單位:由于平均每單位時間經(jīng)過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)為K。因此，在時段[t，t+Δt]，從B側(cè)浸透至A側(cè)的該物質(zhì)的質(zhì)量為：于是有：兩邊除以Δt，并令Δt→0取極限再稍加整理即得：分別表示在初始時辰兩側(cè)溶液的濃度其中〔1〕2)留意到整個容器的溶液中含有該物質(zhì)的質(zhì)量不變,與初始時辰該物質(zhì)的含量一樣，因此從而：加上初值條件：代入式〔1〕得：便可得出CB(t)的變化規(guī)律，從而根據(jù)實驗數(shù)據(jù)進展擬合，估計出參數(shù)K,。三、數(shù)學(xué)模型假設(shè)：1〕薄膜兩側(cè)的溶液一直是均勻的；2〕平均每單位時間經(jīng)過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比。3〕薄膜是雙向同性的即物質(zhì)從膜的任何一側(cè)向另一側(cè)浸透的性能是一樣的。基于假設(shè)和前面的分析，B側(cè)的濃度CB(t)應(yīng)滿足如下微分方程和初始條件：四、求解方法：1.函數(shù)擬合法前面得到的模型是一個帶初值的一階線性微分方程，解之得：問題歸結(jié)為利用CB在時辰tj的丈量數(shù)據(jù)Cj(j=1,2,...,N)來辨識K和。引入從而用函數(shù)CB(t)來擬合所給的實驗數(shù)據(jù)，從而估計出其中的參數(shù)a,b,K。將代入上式有：用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)在任務(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進一步求得：2.非線性規(guī)劃法利用CB在時辰tj的丈量數(shù)據(jù)Cj(j=1,2,...,N)來辨識K和。問題可轉(zhuǎn)化為求函數(shù)即求函數(shù)的最小值點〔K，a，b〕。3.導(dǎo)函數(shù)擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1,2,...,N)假設(shè)知道一組數(shù)據(jù)那么可用最小二乘擬合的方法來求出函數(shù)中的未知參數(shù)K和h。即為求參數(shù)K,a使以下誤差函數(shù)到達最小：該問題等價于用函數(shù)f(K,a,CB)=K(0.01a-0.02CB)來擬合數(shù)據(jù)(j=1,2,...,N)用MATLAB軟件進展計算.％求數(shù)據(jù)點(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〕編寫函數(shù)M-文件baomof.mfunctionf=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中x(1)=K;x(2)=h2)編寫命令M文件(baomo21.m)3)輸出結(jié)果:x=0.10090.014即k=0.1009,h=0.014％作函數(shù)擬合x0=[0.2,0.1];x=lsqcurvefit('baomof',x0,cdata,dcj')4.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中：假設(shè)得到了參數(shù)K的一個較好的近似值K*，那么將CB(t)關(guān)于K在K*處展開，略去K的二次及以上的項得CB(t)的一個近似式經(jīng)過極小化確定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一個新的近似值，用同樣的方法再求新的修正值K。這個過程可以不斷反復(fù)，直到修正值足夠小為止