数模竞赛2007A.doc

中国人口增长预测摘 要：本问题是一个已经相当成熟的数学建模问题，然而，它始终未有定论。我们曾试图遵循前人的足迹，想站在巨人的肩膀上看得更深刻、更广阔，然而，确实深邃的空间越容易使人迷惑。最终，我们决定回到起点，另辟蹊径。我们从人口递推关系：上年人口数+出生数-死亡数迁入迁出数=下年人口数出发，通过上一年的人口总数、各年龄男女比例和死亡率、生育率、城乡迁移率、出生性别比可以方便的解出下一年的人口总数、各年龄男女比例，以此为依据开始建立模型。我们通过设定死亡率、生育率、城乡迁移率和出生性别比四类基本参数，建立起每年的人口特征值（人口总数、各年龄男女比例）的递推模型。而这些基本参数恰恰能深刻反映出如老龄化进程加速、出生人口性别比持续升高，乡村人口城镇化等国情对人口增长及分布的影响。参数估计是该模型精确与否的关键。我们以例证为先判断已知年份死亡率随年龄分布的曲线之间的相似度，认为在中国基本国情和人口国策无重大变动的情况下同一年龄同一类型人口的死亡率维持在均衡水平，再得出6条平均死亡率Expotential曲线。同样的道理，可以得出3条平均生育率Gaussian曲线。城乡迁移率先通过已知数据求得，再同出生性别比一样，进行变化趋势估计，最终认定它们依然处于一个较稳定的波动进程中。利用本模型计算2002年到2060年的全国总人口数，主要数据如下：年份200220032004200520062007200820092010总人口数/万128770 129440130150130890131670132490133330134200135100年份201520202025203020352040204520502055总人口数/万139730144000147240148960148720146140141540136050130060由预测人口变化趋势图可得出以下结论：1、我国总人口在中短期内仍将缓慢增长；2、在2032年左右全国人口总数将达到峰值149130万人；3、之后人口将会转为下降趋势并会持续较长一段时期。结论中2001年到2001年的预测人口数与实际值相比较，误差都在0.25%之内，可见此模型在中短期预测上具有相当高的准确度。模型解答过程中，对死亡率、生育率的曲线拟合效果相当好，极大地减小了误差的累计。关键词：人口递推关系、人口特征值、四类基本参数、人口降低一 问题的重述中国是一个人口大国，人口问题始终是制约我国发展的关键因素之一。根据已有数据，运用数学建模的方法，对中国的人口做出分析和预测是一个十分重要的问题。近年来中国的人口发展出现了一些新的特点，例如，老龄化进程加速、出生人口性别比持续升高，以及乡村人口城镇化等因素，这些都影响着中国人口的增长。对于中国人口问题已有多方面的研究，并积累了大量数据资料。题目给出了2001年到2005年城镇乡男女各年龄人口所占比例及死亡率、2001年到2005年城镇乡男女之间的比例、2001年到2005年城镇乡育龄期各年龄妇女的生育率、1994年到2005年城镇乡出生性别比、1995年到2005年城乡镇育龄期妇女生育率。问题要求从中国的实际情况和人口增长的上述特点出发，参考题目给出的数据（也可以搜索相关文献和补充新的数据），建立中国人口增长的数学模型，并由此对中国人口增长的中短期和长期趋势做出预测。特别要指出模型中的优点与不足之处。二 模型的假设1、不考虑国际迁入迁出对人口增长的影响。2、假设题目给出的生育率已经考虑了一胎双生和多生的因素。3、假设在预测期内国家在控制人口方面的政策、人们的物质生活水平及生育意识等国情不会有太大的变化。4、假设预测期内递推各年人口所用到的城镇乡男女各年龄层的死亡率、城镇乡女性各生育年龄层的生育率、城镇乡的性别出生比、城镇乡年迁入迁出率满足相应的由01至05年的数据拟合的函数曲线。5、再考虑城镇乡的迁入迁出问题时（以城市为例），假设城市每年迁入的人口中各年龄人口的比例与该年城市总人口各年龄人口的比例一致。三 符号说明符号示意备注第i岁0，90第j岁15，49第k年2001，第k年年初全国总人口数以万为单位、第k年年初城市、镇、乡总人口数以万为单位、第k年第i岁城市、镇、乡男性数占当年城市、镇、乡总人口数的比例%、第k年第i岁城市、镇、乡女性数占当年城市、镇、乡总人口数的比例%、第k年第i岁城市、镇、乡男性死亡数占当年城市、镇、乡男性数的比例20012005、第k年第i岁城市、镇、乡女性死亡数占当年城市、镇、乡女性数的比例20012005、第k年第j岁城市、镇、乡女性的生育率%20012005、第k年城市、镇、乡男女出生比例女以100计、第k年的城市、镇、乡迁入迁出率%正值表示迁入20012005第k年中国人口1%调查数据中城市、镇、乡男性的人口数以个为单位第k年中国人口1%调查数据中城市、镇、乡女性的人口数以个为单位拟合的城市、镇、乡男性的平均死亡率拟合的城市、镇、乡女性的平均死亡率拟合的城市、镇、乡女性的平均生育率城市、镇、乡出生性别比平均值城市、镇、乡迁入迁出率平均值四 问题的分析及模型的建立4.1模型的引出英国神父Malthus于1798年提出：人口增长数与原人口数成正比，设定为单位时间段内人口相对增长率r（r为出生率与死亡率之差），从而得到人口增长的Malthus模型。荷兰科学家Verhulst在1840年提出：由于生存资源的有限性，物种成员之间的竞争与约束必定影响他们的数量增长，于是加入约束项，从而得到人口增长的Logistic模型。现当代许多的模型几乎都是在Logistic模型的基础之上发展而来。总结这一类模型，发现它们都是从人口的总数变化着手，试图找到这样一条曲线，使得近些年的人口数尽量的落在这条曲线上，并以之预测未来的人口数变化。然而，人口增长不简单的为一个自然增长，它的变化依赖于各个国家、地区的社会、经济、文化等诸多因素。结合中国具体国情，譬如年龄结构老龄化，男女出生比例不平衡、乡镇人口城市化等，都直接或者间接影响当代或者更远的中国人口的变化。由Logistic模型发展而来的各种模型往往不能深入和独立地反映这些因素的影响，更加不能体现这些因素随国情变化，受国策调控的性质。在Logistic路途上找不到终点的我们，决定回到原点，从人口增长的基本关系出发，即：上年人口数+出生数-死亡数迁入迁出数=下年人口数。我们将全国人口总数按行政级别分为城、镇、乡，以便于考虑城镇乡迁入迁出因素；按年龄分为0岁到90岁，以便于考虑年龄结构的老龄化因素；按性别分开，以便于考虑出生性别比因素。4.2模型的分析我们以城镇乡为主来划分，考虑到三者分析过程极为相似，下面仅以城市人口增长过程为例进行递推分析。4.2.1基本算式第k年年初城市总人口数的计算（k=20012005）分别表示第k年中国人口1%调查数据中城市男性和女性的人口数；分别表示第k年中国人口1%调查数据中镇男性和女性的人口数；分别表示第k年中国人口1%调查数据中乡男性和女性的人口数；于是得到城市人口总数的初值计算式：-即：第k年年初城市总人口数第k年年初全国总人口数第k年1%人口调查数据中城市男女总数占调查总数的比例城市人口总数的递推表示第k年城市0岁人口成活数，即净出生人口数；表示第k年城市人口死亡总数；表示第k年城市人口迁入迁出总数；于是得到如下城市人口总数的递推式：-即：第k+1年年初城市人口总数=第k年年初城市人口总数+第k年净出生人口数+第k年人口死亡总数+第k年迁入迁出总数城市各年龄男、女性人口数的递推表示第k+1年年初第i+1岁城市男性人口数；表示第k年年初第i岁城市男性人口数；表示第k年第i岁城市男性人口死亡数；表示第k年第i岁城市迁入迁出的男性人口数；于是得到同一年龄的城市男性人口数的递推式：-即：第k+1年年初第i+1岁城市男性人口数=第k年年初第i岁城市男性人口数+第k年第i岁城市人口死亡数+第k年第i岁城市迁入迁出的人口数同理：可得同一年龄的城市女性的递推式：-城市年出生人口及男女出生比例的计算表示第k年城市出生人口数，包括死亡的0岁人口数；表示第k年城市女性生育人口数；于是得到城市年出生人口恒等式：-即：第k年城市出生人口数=第k年城市女性生育人口数且有男女出生比例关系式：-全国总人口数的计算对镇、乡人口可以由类似于到的步骤得到相应关于镇、乡人口变化的关系式。从而得第k年的全国总人口数的计算式：-即：第k年年初全国总人口数=第k年年初城市总人口数+第k年年初镇总人口数+第k年年初乡总人口数4.2.2递推性分析及参数的选用我们先作如下讨论：对于第年，我们已知初值；由式可计算出第+1年年初城市总人口数；再由式可计算出第k+1年第i+1岁城市男性数占当年城市总人口数的比例；式可计算出第k+1年第i+1岁城市女性数占当年城市总人口数的比例；如果已知和，加上算出的值，联立式和式可算出；如果再已知，我们就可以类推第k+2年人口及比例了。综上所述：我们设计出基本参数：。初值情况：题目给出的每年人口统计数据包括城镇乡男女不同年龄的人口比例和死亡率，生育期女性的生育率，城镇乡人口总数（须由给出的中国人口1%调查数据进行换算），男女出生比例等。查找人口统计资料易得2001年至2005年的全国人口总数（以万为单位）分别为：127627、128453、129227、129988、130756。所以：由式可分别求出为；30885、33607、33621、33566、36241。再由式可分别求出为0.0925、0.0044、0.0027、0.0834；以2001年作为推算起始年，则初值为。用同样的方法对镇、乡也可以得到相应的初值。使用matlab编程求初值，程序见附录：yuce.m4.2.3参数估计死亡率年龄 函数的建立首先以近五年城市男性各年龄层的死亡率为例，做同一年龄的城市男性的死亡率变化分析。为排除偶然性误差的影响，我们需要对每一年城市男性各年龄层的死亡率进行指数拟和，得到修正的城市男性各年龄层的死亡率，再观察同一年龄的城市男性近五年来的死亡率变化情况。利用matlab的拟合工具箱，求得各散点图的拟合曲线。横坐标为年龄（岁） 纵坐标为死亡率（）对比以上五条曲线，不难看出：这五年来同一年龄的城市男性的死亡率波动较小。事实上，因为国家在控制人口方面的政策和人们物质文化生活水平的增幅在短期内不会有太大的变化，我们有理由相信各年龄层死亡率分布在中短期内不会有太大的变化。因此，作中短期预测，我们将各年龄层死亡率视为一个不变的常数。以下给出由2001年到2005年的城镇乡男女各年龄层的死亡率数据做整体拟和得到按年龄分布的的平均死亡率曲线，并认为以后各年份城镇乡男女各年龄层的死亡率就在这些曲线上求得。下面利用matlab的拟合工具箱，求得各散点图的拟合曲线。再对拟合方程进行编程，程序见附录：swa.m和swv.m城市男性的平均死亡率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）城市女性的平均死亡率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）镇男性的平均死亡率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）镇女性的平均死亡率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）乡男性的平均死亡率曲线及拟和函数： 横坐标为年龄（岁） 纵坐标为死亡率（）乡女性的平均死亡率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）考虑到岁人口死亡率的特殊性，有必要将之独立出来，为了简化讨论，我们认为其死亡率为不随年份变化的定值。以下通过算术平均值来确定这个定值：单位为城市0岁男性的平均死亡率为：=6.116城市0岁女性的平均死亡率为：=6.344镇0岁男性的平均死亡率为：=6.4镇0岁女性的平均死亡率为：=10.99乡0岁男性的平均死亡率为：=19.31乡0岁女性的平均死亡率为：=23.96此项求解程序见附录：yuce.m生育率年龄 函数的建立：类比死亡率的分析，在国家计划生育国策发生大的变化之前我们认定同一生育年龄女性的生育率为一定值。以下给出由2001年到2005年的城镇乡女性各生育年龄层的生育率数据拟和得到的按年龄分布的平均生育率曲线，并认为以后各年份城镇乡女性各生育年龄层的生育率就在这些曲线上求得。下面利用matlab的拟合工具箱，求得各散点图的拟合曲线。再对拟合方程进行编程，程序见附录：spl.m城女的平均生育率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为死亡率（）镇女性的平均生育率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为生育率（）乡女性的平均生育率曲线及拟和函数：横坐标为年龄（岁） 纵坐标为生育率（）性别出生比：由1994年到2005年城镇乡的性别出生比及中国的具体国情，我们仍然有理由认定未来较长期限内性别出生比将维持为一个常数。以下通过1994年到2005年城镇乡的性别出生比的算术平均值来确定：城市出生性别比平均值：=111.829镇出生性别比平均值：=117.9乡出生性别比平均值：=119.65此项求解程序见附录：yuce.m迁入迁出率：由2001年到2004年城镇乡的迁入迁出率及中国的具体国情，我们仍然有理由认定未来较长期限内迁入迁入率将维持为一个常数。以下通过2001年到2004年城镇乡的迁入迁出率的算术平均值来确定：城市迁入迁出率平均值：=0.0458；镇迁入迁出率平均值：=0.0866；乡迁入迁出率平均值：=-0.0202；此项求解程序见附录：yuce.m和qyll.m4.3模型的建立以城市人口数的递推模型为例：利用matlab编程，程序见附录：dt.m、yuce.m和yucetu.m4.4模型的结果利用本模型以2001年为初算值，计算2002年到2060年的全国总人口数，主要数据如下：年份200220032004200520062007200820092010总人口数/万128770 129440130150130890131670132490133330134200135100年份201520202025203020352040204520502055总人口数/万139730144000147240148960148720146140141540136050130060预测人口变化趋势图如下：横坐标为年份（年） 纵坐标为全国总人口数（万人）由上图显示结果，我们有理由推断中国未来人口的发展趋势是有可能下降的。读图可知，在人口增长达到顶峰后会有一个下降的过程。五 模型的检验据调查，2002年到2005年的全国总人口数，下将实际值与预测值作比较： 年份2002200320042005实际值/万128453129227129988130756预测值/万128770129440130150130890相对误差/%0.2470.1650.1250.102对比可知：本模型用于预测短期全国过总人口数的预测具有较高置信度。相关文献表明，中国未来人口增长将在本世纪三四十年代达到顶峰；国家人口发展战略研究报告指出，该峰值将控制在15亿左右。本模型预测的结果是2032年总人口数达到峰值为149130万。与相关预测较为接近，具有良好的被接受性。六 模型的优缺点分析优点：1、Malthus和Logistic两个经典的人口模型都是从人口的总数变化着手。而本模型将全国总人口按行政级别、年龄、性别分类，先对各小类地变化进行分析计算，再将各类人口求和便得到全国总人口的变化情况。2、本模型是从人口变化最基本最简单的关系(人口数+出生数-死亡数迁入迁出数=下年人口数)为基础，显得更为合理。3、本模型将全国总人口按行政级别、年龄、性别分类，各类的变化情况就分别反映了乡村人口城镇化、老龄化进程加速、出生人口性别比持续升高等中国人口发展特点对人口增长的影响。4、本模型的参数都是可以据具体情况而修正的。城镇乡男女各年龄层的死亡率、城镇乡女性各生育年龄层的生育率、城镇乡的性别出生比、城镇乡年迁入迁出率等参数受人口方面的国情的影响。根据模型预测的人口数来判断当人口变化到种状况时国家将哪些宏观调控，从而对上述参数做适当修正。缺点 ：1、2001到2005年五年的数据太少，去估计以后的参数不太精确。2、本模型是一个递推式的模型，当预测期增长时易产生较大的累积误差。七 模型的推广、改进本模型具有十分良好的可推广性。其一在参数估计方法的多样性，较好地参数估计方法可以大大提高模型的精确度。其二在模型自身的改进。可以根据要求预测人口的详细程度适当增加或减少信息量，以达到因时而已、因地制宜的效果。其三在于参数的自变性。国情、政策的变化随时影响参数的变化。最好是能够建立城镇乡男女各年龄层的死亡率、城镇乡女性各生育年龄层的生育率、城镇乡的性别出生比、城镇乡年迁入迁出率等参数关于人口总数的函数关系式。这样，每算出一年的人口数就能算出相应的参数，用于该年人口向下一年人口的递推。这样做能够避免因参数估计不准带来的累积误差。八 参考文献1 江蕾、蒋远馨，中国人口死亡率与经济发展水平关系的实证研究:1952-2002年，中国人口科学，S1期：P1-5，2005。2 刘琼荪、龚劬，数学实验，北京：高等教育出版社，2004。3 乐经良，数学实验，北京：高等教育出版社，1999。4 中华人民共和国国家统计局，中国统计年鉴2005，/cgi-bin/bigate.cgi/b/g/g//tjsj/ndsj/2005/indexch.htm，2007-9-23。九 附录swa.mfunction swa=swa(i,czx)if i=1 if czx=c swa=0.0061;%swac0 elseif czx=z swa=0.0064;%swaz0 elseif czx=x swa=0.0193;%swax0 endelse if czx=c swa=0.03245*exp(0.09675*i)/1000; elseif czx=z swa=0.0428*exp(0.0957*i)/1000; elseif czx=x swa=0.06133*exp(0.09297*i)/1000; endendswv.mfunction swv=swv(i,czx)if i=1 if czx=c swv=0.0063;%swvc0 elseif czx=z swv=0.0110;%swvz0 elseif czx=x swv=0.0239;%swvx0 endelse if czx=c swv=0.00628*exp(0.1124*i)/1000; elseif czx=z swv=0.02791*exp(0.09562*i)/1000; elseif czx=x swv=0.03909*exp(0.09424*i)/1000; endendsyl.mfunction syl=syl(j,czx)if czx=c syl=(0*exp(-(j-26.58)/0.02614)2) + 101.9*exp(-(j-25.08)/3.793)2) + 29.27*exp(-(j-29.65)/5.76)2)/1000;elseif czx=z syl=(43.26*exp(-(j-28.74)/6.011)2) + 119.5*exp(-(j-24.32)/3.272)2) + 34.95*exp(-(j-21.47)/1.505)2)/1000;elseif czx=x syl=(145.5*exp(-(j-22.96)/2.91)2) + 85.5*exp(-(j-27.84)/5.863)2)/1000;endqyll.mfunction qy=qyll(g,czx)if czx=c qy2001=0.028;qyoo=0.01;elseif czx=z qy2001=0.04;qyoo=0.01;elseif czx=x qy2001=-0.01;qyoo=-0.001;endbg=40;%kg=(qyoo-qy2001)/bg;if g=2001+bg qy=(bg+2001)*qy2001-2001*qyoo+(qyoo-qy2001)*g)/bg;else qy=qyoo;enddt.mfunction tz2=dt(tz,czx,qyl,csb)x=tz(1,1);na=tz(:,2);nv=tz(:,3); Esw=0; for i=1:90 Esw=Esw+na(i)*swa(i,czx)+nv(i)*swv(i,czx); end %Esw=Esw x2=x*(1+na(1)*(1-swa(1,czx)+nv(1)*(1-swv(1,czx)-Esw+qyl); %x2=x2 for i=1:90 na2(i+1)=na(i)*(1-swa(i,czx)+qyl)*x/x2; nv2(i+1)=nv(i)*(1-swv(i,czx)+qyl)*x/x2; end Esy=0; for j=15:49 Esy=Esy+x2*nv2(j)*syl(j,czx); end %Esy=Esy na2(1)=csb*Esy/(csb+100)*x2); nv2(1)=100*Esy/(csb+100)*x2); %na2(1),nv2(1)tz2(1,1)=x2;for h=1:90 tz2(h,2)=na2(h); tz2(h,3)=nv2(h);end%na2=na2(1:91),nv2=nv2(1:91)yuce.mfunction RenKou=yuce(N)tzc=zeros(92,3);tzz=zeros(92,3);tzx=zeros(92,3);tz0=load(tz0.txt);tzc(1,1)=tz0(1,1);tzz(1,1)=tz0(1,4);tzx(1,1)=tz0(1,7);for h=1:91 tzc(h,2)=tz0(h,2);tzz(h,2)=tz0(h,5);tzx(h,2)=tz0(h,8); tzc(h,3)=tz0(h,3);tzz(h,3)=tz0(h,6);tzx(:,3)=tz0(h,9);endif N2001 RenKouShu=Error;elseif N=2001 RenKouShu=tzc(1,1)+tzz(1,1)+tzx(1,1); RenKou(1)=RenKouShu; %tzc(1,1),tzz(1,1),tzx(1,1)else%-%X2001=load(tz2001.txt);X2001c=X2001(7,14)*(X2001(1,14)+X2001(2,14)/(X2001(1,14)+X2001(2,14)+X2001(3,14)+X2001(4,14)+X2001(5,14)+X2001(6,14);X2001z=X2001(7,14)*(X2001(3,14)+X2001(4,14)/(X2001(1,14)+X2001(2,14)+X2001(3,14)+X2001(4,14)+X2001(5,14)+X2001(6,14);X2001x=X2001(7,14)*(X2001(5,14)+X2001(6,14)/(X2001(1,14)+X2001(2,14)+X2001(3,14)+X2001(4,14)+X2001(5,14)+X2001(6,14);X2002=load(tz2002.txt);X2002c=X2002(7,14)*(X2002(1,14)+X2002(2,14)/(X2002(1,14)+X2002(2,14)+X2002(3,14)+X2002(4,14)+X2002(5,14)+X2002(6,14);X2002z=X2002(7,14)*(X2002(3,14)+X2002(4,14)/(X2002(1,14)+X2002(2,14)+X2002(3,14)+X2002(4,14)+X2002(5,14)+X2002(6,14);X2002x=X2002(7,14)*(X2002(5,14)+X2002(6,14)/(X2002(1,14)+X2002(2,14)+X2002(3,14)+X2002(4,14)+X2002(5,14)+X2002(6,14);X2003=load(tz2003.txt);X2003c=X2003(7,14)*(X2003(1,14)+X2003(2,14)/(X2003(1,14)+X2003(2,14)+X2003(3,14)+X2003(4,14)+X2003(5,14)+X2003(6,14);X2003z=X2003(7,14)*(X2003(3,14)+X2003(4,14)/(X2003(1,14)+X2003(2,14)+X2003(3,14)+X2003(4,14)+X2003(5,14)+X2003(6,14);X2003x=X2003(7,14)*(X2003(5,14)+X2003(6,14)/(X2003(1,14)+X2003(2,14)+X2003(3,14)+X2003(4,14)+X2003(5,14)+X2003(6,14);X2004=load(tz2004.txt);X2004c=X2004(7,14)*(X2004(1,14)+X2004(2,14)/(X2004(1,14)+X2004(2,14)+X2004(3,14)+X2004(4,14)+X2004(5,14)+X2004(6,14);X2004z=X2004(7,14)*(X2004(3,14)+X2004(4,14)/(X2004(1,14)+X2004(2,14)+X2004(3,14)+X2004(4,14)+X2004(5,14)+X2004(6,14);X2004x=X2004(7,14)*(X2004(5,14)+X2004(6,14)/(X2004(1,14)+X2004(2,14)+X2004(3,14)+X2004(4,14)+X2004(5,14)+X2004(6,14);X2005=load(tz2005.txt);X2005c=X2005(7,14)*(X2005(1,14)+X2005(2,14)/(X2005(1,14)+X2005(2,14)+X2005(3,14)+X2005(4,14)+X2005(5,14)+X2005(6,14);X2005z=X2005(7,14)*(X2005(3,14)+X2005(4,14)/(X2005(1,14)+X2005(2,14)+X2005(3,14)+X2005(4,14)+X2005(5,14)+X2005(6,14);X2005x=X2005(7,14)*(X2005(5,14)+X2005(6,14)/(X2005(1,14)+X2005(2,14)+X2005(3,14)+X2005(4,14)+X2005(5,14)+X2005(6,14);Esw2001c=0;Esw2001z=0;Esw2001x=0;for i=1:91 Esw2001c=Esw2001c+X2001(i,2)*X2001(i,3)+X2001(i,4)*X2001(i,5); Esw2001z=Esw2001z+X2001(i,6)*X2001(i,7)+X2001(i,8)*X2001(i,9); Esw2001x=Esw2001x+X2001(i,10)*X2001(i,11)+X2001(i,12)*X2001(i,13);endqyl2001c=X2002c/X2001c-1-X2001(1,2)*(1-X2001(1,3)+X2001(1,4)*(1-X2001(1,5)+Esw2001c;qyl2001z=X2002z/X2001z-1-X2001(1,6)*(1-X2001(1,7)+X2001(1,8)*(1-X2001(1,9)+Esw2001z;qyl2001x=X2002x/X2001x-1-X2001(1,10)*(1-X2001(1,11)+X2001(1,12)*(1-X2001(1,13)+Esw2001x;Esw2002c=0;Esw2002z=0;Esw2002x=0;for i=1:90 Esw2002c=Esw2002c+X2002(i,2)*X2002(i,3)+X2002(i,4)*X2002(i,5); Esw2002z=Esw2002z+X2002(i,6)*X2002(i,7)+X2002(i,8)*X2002(i,9); Esw2002x=Esw2002x+X2002(i,10)*X2002(i,11)+X2002(i,12)*X2002(i,13);endqyl2002c=X2003c/X2002c-1-X2002(1,2)*(1-X2002(1,3)+X2002(1,4)*(1-X2002(1,5)+Esw2002c;qyl2002z=X2003z/X2002z-1-X2002(1,6)*(1-X2002(1,7)+X2002(1,8)*(1-X2002(1,9)+Esw2002z;qyl2002x=X2003x/X2002x-1-X2002(1,10)*(1-X2002(1,11)+X2002(1,12)*(1-X2002(1,13)+Esw2002x;Esw2003c=0;Esw2003z=0;Esw2003x=0;for i=1:90 Esw2003c=Esw2003c+X2003(i,2)*X2003(i,3)+X2003(i,4)*X2003(i,5); Esw2003z=Esw2003z+X2003(i,6)*X2003(i,7)+X2003(i,8)*X2003(i,9); Esw2003x=Esw2003x+X2003(i,10)*X2003(i,11)+X2003(i,12)*X2003(i,13);endqyl2003c=X2004c/X2003c-1-X2003(1,2)*(1-X2003(1,3)+X2003(1,4)*(1-X2003(1,5)+Esw2003c;qyl2003z=X2004z/X2003z-1-X2003(1,6)*(1-X2003(1,7)+X2003(1,8)*(1-X2003(1,9)+Esw2003z;qyl2003x=X2004x/X2003x-1-X2003(1,10)*(1-X2003(1,11)+X2003(1,12)*(1-X2003(1,13)+Esw2003x;Esw2004c=0;Esw2004z=0;Esw2004x=0;for i=1:90 Esw2004c=Esw2004c+X2004(i,2)*X2004(i,3)+X2004(i,4)*X2004(i,5); Esw2004z=Esw2004z+X2004(i,6)*X2004(i,7)+X2004(i,8)*X2004(i,9); Esw2004x=Esw2004x+X2004(i,10)*X2004(i,11)+X2004(i,12)*X2004(i,13);endqyl2004c=X2005c/X2004c-1-X2004(1,2)*(1-X2004(1,3)+X2004(1,4)*(1-X2004(1,5)+Esw2004c;qyl2004z=X2005z/X2004z-1-X2004(1,6)*(1-X2004(1,7)+X2004(1,8)*(1-X2004(1,9)+Esw2004z;qyl2004x=X2005x/X2004x-1-X2004(1,10)*(1-X2004(1,11)+X2004(1,12)*(1-X2004(1,13)+Esw2004x; qylc=(qyl2001c+qyl2002c+qyl2003c+qyl2004c)/4;qylz=(qyl2001z+qyl2002z+qyl2003z+qyl2004z)/4;qylx=(qyl2001x+qyl2002x+qyl2003x+qyl2004x)/4;swac0=(X2001(1,3)+X2002(1,3)+X2003(1,3)+X2004(1,3)+X2005(1,3)/5;swvc0=(X2001(1,5)+X2002(1,5)+X2003(1,5)+X2004(1,5)+X2005(1,5)/5;swaz0=(X2001(1,7)+X2002(1,7)+X2003(1,7)+X2004(1,7)+X2005(1,7)/5;swvz0=(X2001(1,9)+X2002(1,9)+X2003(1,9)+X2004(1,9)+X2005(1,9)/5;swax0=(X2001(1,11)+X2002(1,11)+X2003(1,11)+X2004(1,11)+X2005(1,11)/5;swvx0=(X2001(1,13)+X2002(1,13)+X2003(1,13)+X2004(1,13)+X2005(1,13)/5; csbl=load(chushengbi.txt);csbc=mean(csbl(:,1);csbz=mean(csbl(:,2);csbx=mean(csbl(:,3);%swac0,swvc0,swaz0,swvz0,swax0,swvx0,csbc,csbz,csbx,qylc,qylz,qylx%qyl2001c,qyl2002c,qyl2003c,qyl2004c%qyl2001z,qyl2002z,qyl2003z,qyl2004z%qyl2001x,qyl2002x,qyl2003x,qyl2004x%-% czx=Null;RenKou=zeros(1,N-2000);%RenKou RenKou(1)=tzc(1,1)+tzz(1,1)+tzx(1,1); for k=2001:N-1 czx=c;qyl=qyll(k,czx);csb=csbc;tzc=dt(tzc,czx,qyl,csb); czx=z;qyl=qyll(k,czx);csb=csbz;tzz=dt(tzz,czx,qyl,csb); czx=x;qyl=qyll(k,czx);csb=csbx;tzx=dt(tzx,czx,qyl,csb); RenKou(k-1999)=tzc(1,1)+tzz(1,1)+tzx(1,1); end %tzc(1,1),tzz(1,1),tzx(1,1) RenKouShu=tzc(1,1)+tzz(1,1)+tzx(1,1);end%RenKouyucetu.mnian=2001:2055;ren=yuce(2055);plot(nian,ren,o);chushengbi.txt114.52124.33116.15111.92115.62117.75111.68111.68117.7108.81125.9118.85110.68108.73119.98110.27118.4122.03113 116.3119.3109.28116.02117.59111.37123.12122.11112.06110.97120.9114.44126.9122.21113.92117.21121.21tz0.txt308850.00470.0042165480.00590.0052801940.00670.005800.00440.00400.00580.004600.00580.004800.00420.003900.00520.004300.00630.004900.00480.004200.00610.005400.00690.005900.00490.004400.00650.005500.00730.006200.00510.004400.00720.005700.00770.006800.00530.004900.00720.006200.00840.007400.00510.00500.00720.006200.00830.007200.00550.00500.00730.006400.00910.008200.00570.005200.00770.006600.00970.008500.00610.005700.00860.007500.01070.009600.0070.006500.00990.008600.01270.011300.00720.006600.00990.008300.01210.010700.00730.00700.00950.008700.01120.01