1、馬爾薩斯數據擬合實驗實驗內容實驗一：數據擬合Malthus人口指數增長模型中的參數1790-1980年間美國每隔10年的人口記錄如表5.3：表5.3年份1790180018101820183018401850人口（ X106）3.95.37.29.612.917.123.2年份1860187018801890190019101920人口（ X106）31.438.650.262.976.092.0106.5年份193019401950196019701980人口（ X10 6）123.2131.7150.7179.3204.0226.5用以上數據檢驗馬爾薩斯（Malthus）人口指數增長模型，

2、根據檢驗結果進一步討論馬爾薩斯人口模型 的改進。試驗二：經濟增長模型增加生產、發展經濟所依靠的主要因素有增加投資、增加勞動力以及技術革新等，在研究國名經濟產值 與這些因素的數量關系時，由于技術水平不像資金、勞動力那樣容易定量化，作為初步的模型，可認為技 術水平不變，只討論不像資金、勞動力之間的關系。在科學技術發展不快時，如資本主義經濟發展的前期， 這種模型是有意義的。用Q,K,L分別表示產值、資金、勞動力，要尋求的數量關系Q（K，L）。經過簡化假設與分析，在經濟 學中，推導出一個著名的Cobb-Douglas生產函數Q（K, L ）= aK aLP,0 , P 1（*）式中a, P,a要由緊急

3、統計數據確定。現有美國馬薩諸賽州1900-1926年上述三個經濟指數的統計數據。如表5.5,試用數據擬合的方法，求出（*）式中的參數a, P,a。表5.5tQKLtQKL19001.051.041.0519142.013.241.6519011.181.061.0819152.003.241.6219021.291.161.1819162.092.611.8619031.301.221.2219171.964.101.9319041.301.271.1719182.204.361.9619051.421.371.3019192.124.771.9519061.501.441.3919202.1

4、64.751.9019071.521.531.4719212.084.541.5819081.461.571.3119222.244.541.6719091.602.051.4319232.564.581.8219101.692.511.5819242.344.581.6019111.812.631.5919252.454.581.6119121.932.741.6619262.584.541.6419131.952.821.68三、實驗環境Windows操作系統；MATLAB 7.0.四、實驗過程一、問題：實驗一二、問題分析因為假設人口增長率為常數，且已知初始時刻的人口數，可根據t時刻人口變

5、化用差商代替微商。當 菩T 0時，可列出關于t時刻人口變化的微分方程。dx=TX dt。根據該微分方程求出解，再用數據擬合模型球參數。x（0 ）二 xI 0三、假設1）假設人口增長率為常數，記為r.2）記時刻t的人口為x（t）（即x（t）為模型的狀態變量）。3）假設初始時刻人口為x0。四、模型建立1）求解該微分方程得到x（t）的表達式為xD = xo *e，*t則可轉化為求關于x（t）在t時刻的測量數據來辨識參數對應的數學模型球函數的最小二乘解。其中a=0,b=x0，則需求r的最小值點。2）運用最小二函數擬合求出所需的r值。五、模型求解編寫 M 文件 curvefun.mfunction f=

6、curvefun(x,tdata)f=1e+05*39*exp(x*tdata);編寫(testl.m)tdata=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;x0=0.2;y=lsqcurvefit(curvefun,x0,tdata,cdata)即求出結果為：y = 0.0211即 r= 0.0211根據下述程序作出模型圖像：tdata=linspace(0,200,20);cdata=1e+05.*39

7、53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;plot(tdata,cdata)hold onf=1e + 05*39*exp(0.0211 *tdata);plot(tdata,f,r),title( 馬爾薩斯模型)模型圖像：六、模型評價可以看出使用一般曲線擬合，出來走勢基本一致，但是還有一定誤差，并且由圖可以看出該曲線接近一條開口向上的二次曲線。因此可以嘗試使用二次多項式對馬爾薩斯模型進行擬合。將會大大簡化求解的過程。七、模型改進即推廣則二次曲線擬合求解如下：編寫程序：tdat

8、a=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2265;p=polyfit(tdata,cdata,2)得到結果：p =1.0e+006 *0.0059-0.08265.4634編程作圖：tdata=linspace(0,200,20);cdata=1e+05.*39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232 1317 1507 1793 2040 2

9、265;plot(tdata,cdata)hold onf=1.0e+006 * 5.4634.*tdata.A2-1.0e+006 *0.0826.*tdata+1.0e+006 * 0.0059；plot(tdata,f,r),title( 馬爾薩斯模型 1)則可根據上述程序作出二次曲線與點線圖如下：可見該二次曲線與馬爾薩斯模型基本上是完全重合的。因此兒二次曲線擬合是優于一般的曲線擬合 的，對于馬爾薩斯模型來說。、問題：實驗二二、問題分析因為已知生產函數為Q(K L，且已知初始時刻的相關Q、K、L，可根據上式進行數據擬合，但是需將其轉化為2維數組計算即將L、K合成為一個2維數組data再從

10、data中將K、L取出， 則可利用一維擬合求出所要的解。三、假設Q,K,L分別表示產值、資金、勞動力認為技術水平不變，只討論不像資金、勞動力之間的關系四、模型建立五、模型求解編程如下：編寫M文件curveu.mfunction q=curvev(x,data)K=data(1,:);L=data(2,:);q=x(1).*(K.八x(2).*(L.八x(3);編寫程序(test2.m)L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160

11、 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;data=K;L;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224256 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata

12、)解得答案為：x =1.2246S4612a = 1.2246, = 0.4612,代.1276即 匕-0.1276根據程序：L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;K,L=mes

13、hgrid(K,L);Y=1.2246.*(K.A(0.4612).*(L.A(-0.1276);mesh(K,L,Y);做出圖像為：六、模型評價根據程序：L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458

14、458 454;t=1900:1926;data=K;L;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,+,t,qdata1,o)做出擬合值對比圖像如下：2.6-+2.42.2+ + + +1.81.6+ +41.41.219001905

15、19101915192019251930七、模型改進即推廣經過以上計算出的 為小于0的數，但是從對比圖像可以看出擬合效果并不理想。這是因為不同歷史時期 科技發展狀態不同，而不同科技發展水平條件下的道格拉斯函數的系數是不同的，因此我們擬合的函數需 要根據不同時期科技發展水平整合。而在不知曉科技發展水平的情況下，只能通過已知數據分析，因此我 做出下圖對比Q、K、L增長情況。t=1900:1926;L=1e-2*105 108 118 122 117 130 139 147 131 143 158 159 166 168 165 162 186 193 196 195 190 158 167 182

16、 160 161 164;K=1e-2*104 106 116 122 127 137 144 153 157 205 251 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454;qdata=1e-2*105 118 129 130 130 142 150 152 146 160 169 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;plot(t,K)hold onplot(t,L)hold onplot(t,qdata)知大致可以把數據

17、分為3段，即1900-1906,然后是1907-1910,1911-1926,3段，再次做擬合:編寫M文件curveu.m function q=curvev(x,data) K=data(1,:);L=data(2,:);q=x(1).*(K.八x(2).*(L.八x(3);編寫程序(test2.m)i . L=1e-2*105 108 118 122 117 130 139;K=1e-2*104 106 116 122 127 137 144; data=K;L;qdata=1e-2*105 118 129 130 130 142 150; x0=0.2,0.5,0.5;x=lsqcurv

18、efit(curvev,x0,data,qdata) 解之得：x 1=1.06100.31010.7404ii.L=1e-2* 147 131 143 158;K=1e-2* 153 157 205 251;data=K;L;qdata=1e-2* 152 146 160 169;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata) 解之得：x 2=1.24240.17950.3215iii.L=1e-2* 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2

19、* 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454; data=K;L;qdata=1e-2* 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258; x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之的：x3 =1.48000.4116-0.3425又由上x1,x2，已滿足條件，而x3不滿足，則可對x3段再進行分段。則觀察1911-1926圖像如下：t=1911:1926;L=1e-2*

20、 159 166 168 165 162 186 193 196 195 190 158 167 182 160 161 164;K=1e-2* 263 274 282 324 324 361 410 436 477 475 454 454 458 458 458 454; qdata=1e-2* 181 193 195 201 200 209 196 220 212 216 208 224 256 234 245 258;plot(t,K)hold onplot(t,L)hold onplot(t,qdata)知需將 1910-1926 再分為 1911-1920,1921-1923,192

21、4-1926三段，再次擬合:L=1e-2* 159 166 168 165 162 186 193 196 195 190;K=1e-2* 263 274 282 324 324 361 410 436 477 475;data=K;L;qdata=1e-2* 181 193 195 201 200 209 196 220 212 216;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之得：x 31=1.53370.21700.0023又因為多次實驗，1920-1926這段時間內的點并不一定滿足所求的道格拉斯公式，因此將其視為奇異點。L=

22、1e-2* 158 167 182;K=1e-2* 454 454 458;data=K;L;qdata=1e-2* 208 224 256;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)L=1e-2* 160 161 164;K=1e-2* 458 458 454;data=K;L;qdata=1e-2* 234 245 258;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)解之得：x32 =0.01383.18210.6371則綜上可分為3段函數擬合，所求出的a,a, P如下:a =

23、 1.0610a1 = 0.31011 (1900 -1906)c1p = 0.740412 =a2P2a3a3p3T.2424=0.1795(1907 -1910)=0.3215=1.5337=0.2170 ”911-1920)=0.0023因為1919以后因為科技發展加快則道格拉斯模型不再適用，因此則會有1919之后的點為奇異點，與道格 拉斯模型擬合并不符合。則再次根據下列程序t=1900:1906L=1e-2*105 108 118 122 117 130 139;K=1e-2*104 106 116 122 127 137 144;data=K;L;qdata=1e-2*105 118 129 130 130 142 150;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,data,qdata)qdata1=curveu(x,data);plot(t,qdata,+,t,qdata1,o)hold ont=1907:1910L=1e-2* 147 131 143 158;K=1e-2* 153 157 205 251;data=K;L;qdata=1e-2* 152 146 160 169;x0=0.2,0.5,0.5;x=lsqcurvefit(curvev,x0,d