用MATLAB進(jìn)行數(shù)據(jù)擬合PPT學(xué)習(xí)教案

1、會(huì)計(jì)學(xué)1用用MATLAB進(jìn)行數(shù)據(jù)擬合進(jìn)行數(shù)據(jù)擬合擬擬 合合2.2.擬合的基本原理擬合的基本原理1. 擬合問(wèn)題引例擬合問(wèn)題引例第1頁(yè)/共65頁(yè)擬擬 合合 問(wèn)問(wèn) 題題 引引 例例 1 1溫度溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻電阻R( ) 765 826 873 942 1032已知熱敏電阻數(shù)據(jù)已知熱敏電阻數(shù)據(jù)：求求60600C時(shí)的電阻時(shí)的電阻R。2040608010070080090010001100 設(shè)設(shè) R=at+ba,b為待定系數(shù)為待定系數(shù)第2頁(yè)/共65頁(yè)擬擬 合合 問(wèn)問(wèn) 題題 引引 例例 2 2 t (h) 0.25 0.5 1 1.5 2 3 4 6 8

2、c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射注射300mg)求血藥濃度隨時(shí)間的變化規(guī)律求血藥濃度隨時(shí)間的變化規(guī)律c(t).作半對(duì)數(shù)坐標(biāo)系作半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形下的圖形為待定系數(shù)kcectckt,)(002468100101102MATLAB(aa1)第3頁(yè)/共65頁(yè)曲曲 線(xiàn)線(xiàn) 擬擬 合合 問(wèn)問(wèn) 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n個(gè)點(diǎn)個(gè)點(diǎn)（xi,yi) i=1,

3、n, 尋求一個(gè)函數(shù)（曲線(xiàn)）尋求一個(gè)函數(shù)（曲線(xiàn)）y=f(x), 使使 f(x) 在某種準(zhǔn)則下與所在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線(xiàn)擬合得最好。有數(shù)據(jù)點(diǎn)最為接近，即曲線(xiàn)擬合得最好。 +xyy=f(x)(xi,yi)i i 為點(diǎn)為點(diǎn)（xi,yi) 與與曲線(xiàn)曲線(xiàn) y=f(x) 的距離的距離第4頁(yè)/共65頁(yè)擬合與插值的關(guān)系擬合與插值的關(guān)系 函數(shù)插值與曲線(xiàn)擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作函數(shù)插值與曲線(xiàn)擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。的。 實(shí)例：實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到

4、X和 f之間的關(guān)系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)問(wèn)題：?jiǎn)栴}：給定一批數(shù)據(jù)點(diǎn)，需確定滿(mǎn)足特定要求的曲線(xiàn)或曲面解決方案：解決方案：若不要求曲線(xiàn)（面）通過(guò)所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，這就是數(shù)據(jù)擬合數(shù)據(jù)擬合，又稱(chēng)曲線(xiàn)擬合或曲面擬合。若要求所求曲線(xiàn)（面）通過(guò)所給所有數(shù)據(jù)點(diǎn)，就是插值問(wèn)題插值問(wèn)題；第5頁(yè)/共65頁(yè)最臨近插值、線(xiàn)性插值、樣條插值與曲線(xiàn)擬合結(jié)果最臨近插值、線(xiàn)性插值、樣條插值與曲線(xiàn)擬合結(jié)果：0246810121416180510152025已知數(shù)據(jù)點(diǎn)spl

5、ine三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)linest三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)nearest三次多項(xiàng)式插值第6頁(yè)/共65頁(yè)曲線(xiàn)擬合問(wèn)題最常用的解法曲線(xiàn)擬合問(wèn)題最常用的解法線(xiàn)性最小二乘法的基本思路線(xiàn)性最小二乘法的基本思路第一步: :先選定一組函數(shù)先選定一組函數(shù) r1(x), r2(x), rm(x), mn, 令令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1）其中其中 a1,a2, am 為待定系數(shù)。為待定系數(shù)。 第二步: 確定確定a1,a2, am 的準(zhǔn)則（最小二乘準(zhǔn)則）：的準(zhǔn)則（最

6、小二乘準(zhǔn)則）：使使n個(gè)點(diǎn)個(gè)點(diǎn)（xi,yi) 與與曲線(xiàn)曲線(xiàn) y=f(x) 的距離的距離 i 的平方和最小的平方和最小 。記記 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 問(wèn)題歸結(jié)為，求問(wèn)題歸結(jié)為，求 a1,a2, am 使使 J(a1,a2, am) 最小。最小。第7頁(yè)/共65頁(yè)線(xiàn)性最小二乘法的求解：預(yù)備知識(shí)線(xiàn)性最小二乘法的求解：預(yù)備知識(shí)超定方程組超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組)( 221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaarrrrr

7、rR112111211,其中其中超定方程一般是不存在解的矛盾方程組。超定方程一般是不存在解的矛盾方程組。 如果有向量如果有向量a使得使得 達(dá)到最小，達(dá)到最小，則稱(chēng)則稱(chēng)a為上述為上述超定方程的最小二乘解超定方程的最小二乘解。 212211)(imniimiiyararar第8頁(yè)/共65頁(yè)線(xiàn)性最小二乘法的求解線(xiàn)性最小二乘法的求解 定理：定理：當(dāng)當(dāng)R RT TR R可逆時(shí)，超定方程組（可逆時(shí)，超定方程組（3 3）存在最小二乘解）存在最小二乘解，且即為方程組，且即為方程組 R RT TRa=RRa=RT Ty y的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y 所以，曲線(xiàn)擬合的最小

8、二乘法要解決的問(wèn)題，實(shí)際上就是所以，曲線(xiàn)擬合的最小二乘法要解決的問(wèn)題，實(shí)際上就是求以下超定方程組的最小二乘解的問(wèn)題。求以下超定方程組的最小二乘解的問(wèn)題。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）第9頁(yè)/共65頁(yè)線(xiàn)性最小二乘擬合線(xiàn)性最小二乘擬合 f(x)=a1r1(x)+ +amrm(x)中函數(shù)中函數(shù)rr1 1(x), (x), r rm m(x)(x)的選取的選取 1. 1. 通過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定通過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/

9、xf=aebxf=ae-bx 2. 2. 將數(shù)據(jù)將數(shù)據(jù) (xi,yi) i=1, n 作圖，通過(guò)直觀判斷確定作圖，通過(guò)直觀判斷確定 f(x)：第10頁(yè)/共65頁(yè)用用MATLAB解擬合問(wèn)題解擬合問(wèn)題1 1、線(xiàn)性最小二乘擬合、線(xiàn)性最小二乘擬合2 2、非線(xiàn)性最小二乘擬合、非線(xiàn)性最小二乘擬合第11頁(yè)/共65頁(yè)用用MATLAB作線(xiàn)性最小二乘擬合作線(xiàn)性最小二乘擬合1. 1. 作多項(xiàng)式作多項(xiàng)式f(x)=a1xm+ +amx+am+1擬合擬合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2. 2. 對(duì)超定方程組對(duì)超定方程組)(11nmyaRnmmn可得最小二乘意義下的解。可得最小二乘意義

10、下的解。，用用yRa3.3.多項(xiàng)式在多項(xiàng)式在x x處的值處的值y y可用以下命令計(jì)算：可用以下命令計(jì)算： y=polyvaly=polyval（a a，x x）輸出擬合多項(xiàng)式系數(shù)輸出擬合多項(xiàng)式系數(shù)a=a1, am , am+1 (數(shù)組數(shù)組) )）輸入同長(zhǎng)度輸入同長(zhǎng)度的數(shù)組的數(shù)組X，Y擬合多項(xiàng)擬合多項(xiàng)式次數(shù)式次數(shù)第12頁(yè)/共65頁(yè)即要求即要求 出二次多項(xiàng)式出二次多項(xiàng)式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例1 對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合x(chóng)i0.10.20.40.50.60.70.80.91yi1.9

11、783.286.167.347.669.589.489.3011.2第13頁(yè)/共65頁(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=Ry11 11211121xxxxR此時(shí)MATLAB(zxec1)解法解法1 1用解超定方程的方法用解超定方程的方法2）計(jì)算結(jié)果）計(jì)算結(jié)果： = -9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf第14頁(yè)/共65頁(yè)第15頁(yè)/共65頁(yè)1）輸入以

12、下命令：）輸入以下命令： 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)和擬合曲線(xiàn)的圖形作出數(shù)據(jù)點(diǎn)和擬合曲線(xiàn)的圖形2）計(jì)算結(jié)果：）計(jì)算結(jié)果： = -9.8108 20.1293 -0.0317解法解法2用多項(xiàng)式擬合的命令用多項(xiàng)式擬合的命令MATLAB(zxec2)00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf第16頁(yè)/共65頁(yè)例

13、例2 已知觀測(cè)數(shù)據(jù)點(diǎ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)次多項(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.2plot(x,y,k.,markersize,25)axis(0 1.3 -2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)編寫(xiě)編寫(xiě)Matlab程序如下程

14、序如下:第17頁(yè)/共65頁(yè)t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)hold onplot(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.2plot(x,y,k.,markersize,25)axis(0 1.3 -2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)第18頁(yè)/共65頁(yè)例例3 用切削機(jī)床進(jìn)行金屬品加工時(shí)用切削機(jī)床進(jìn)行金屬品加工時(shí), 為

15、了適當(dāng)?shù)卣{(diào)整為了適當(dāng)?shù)卣{(diào)整機(jī)床機(jī)床, 需要測(cè)定刀具的磨損速度需要測(cè)定刀具的磨損速度. 在一定的時(shí)間測(cè)量刀在一定的時(shí)間測(cè)量刀具的厚度具的厚度, 得數(shù)據(jù)如表所示得數(shù)據(jù)如表所示:切削時(shí)間切削時(shí)間 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削時(shí)間切削時(shí)間 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm第19頁(yè)/共65頁(yè)解解: 描出散點(diǎn)圖描出散點(diǎn)圖, 在命令窗口輸入在命令窗口輸入:t=0:1:16y=30.0 29.1 28.4 28.1 28

16、.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)第20頁(yè)/共65頁(yè)解解: 描出散點(diǎn)圖描出散點(diǎn)圖, 在命令窗口輸入在命令窗口輸入:t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)a = -0.3012 29.3804hold onplot(t, y1), hold offa=polyfit(t,y,1)y1=-0.3012*t+29.380

17、4第21頁(yè)/共65頁(yè)例例3 用切削機(jī)床進(jìn)行金屬品加工時(shí)用切削機(jī)床進(jìn)行金屬品加工時(shí), 為了適當(dāng)?shù)卣{(diào)整為了適當(dāng)?shù)卣{(diào)整機(jī)床機(jī)床, 需要測(cè)定刀具的磨損速度需要測(cè)定刀具的磨損速度. 在一定的時(shí)間測(cè)量刀在一定的時(shí)間測(cè)量刀具的厚度具的厚度, 得數(shù)據(jù)如表所示得數(shù)據(jù)如表所示:切削時(shí)間切削時(shí)間 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削時(shí)間切削時(shí)間 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm擬合曲線(xiàn)為擬合曲線(xiàn)為: y=-0.3012t+29.380

18、4第22頁(yè)/共65頁(yè)例例4 一個(gè)一個(gè)15.4cm30.48cm的混凝土柱在加壓實(shí)驗(yàn)中的的混凝土柱在加壓實(shí)驗(yàn)中的應(yīng)力應(yīng)力-應(yīng)變關(guān)系測(cè)試點(diǎn)的數(shù)據(jù)如表所示應(yīng)變關(guān)系測(cè)試點(diǎn)的數(shù)據(jù)如表所示1.55 2/ N/m 2/ / N/m 6500 10 33.103 10 2.4761000 10 32.465 10 2. 9361500 10 31.953 10 3. 0362000 10 31.517 10 已知應(yīng)力已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述, 即假設(shè)即假設(shè)21 kke式中式中, 表示應(yīng)力表示應(yīng)力, 單位是單位是 N/m2; 表示應(yīng)變表示應(yīng)變. 2.896237

19、5 10 31.219 10 第23頁(yè)/共65頁(yè)已知應(yīng)力已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述, 即假設(shè)即假設(shè)21 kke式中式中, 表示應(yīng)力表示應(yīng)力, 單位是單位是 N/m2; 表示應(yīng)變表示應(yīng)變. 解解 選取選取指數(shù)函數(shù)指數(shù)函數(shù)作擬合時(shí)作擬合時(shí), 在擬合前需作變量代換在擬合前需作變量代換, 化為化為 k1, k2 的線(xiàn)性函的線(xiàn)性函數(shù)數(shù).于是于是,12lnln kk令令0211ln,ln zakak即即01 zaa第24頁(yè)/共65頁(yè)在命令窗口輸入在命令窗口輸入:x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*1.0e-6

20、2375*1.0e-6y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=1.55 2.47 2.93 3.03 2.89plot(x,w,*)y1=exp(8.3009)*x.*exp( -494.5209*x)plot(x,w,*,x,y1,r-)第25頁(yè)/共65頁(yè)已知應(yīng)力已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述應(yīng)變關(guān)系可以用一條指數(shù)曲線(xiàn)來(lái)描述, 即假設(shè)即假設(shè)21 kke式中式中, 表示應(yīng)力表示應(yīng)力, 單位是單位是 N/m2; 表

21、示應(yīng)變表示應(yīng)變. 擬合曲線(xiàn)為擬合曲線(xiàn)為:3 -494.52094.0275 10 e0211ln,ln, zakak01 zaa0211-494.5209,ln8.3009, akak3124.0275 10 ,494.5209kk令令則則求求得得于于是是第26頁(yè)/共65頁(yè)在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線(xiàn)有在實(shí)際應(yīng)用中常見(jiàn)的擬合曲線(xiàn)有:01ya xa直直線(xiàn)線(xiàn)101 nnnya xa xa多項(xiàng)式多項(xiàng)式一般一般 n=2, 3, 不宜過(guò)高不宜過(guò)高.01ayax雙曲線(xiàn)雙曲線(xiàn)(一支一支) bxyae指數(shù)曲指數(shù)曲線(xiàn)線(xiàn)第27頁(yè)/共65頁(yè)27754,1521 aa從而得到人口數(shù)與年份的函數(shù)關(guān)系為從而得到人口數(shù)與年份

22、的函數(shù)關(guān)系為2775415 xy把把x=1999代如，估算出代如，估算出1999年的人口數(shù)為年的人口數(shù)為 y=1252.1（百萬(wàn)）（百萬(wàn)）12.52億億1999年實(shí)際人口數(shù)量為年實(shí)際人口數(shù)量為.億。億。線(xiàn)性預(yù)測(cè)模型線(xiàn)性預(yù)測(cè)模型第28頁(yè)/共65頁(yè) 英國(guó)人口學(xué)家英國(guó)人口學(xué)家MalthusMalthus根據(jù)百余年的人口統(tǒng)計(jì)資根據(jù)百余年的人口統(tǒng)計(jì)資料，于料，于17981798年提出了著名的人口自然增長(zhǎng)的年提出了著名的人口自然增長(zhǎng)的指數(shù)增指數(shù)增長(zhǎng)模型長(zhǎng)模型。 btaetx 0)0(,xxrxdtdxt rextx0)( 基本假設(shè)基本假設(shè) : : 人口人口( (相對(duì)相對(duì)) )增長(zhǎng)率增長(zhǎng)率 r 是常數(shù)是常數(shù)

23、設(shè)設(shè)x(t) 為為t時(shí)刻時(shí)刻的的人口人口數(shù)量數(shù)量,t=0時(shí)人口數(shù)為時(shí)人口數(shù)為x0 btatx ln指數(shù)增長(zhǎng)模型指數(shù)增長(zhǎng)模型實(shí)際中，常用實(shí)際中，常用第29頁(yè)/共65頁(yè)1. 1. 由前由前100100年的數(shù)據(jù)求出美國(guó)的人口增長(zhǎng)年的數(shù)據(jù)求出美國(guó)的人口增長(zhǎng)Malthus模型模型。2. 2. 預(yù)測(cè)后預(yù)測(cè)后100100年（每隔年（每隔1010年）的人口狀況。年）的人口狀況。3. 3. 根據(jù)預(yù)測(cè)的人口狀況和實(shí)際的人口數(shù)量根據(jù)預(yù)測(cè)的人口狀況和實(shí)際的人口數(shù)量, ,討論人討論人口模型的改進(jìn)情況?？谀Ｐ偷母倪M(jìn)情況。美國(guó)美國(guó)17901790年年19801980年每隔年每隔1010年的人口記錄年的人口記錄226.520

24、4.0179.3150.7131.7123.2106.592.076.062.9人口人口( (百萬(wàn)百萬(wàn)) )1980197019601950194019301920191019001890年份年份50.238.631.423.217.112.99.67.25.33.9人口人口( (百萬(wàn)百萬(wàn)) )1880187018601850184018301820181018001790年份年份例例第30頁(yè)/共65頁(yè)解：解： btaetx 設(shè)設(shè)btax ln2101)ln(),(iiixbtabaJ 取得最小值取得最小值. . 其中其中, , iitxx 表示人口數(shù)量表示人口數(shù)量。it表示年份表示年份, ,

25、解方程組解方程組: 0ln210ln10211011012101101101iiiiiiiiiiitxbtatbJxbtaaJ即得參數(shù)即得參數(shù)ba,的值的值.ba,使得使得問(wèn)題轉(zhuǎn)化為求參數(shù)問(wèn)題轉(zhuǎn)化為求參數(shù)第31頁(yè)/共65頁(yè) prog41.m %prog41.m % This program is to predict the number of % This program is to predict the number of population %population %format longformat longt1=1790;1800;1810;1820;1830;1840;1850

26、;1860;18t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;70;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;19t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;70;1980;x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;5x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;0.2;x2=62.9;76.0;92.0;106.5;123

27、.2;131.7;150.7;179x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;.3;204.0;226.5;lnx1=log(x1); lnx2=log(x2);lnx1=log(x1); lnx2=log(x2);第32頁(yè)/共65頁(yè)第33頁(yè)/共65頁(yè)a= -49.79535457790735b=0.02859807120038仿真結(jié)果表明仿真結(jié)果表明： 人口增加的人口增加的指數(shù)模型在短指數(shù)模型在短期內(nèi)基本上能期內(nèi)基本上能比較準(zhǔn)確地反比較準(zhǔn)確地反映人口自然增映人口自然增長(zhǎng)的規(guī)律，但長(zhǎng)的規(guī)律，但長(zhǎng)期預(yù)測(cè)誤差長(zhǎng)期預(yù)測(cè)誤差

28、很大，需要修很大，需要修正預(yù)測(cè)模型。正預(yù)測(cè)模型。擬合曲線(xiàn)擬合曲線(xiàn)原始數(shù)據(jù)曲線(xiàn)原始數(shù)據(jù)曲線(xiàn)第34頁(yè)/共65頁(yè)人口增長(zhǎng)到一定數(shù)量后，增長(zhǎng)率下降的原因：人口增長(zhǎng)到一定數(shù)量后，增長(zhǎng)率下降的原因：資源、環(huán)境等因素對(duì)人口增長(zhǎng)的阻滯作用資源、環(huán)境等因素對(duì)人口增長(zhǎng)的阻滯作用且阻滯作用隨人口數(shù)量增加而變大且阻滯作用隨人口數(shù)量增加而變大假設(shè)假設(shè)r固有增長(zhǎng)率固有增長(zhǎng)率(x很小時(shí)很小時(shí))k人口容量（資源、環(huán)境能容納的最大數(shù)量）人口容量（資源、環(huán)境能容納的最大數(shù)量）)1()(kxrxr r是是x的減函數(shù)的減函數(shù)krs 0)( kr)1()(mxxrxxxrdtdx 000)()(xexkkxtxtr 例的例的Logis

29、ticLogistic模型留給同學(xué)們練習(xí)模型留給同學(xué)們練習(xí)第35頁(yè)/共65頁(yè)% This program is to predict the model of population by % This program is to predict the model of population by 4-degree polynomial% 4-degree polynomial% %prog42.m%prog42.m%format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t1=1790;1800;1

30、810;1820;1830;1840;1850;1860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t=t1;t2;t=t1;t2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7

31、;179.3;204.0;2P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;26.5;P=P1;P2;P=P1;P2;n=4; n=4; % The degree of the fitting polynomial% The degree of the fitting polynomial%a,s=polyfit(t1,P1,n);a,s=polyfit(t1,P1,n);y=polyval(a,t);y=polyval(a,t);% % a a is the coefficients vector from n-degr

32、ee to 0- is the coefficients vector from n-degree to 0-degree%degree%plot(t,P,rplot(t,P,r* *-,t,y,b+-);-,t,y,b+-);23第36頁(yè)/共65頁(yè)a =1.0e+006 * -0.00000000000014 0.00000000107892 -0.00000304878595 0.00381927346813 -1.79012132225427仿真結(jié)果表明仿真結(jié)果表明, , 人口增加的模型用多項(xiàng)式擬合能人口增加的模型用多項(xiàng)式擬合能比較準(zhǔn)確地反映人口自然增長(zhǎng)的規(guī)律，對(duì)長(zhǎng)期預(yù)比較準(zhǔn)確地反映人口

33、自然增長(zhǎng)的規(guī)律，對(duì)長(zhǎng)期預(yù)測(cè)具有指導(dǎo)意義。測(cè)具有指導(dǎo)意義。第37頁(yè)/共65頁(yè)例例2: 2: 海底光纜線(xiàn)長(zhǎng)度預(yù)測(cè)模型海底光纜線(xiàn)長(zhǎng)度預(yù)測(cè)模型某一通信公司在一次某一通信公司在一次施工中施工中, ,需要在水面寬需要在水面寬為為20m20m的河溝底沿直線(xiàn)的河溝底沿直線(xiàn)走向鋪設(shè)一條溝底光走向鋪設(shè)一條溝底光纜纜. .在鋪設(shè)光纜之前需在鋪設(shè)光纜之前需要對(duì)溝底的地形做初要對(duì)溝底的地形做初Bix2468101214161820986420ADCih0 x20 x探測(cè)到一組等分點(diǎn)位置的深度數(shù)據(jù)如下表所示探測(cè)到一組等分點(diǎn)位置的深度數(shù)據(jù)如下表所示. .25步探測(cè)步探測(cè), ,從而估計(jì)所需光纜的長(zhǎng)度從而估計(jì)所需光纜的長(zhǎng)度,

34、,為工程預(yù)算為工程預(yù)算提供依據(jù)提供依據(jù). .基本情況如圖所示基本情況如圖所示. .第38頁(yè)/共65頁(yè)10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度深度(m)9876543210分點(diǎn)分點(diǎn)2121個(gè)等分點(diǎn)處的深度個(gè)等分點(diǎn)處的深度(1) (1) 預(yù)測(cè)通過(guò)這條河溝所需光纜長(zhǎng)度的近似值預(yù)測(cè)通過(guò)這條河溝所需光纜長(zhǎng)度的近似值. .(2) (2) 作出鋪設(shè)溝底光纜的曲線(xiàn)圖作出鋪設(shè)溝底光纜的曲線(xiàn)圖. .第39頁(yè)/共65頁(yè)02

35、4681012141618207891011121314解：解： 用用12次多項(xiàng)式函數(shù)擬合光纜走勢(shì)的曲線(xiàn)圖如下次多項(xiàng)式函數(shù)擬合光纜走勢(shì)的曲線(xiàn)圖如下仿真結(jié)果表仿真結(jié)果表明明,擬合曲擬合曲線(xiàn)能較準(zhǔn)確線(xiàn)能較準(zhǔn)確地反映光纜地反映光纜的走勢(shì)圖的走勢(shì)圖.The length of the label is L= 26.3809 (m)假設(shè)所鋪設(shè)的光纜足夠柔軟假設(shè)所鋪設(shè)的光纜足夠柔軟,在鋪設(shè)過(guò)程中光纜觸地在鋪設(shè)過(guò)程中光纜觸地走勢(shì)光滑走勢(shì)光滑,緊貼地面緊貼地面,并且忽略水流對(duì)光纜的沖擊并且忽略水流對(duì)光纜的沖擊.第40頁(yè)/共65頁(yè)% prog45.m This program is to fit the dat

36、a by polynomial % prog45.m This program is to fit the data by polynomial %format longformat longt=linspace(0,20,21);t=linspace(0,20,21);x=linspace(0,20,100);x=linspace(0,20,100);P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,1

37、1.29,10.22,9.15,7.90,7.95,8.86,9.81,1028,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93;.80,10.93;a,s=polyfit(t,P,12);a,s=polyfit(t,P,12);yy=polyval(a,x);yy=polyval(a,x);disp(yy=);disp(yy);disp(yy=);disp(yy);plot(x,yy,rplot(x,yy,r* *-,t,P,b+-);-,t,P,b+-);L=0;L=0;forfor i=2:100 i=2:100

38、L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2); L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);endenddisp(The length of the label is L=);disp(L);disp(The length of the label is L=);disp(L);第41頁(yè)/共65頁(yè)1. lsqcurvefit1. lsqcurvefit已知數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)： xdata=xdata=（xdata1，xdata2，xdataxdatan n），）， ydata=ydata=（ydataydata1 1，ydatayda

39、ta2 2，ydataydatan n） 用用MATLAB作非線(xiàn)性最小二乘擬合作非線(xiàn)性最小二乘擬合 Matlab Matlab的提供了兩個(gè)求非線(xiàn)性最小二乘擬合的函數(shù)：的提供了兩個(gè)求非線(xiàn)性最小二乘擬合的函數(shù)：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。兩個(gè)命令都要先建立。兩個(gè)命令都要先建立M-M-文件文件fun.mfun.m，在其中定義函數(shù)，在其中定義函數(shù)f(x)f(x)，但兩者定義，但兩者定義f(x)f(x)的方式是不同的的方式是不同的, ,可可參考例題參考例題.最小 ),(21niiiydataxdataxF lsqcurvefitlsqcurvefi

40、t用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù)F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F(xiàn) F（x x，xdataxdatan n）T T中的參變量中的參變量x(x(向量向量),),使得使得 第42頁(yè)/共65頁(yè) 輸入格式為輸入格式為: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefit (fun,x0,xdata,ydata,options,g

41、rad); (4) x, options = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, options,funval = lsqcurvefit (fun,x0,xdata,ydata,); (6) x, options,funval, Jacob = lsqcurvefit (fun,x0,xdata,ydata,);fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)F(x,xdata) 的的M-文件文件, 自變量為自變量為x和和xdata說(shuō)明：x = lsqcurvefit (fun,x0,xdata,ydata,options);迭代初值迭代

42、初值已知數(shù)據(jù)點(diǎn)已知數(shù)據(jù)點(diǎn)選項(xiàng)見(jiàn)無(wú)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化約束優(yōu)化第43頁(yè)/共65頁(yè) lsqnonlin用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù) f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的參量中的參量x x，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22221)()()()()(xfxfxfxfxfnT2. lsqnonlin已知數(shù)據(jù)點(diǎn)：已知數(shù)據(jù)點(diǎn)： xdata=xdata=（xdata1，xdata2，xdataxdat

43、an n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）第44頁(yè)/共65頁(yè)輸入格式為：輸入格式為： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0，options）； 3） x= lsqnonlin （fun，x0，options，grad）； 4） x，options= lsqnonlin （fun，x0，）； 5） x，options，funval= lsqnonlin （fun，x0，）；說(shuō)明：x= lsqnonlinlsqnonlin （fun，x0，options）；）；

44、fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)f(x)的的M-文件文件，自變量為自變量為x迭代初值迭代初值選項(xiàng)見(jiàn)無(wú)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化約束優(yōu)化第45頁(yè)/共65頁(yè)例例6 已知觀測(cè)數(shù)據(jù)點(diǎ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ù)求三個(gè)參數(shù) a, b, c的值的值, 使得曲線(xiàn)使得曲線(xiàn) f(x)=aex+bx2+cx3 與與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.首先編寫(xiě)存儲(chǔ)擬合函數(shù)的函數(shù)文件首先編寫(xiě)存儲(chǔ)擬合函數(shù)的函數(shù)文件.fu

45、nction f=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3保存為文件保存為文件 nihehanshu.m第46頁(yè)/共65頁(yè)例例6 已知觀測(cè)數(shù)據(jù)點(diǎ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ù)求三個(gè)參數(shù) a, b, c的值的值, 使得曲線(xiàn)使得曲線(xiàn) f(x)=aex+bx2+cx3 與與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.編寫(xiě)下面的程序調(diào)用擬合函

46、數(shù)編寫(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)第47頁(yè)/共65頁(yè)編寫(xiě)下面的程序調(diào)用擬合函數(shù)編寫(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,xda

47、ta,ydata)程序運(yùn)行后顯程序運(yùn)行后顯示示x = 3.0022 4.0304 0.9404resnorm = 0.0912第48頁(yè)/共65頁(yè)例例6 已知觀測(cè)數(shù)據(jù)點(diǎ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ù)求三個(gè)參數(shù) a, b, c的值的值, 使得曲線(xiàn)使得曲線(xiàn) f(x)=aex+bx2+cx3 與與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.說(shuō)明說(shuō)明: 最小二乘意義上的最佳擬合函數(shù)為最小二乘意義上的最佳擬合函數(shù)為f(x)= 3

48、ex+ 4.03x2 + 0.94 x3.此時(shí)的殘差是此時(shí)的殘差是: 0.0912.第49頁(yè)/共65頁(yè)f(x)= 3ex+ 4.03x2 + 0.94 x3.擬合函數(shù)為擬合函數(shù)為:第50頁(yè)/共65頁(yè)最小二乘擬合：最小二乘擬合：lsqcurvefit()例例7 假設(shè)有一組實(shí)測(cè)數(shù)據(jù)假設(shè)有一組實(shí)測(cè)數(shù)據(jù)xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 yi 2.3201 2.647 2.9707 3.2885 3.6008 3.909 4.2147 4.5191 4.8232 5.1275 假設(shè)已知該數(shù)據(jù)可能滿(mǎn)足的原型函數(shù)為假設(shè)已知該數(shù)據(jù)可能滿(mǎn)足的原型函數(shù)為debx

49、axycx2試求出滿(mǎn)足下面數(shù)據(jù)的最小二乘解的試求出滿(mǎn)足下面數(shù)據(jù)的最小二乘解的a,b,c,d的值的值。 第51頁(yè)/共65頁(yè)100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc210102. 0),(minjjktcbeakbaFj 例例8 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a，b，kktbeatc2 . 0 . 0)(該問(wèn)題即解最優(yōu)化問(wèn)題：該問(wèn)題即解最優(yōu)化問(wèn)題：第52頁(yè)/共65頁(yè)MATLAB(fzxec1) 1 1）編寫(xiě)）編寫(xiě)M-M-文件文件 curve

50、fun1.mcurvefun1.m function f=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:1000tdata=100:100:1000cdata=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0

51、.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tktktbeabea),(10102. 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit第53頁(yè)/共65頁(yè)3 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059

52、0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542 x = 0.0063 -0.0034 0.25424）結(jié)論）結(jié)論：a=0.0063, b=-0.0034, k=0.2542第54頁(yè)/共65頁(yè)MATLAB(fzxec2)Tktktcbeacbea),(102. 0102. 0101 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）1）編寫(xiě)編寫(xiě)M-M-文件文件 curvefun2.

53、mcurvefun2.m function f=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ù)函數(shù)curvefun2的自變量是的自變量是x，cdata和和tdata是已知參數(shù)，故應(yīng)是已知參數(shù)，故應(yīng)將將cdata tdata的值寫(xiě)在的值寫(xiě)在curvefun2.m中中第55頁(yè)/共65頁(yè)3 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為 f =1.0e-003 f =1.0e-003 * *(0.2322 -0.1243 -0.2495 -0.2413 (0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0