1、基于常微分方程的中国人口增长预测摘要人口问题是制约我国经济和社会发展的一个很重要的方面，因此对人口的预测是一个由来已久的研究课题，但传统方法往往是使用单一的模型进行计算，因而考虑的影响因素不够全面。事实上，人口的增长的影响因素是多方面的，如出生率、死亡率、城市化进程、国家生育政策、经济发达程度、生活方式、传统习惯和自然灾害等。结合我国的发展实际，考虑其中主要的影响因素建立了我们的模型。本文是基于微分方程的思想建立的人口预测模型，以附录 2 的人口普查数据为基础数据。在对数据进行评估和修订的基础上，建立女性的生育儿童的函数，利用各个年龄段的生育率求出婴儿的出生率。之后，考虑各个年龄段的人口分布，
2、利用人口年龄密度函数、对应年龄人的死亡率和新生儿的出生率来算次年的人口总数。建立了人口总数、各个年龄段人口密度、自然增长率等关于年龄与时间的数学函数模型。在此基础上，考虑到国家政策调控、农村人口城市化、社会老龄化、科技医疗水平进步以及可能的如艾滋病等疾病对出生率、死亡率等的影响，并采用添加干预因子的方法对模型进行了矫正分析，以体现这些因素对人口增长趋势的影响。采取上述人口预测的建模思想，我们最终得到了中国人口增长的数学模型，利用附录 2 中人口数据进行模拟、预测计算，得到了我国人口增长预测的一些结果：对我国人口进行中短期和长期预测，我国人口在 XX 年达到 13.565 亿；在 2020 年达
3、到 14.767 亿；在 2040 年达到 14.843 亿。其中我国在 2028 年左右会达到人口最高蜂，约为 15.17 亿人左右。在对模型的分析中，我们看到我国将出现老龄化趋势。在 XX 年后，老龄化趋势加快，预计在 XX 年达到 12%，此后我国将进入人口老龄化迅速发展时期。预测显示从 XX 到2040 年间里，中国老年人比例会增加一倍左右。关键词：干预因子，人口年龄分布，微分方程，老龄化，人口结构，出生率一、问题的重述我国是一个人口大国，人口问题对我国的经济发展有着很大的影响，是制约我国经济发展的关键因素之一。根据最新的“XX 年全国%1 人口抽样调查公报”指出，我国人口和计划生育工
4、作进入稳定低生育水平、统筹解决人口问题、促进人的全面发展的新阶段。要统筹解决人口与社会发展的各种问题，必须对目前和将来一段时间内的人口问题进行科学的判断，如城市化问题、老龄化问题、出生人口性别比持续升高问题。因此要制度科学的人口政策需要对未来中国的人口变化进行科学预测。关于中国人口的预测，目前已经有了许多的模型。根据题意，本文从中国的实际情况和人口增长的特点，建立对中国人口增长的中短期和长期的模型预测。二、问题的分析由于影响人口增长的因素比较复杂，从中国的实际情况和人口增长的特点如果一开始把众多因素都考虑进去，将会对模型的求解很困难。所以我们首先考虑影响人口变化的一些基本因素，如出生率、死亡率
5、、人口年龄的结构密度等。其中人口的年龄结构将在未来一段时间内，产生很大的影响。在人口预测中人口按年龄分布状况是十分重要的，因为不同的年龄人的生育率和死亡率有很大的差别，未来的人口发展状况将大不一样。下图是我们根据附表 2的数据计算出的 01-05 年人口年龄密度图。我们从图 2-1 中可看到 01-05 年年龄结构中，老年人数量有增大的趋势。对此我们可预见我国人口结构有老龄化趋势。在模型的建立中我们可看到这一点。图 2-1在解决变化率的问题时，常常引用微分方程的方法1。在此基础上，我们考虑人口按年龄的分布，即除了时间变量外，年龄是另一自变量，来建立微分方程模型。考虑各种微分方程的方法，指数增长
6、模型和阻滞增长模型(Logistic 人口增长模型)显然考虑的因素太少，很难做出合适的模型。宋健模型中考虑的问题比较全面，可以取其优秀的地方学习。可是宋健模型假定死亡率不随时间 t 的变化而改变，未考虑各种因素对死亡率的影响。由此，我们考虑到在宋健模型的基础上加以深化，考虑人口结构、城市化、老龄化、社会发展等因素对死亡率、生育率的影响，用变化的死亡率和生育率来预计未来的人口结构等。从而实现中国人口增长的中短期和长期的模型预测。三、模型的基本假设和符号说明模型的假设1 我们研究的人口是一个整体系统2 影响人口的变化因素只有时间、婴儿的出生、人口的死亡、人口的迁移3 期间不考虑突发事件对人口变化的
7、影响4 系统是封闭的，不考虑出国与海归人士对我国人口的影响模型符号的说明和名词解释N(t) t 时刻该地区人口的总数rm 人的最高寿命r 表示年龄t 表示时间F(r,t) 表示在时刻 t 该地区一切年龄小于 r 的人数p(r,t) 表示在时刻 t，年龄为 r 的人数 ，M(r,t) 表示在时刻 t 该地区年龄 r 为的死亡数K(r,t) 表示女性在人口中占的比例b(r,t) 表示 t 时刻平均每个 r 岁的女性的生育数育龄区表示时刻 t 单位时间内平均每个育龄女性的生育数h 表示一个女性在 r 岁时的生育概率u 表示在时刻 t 年龄 r 的人的死亡率四、模型的建立与求解1 模型的建立我们记 u
8、(r,t)为时刻 t 年龄 r 的人的死亡率，其含义是，updr 表示时刻 t 年龄在r，r+dr内单位时间死亡的人数。为了得到 p 满足的方程，考察时刻 t 年龄在r,r+dr内的人到时刻 t+dt 的情况。他们中活着的那一部分人的年龄变为r+dr1，r+dr+dr1，这里 dr1=dt。而在 dt 这段时间内死亡的人数为 up(r,t)drdt.于是上式可写为：两边同除以 drdt，于是：=-u(r,t)p(r,t)取极限：=-up， 这即是人口密度函数 p 的一阶偏微分方程。对方程有两个定解条件：1 初始密度函数记为 p2 单位时间内出生的婴儿数记为 p(0,t)=f(t)(2)出生儿童
9、的模型在育龄区间r1，r2内为 t 时刻平均每个 r 岁的女性的生育数则在 t 时刻，出生儿童的总数为：方程就是我们建立人口预测模型的基础。2 模型的求解当 u 不依赖于 t，即在社会安定的局面下和不太长的时间内，死亡率大致与时间无关，于是对方程求得：由方程我们利用已知的 01-05 年的数据对其交替求解设 t=1，2，3，4，5 令 XX 年为 t=1，此后依次类推则根据附表 2 中拟和出年龄密度分布p(r,1),p(r,2),p(r,3),P(r,4),p(r,5)因此，我们可以依次推出t=6，t=7,t=8 ， 依次计算，就可以预测以后的人口总数和年龄结构。考虑到其他一些因素的影响，在具
10、体的计算过程中加入干预因子优化模型。例如，考虑到城市化过程对人口分布的影响，引入分布因子的概念，根据目前的城市化进程和国家制定的城市化政策规划，我们设定了如下变化的因子：io1=0.002*t+0.25;io2=0.003*t+0.15;io3=-0.007*t+0.6;图 4-1老龄化对死亡率的影响:lp=0.000019212*t;科学文化发展对死亡率降低到影响:sw=0.000034523*t;具体的见附录3、 模型结论的分析通过已知的数据和 matlab 编程2，通过建立的模型我们得出我国未来的人口总数变化曲线，1 、人口规模与增长：图 4-1 即是我国未来一段时间内，人口的变化趋势。
11、国家统计局根据 XX 年全国人口 1%抽样调查结果推算 XX 年底人口为 130756 万人。根据我们的模型预测结果，我国总人口在 XX 年将达到 13.565 亿人。预计到 2028 年将达到峰值，约为 15.17 亿，此后开始缓慢下降。如果长期保持 1.8的生育水平，则到 2040 年达到 14.843 亿。图 4-2 我国未来 50 年人口变化曲线其预测的中国人口发展趋势为：2、人口年龄结构的变化：由于生育率的快速下降和长期的低生育水平，我国未来年龄结构变化将体现出两大特征3：显著的的人口红利和人口老龄化。人口红利是指人口转变过程中出现的劳动年龄人口比例上升带来的经济增长效应。它的时间长
12、度与程度取决于死亡率、生育率的下降的速度和时间差。他将在未来对我国的就业带来很大的压力和挑战。我国未来发展的另一大挑战是持续、快速的人口老龄化。无论是老年人规模还是其比例都在迅速的上升。下面的图我们可以看出。人口的老龄化将导致抚养比不断的提高，对社会的保障体系和公共服务体系的压力加大，并影响到社会代际关系的和谐。同时由于城镇化的进程加快，进城打工人数的增加，将使得农村的老龄化形式更为严峻。五、模型的评价1、模型的优点：本模型在建立的过程中，考虑了多个因素比以往的人口预测模型如指数增长模型、阻滞增长模型有了很大的改进。由于我们是根据人口密度推总人口分布的，人口与时间和年龄都有联系，所以我们可以推
13、出将来的人口的年龄结构。此外，我们的模型还涉及到女性性别比、总和生育率，这对我国的计划生育政策提供依据和手段。对今后一段时间内的预测很好的吻合了我国最新的国家人口战略发展研究报告。2、模型存在的问题：从控制论的观点看，在方程描述的人口系统中 p(r,t)可视为状态变量，p(0,t)=f(t)视为控制变量，是分布参数系统的边界控制函数。式表明控制输入中含有状态变量，形成反馈，视为反馈增益。并且通常是正反馈，即人口密度函数 p(r,t)增加，通过婴儿出生率 f(t)又使 p(r,t)进一步增大。方程式中因子 f(t-r)表明这种反馈有很大的滞后作用，所以一旦政策失误，使 p(r,t)在一段时间内增
14、长得过多过快，再想同过控制手段和 h(r,t)把人口增长降下来，就很困难，可能需要几代人。显然，随着科学的发展，统计预测的加强，政策方面的大失误，我们暂不考虑。参考文献1姜启源 谢金星等，数学模型，高等教育出版，XX2苏金明 阮沈勇 matlab 实用教程，XX3陈伟 中国未来人口发展趋势：XX2050 人口研究 XX.74刘祯 张秋辉 我国人口发展现状研究 科技统计 河南科技 XX.5源程序附录%pp.mload a01-05.matfor r=1:91p1(r,1)=(aXX(r,2)+aXX(r,4);p2(r,1)=(aXX(r,6)+aXX(r,8);p3(r,1)=(aXX(r,1
15、0)+aXX(r,12);p(r,1)=(p1(r,1)*(z01_05(2,1)+z01_05(3,1)+p2(r,1)*(z01_05(4,1)+z01_05(5,1)+p3(r,1)*(z01_05(6,1)+z01_05(7,1)/(z01_05(2,1)+z01_05(3,1)+z01_05(4,1)+z01_05(5,1)+z01_05(6,1)+z01_05(7,1);p1(r,2)=(aXX(r,2)+aXX(r,4);p2(r,2)=(aXX(r,6)+aXX(r,8);p3(r,2)=(aXX(r,10)+aXX(r,12);p(r,2)=(p1(r,2)*(z01_05(
16、2,2)+z01_05(3,2)+p2(r,2)*(z01_05(4,2)+z01_05(5,2)+p3(r,2)*(z01_05(6,2)+z01_05(7,2)/(z01_05(2,2)+z01_05(3,2)+z01_05(4,2)+z01_05(5,2)+z01_05(6,2)+z01_05(7,2);p1(r,3)=(aXX(r,2)+aXX(r,4);p2(r,3)=(aXX(r,6)+aXX(r,8);p3(r,3)=(aXX(r,10)+aXX(r,12);p(r,3)=(p1(r,3)*(z01_05(2,3)+z01_05(3,3)+p2(r,3)*(z01_05(4,3)
17、+z01_05(5,3)+p3(r,3)*(z01_05(6,3)+z01_05(7,3)/(z01_05(2,3)+z01_05(3,3)+z01_05(4,3)+z01_05(5,3)+z01_05(6,3)+z01_05(7,3);p1(r,4)=(aXX(r,2)+aXX(r,4);p2(r,4)=(aXX(r,6)+aXX(r,8);p3(r,4)=(aXX(r,10)+aXX(r,12);p(r,4)=(p1(r,4)*(z01_05(2,4)+z01_05(3,4)+p2(r,4)*(z01_05(4,4)+z01_05(5,4)+p3(r,4)*(z01_05(6,4)+z01
18、_05(7,4)/(z01_05(2,4)+z01_05(3,4)+z01_05(4,4)+z01_05(5,4)+z01_05(6,4)+z01_05(7,4);p1(r,5)=(aXX(r,2)+aXX(r,4);p2(r,5)=(aXX(r,6)+aXX(r,8);p3(r,5)=(aXX(r,10)+aXX(r,12);p(r,5)=(p1(r,5)*(z01_05(2,5)+z01_05(3,5)+p2(r,5)*(z01_05(4,5)+z01_05(5,5)+p3(r,5)*(z01_05(6,5)+z01_05(7,5)/(z01_05(2,5)+z01_05(3,5)+z01
19、_05(4,5)+z01_05(5,5)+z01_05(6,5)+z01_05(7,5);end;krt=0.460012;bata=1.8;swl=sum(aXX(:,2).*aXX(:,3)+sum(aXX(:,4).*aXX(:,5);%t=0:50;%io1=0.002*t+0.25;io2=0.003*t+0.15;io3=-0.007*t+0.6;%lp=0.000019212*t;%sw=0.000034523*t;%for t=0:50;f(t+1)=0;for r=0:90if t4p(r+1,t+1)=p(r-t+1,1)*exp(-jf1);end;end;if trjf
20、2=(sum(aXX(1:(r+1),2).*aXX(1:(r+1),3)+sum(aXX(1:(r+1),4).*aXX(1:(r+1),5)*io1(t+1)+(sum(aXX(1:(r+1),6).*aXX(1:(r+1),7)+sum(aXX(1:(r+1),8).*aXX(1:(r+1),9)*io2(t+1)+(sum(aXX(1:(r+1),10).*aXX(1:(r+1),11)+sum(aXX(1:(r+1),12).*aXX(1:(r+1),13)*io3(t+1)/100000-lp(t+1)+sw(t+1);if t4p(r+1,t+1)=f(t-r+1)*exp(-j
21、f2);end;end;end;for j=15:49f(t+1)=f(t+1)+(aXX(j+1,15)*krt*p(j+1,t+1)/100000;end;f(t+1)=bata*f(t+1);if f(t+1)0.002f(t+1)=0.002;end;zrzzl(t+1)=f(t+1)-swl/100000-lp(t+1)+sw(t+1);end;rk(1)=127627;rk(2)=128453;rk(3)=129227;rk(4)=129988;rk(5)=130756;rk(6)=131448;for i=6:50rk(i+1)=rk(i)*(1+zrzzl(i+1);end;t=0:50;plot(t,rk,*);grid on;hold on;
