马尔萨斯数据拟合实验

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