基于常微分方程的中国人口增长预测

基于常微分方程的中国人口增长预测
 

摘要

人口问题是制约我国经济和社会发展的一个很重要的方面，因此对人口的预测是一个由来已久的研究课题，但传统方法往往是使用单一的模型（如灰色系统模型、自回归模型和指数模型等）进行计算，因而考虑的影响因素不够全面。事实上，人口的增长的影响因素是多方面的，如出生率、死亡率、城市化进程、国家生育政策、经济发达程度、生活方式、传统习惯和自然灾害等。结合我国的发展实际（中国是一个经济高速发展的发展中国家，人口多、区域经济发展不平衡、人口老龄化加剧和相当长一段时间里实行计划生育的国策等），考虑其中主要的影响因素（如人口结构、人口年龄分布以及老龄化等）建立了我们的模型。

    本文是基于微分方程的思想建立的人口预测模型，以附录2的人口普查数据为基础数据。在对数据进行评估和修订的基础上，建立女性的生育儿童的函数，利用各个年龄段的生育率求出婴儿的出生率。之后，考虑各个年龄段的人口分布，利用人口年龄密度函数、对应年龄人的死亡率和新生儿的出生率来算次年的人口总数。建立了人口总数、各个年龄段人口密度、自然增长率等关于年龄与时间的数学函数模型。在此基础上，考虑到国家政策调控、农村人口城市化、社会老龄化、科技医疗水平进步以及可能的如艾滋病等疾病对出生率、死亡率等的影响，并采用添加干预因子的方法对模型进行了矫正分析，以体现这些因素对人口增长趋势的影响。

    采取上述人口预测的建模思想，我们最终得到了中国人口增长的数学模型，利用附录2中人口数据进行模拟、预测计算，得到了我国人口增长预测的一些结果：对我国人口进行中短期和长期预测，我国人口在2010年达到13.565 亿（短期）；在2020年达到 14.767亿（中期）；在2040年达到14.843亿（长期）。其中我国在2028年左右会达到人口最高蜂，约为15.17亿人左右。

在对模型的分析中，我们看到我国将出现老龄化（年龄大于60）趋势。在2006年后，老龄化趋势加快，预计在2015年达到12%，此后我国将进入人口老龄化迅速发展时期。预测显示从2015到2040年间里，中国老年人比例会增加一倍左右。

关键词：干预因子，人口年龄分布，微分方程，老龄化，人口结构，出生率

 一、问题的重述

我国是一个人口大国，人口问题对我国的经济发展有着很大的影响，是制约我国经济发展的关键因素之一。根据最新的“2005年全国%1人口抽样调查公报”指出，我国人口和计划生育工作进入稳定低生育水平、统筹解决人口问题、促进人的全面发展的新阶段。要统筹解决人口与社会发展的各种问题，必须对目前和将来一段时间内的人口问题进行科学的判断，如城市化问题、老龄化问题、出生人口性别比持续升高问题。因此要制度科学的人口政策需要对未来中国的人口变化进行科学预测。

    关于中国人口的预测，目前已经有了许多的模型。根据题意，本文从中国的实际情况和人口增长的特点，建立对中国人口增长的中短期和长期的模型预测。

二、问题的分析  

由于影响人口增长的因素比较复杂，从中国的实际情况和人口增长的特点如果一开始把众多因素都考虑进去，将会对模型的求解很困难。所以我们首先考虑影响人口变化的一些基本因素，如出生率、死亡率、人口年龄的结构密度等。其中人口的年龄结构将在未来一段时间内，产生很大的影响。在人口预测中人口按年龄分布状况是十分重要的，因为不同的年龄人的生育率和死亡率有很大的差别，未来的人口发展状况将大不一样。下图是我们根据附表2的数据计算出的01-05年人口年龄密度图。

我们从图2-1中可看到01-05年年龄结构中，老年人数量有增大的趋势。对此我们可预见我国人口结构有老龄化趋势。在模型的建立中我们可看到这一点。

图2-1

在解决变化率的问题时，常常引用微分方程的方法[1]。在此基础上，我们考虑人口按年龄的分布，即除了时间变量外，年龄是另一自变量，来建立微分方程模型。考虑各种微分方程的方法，指数增长模型和阻滞增长模型(Logistic人口增长模型)显然考虑的因素太少，很难做出合适的模型。宋健模型中考虑的问题比较全面，可以取其优秀的地方学习。可是宋健模型假定死亡率不随时间t的变化而改变，未考虑各种因素对死亡率的影响。

由此，我们考虑到在宋健模型的基础上加以深化，考虑人口结构、城市化、老龄化、社会发展等因素对死亡率、生育率的影响，用变化的死亡率和生育率来预计未来的人口结构等。从而实现中国人口增长的中短期和长期的模型预测。

三、模型的基本假设和符号说明
（1）模型的假设

1  我们研究的人口是一个整体系统

2  影响人口的变化因素只有时间、婴儿的出生、人口的死亡、人口的迁移

3  期间不考虑突发事件（大规模的流行疾病、战争等）对人口变化的影响

4  系统是封闭的，不考虑出国与海归人士对我国人口的影响

（2）模型符号的说明和名词解释

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）     表示时刻t单位时间内平均每个育龄女性的生育数

h（r,t）    表示一个女性在r岁时的生育概率

u（r,t）    表示在时刻t年龄r的人的死亡率

四、模型的建立与求解
 

1 模型的建立

我们记u(r,t)为时刻t年龄r的人的死亡率，其含义是，u（r，t）p（r，t）dr表示时刻t年龄在[r，r+dr]内单位时间死亡的人数。

为了得到p（r,t）满足的方程，考察时刻t年龄在[r,r+dr]内的人到时刻t+dt的情况。他们中活着的那一部分人的年龄变为[r+dr1，r+dr+dr1]，这里dr1=dt。而在dt这段时间内死亡的人数为u（r,t）p(r,t)drdt.

于是

上式可写为：两边同除以drdt，于是：

 =-u(r,t)p(r,t)

 取极限：

=-u（r,t）p（r,t），                       这即是人口密度函数p（r,t）的一阶偏微分方程。

对方程（1）有两个定解条件：

1初始密度函数记为 p（r,0）

2单位时间内出生的婴儿数记为p(0,t)=f(t)

     (2)

出生儿童的模型

在育龄区间[r1，r2]内

 为t时刻平均每个r岁的女性的生育数

则在t时刻，出生儿童的总数为：

（2）（3）方程就是我们建立人口预测模型的基础。

2 模型的求解

当u（r,t）不依赖于t，即在社会安定的局面下和不太长的时间内，死亡率大致与时间无关，于是对方程（2）求得：

 由方程（3）（4）我们利用已知的01-05年的数据对其交替求解

设t=1，2，3，4，5       令2001年为t=1，此后依次类推

则根据附表2中拟和出年龄密度分布

p(r,1),p(r,2),p(r,3),P(r,4),p(r,5)

因此，我们可以依次推出

t=6，

t=7,

t=8 ，  ……

依次计算，就可以预测以后的人口总数和年龄结构。

考虑到其他一些因素的影响，在具体的计算过程中加入干预因子优化模型。例如，考虑到城市化过程对人口分布的影响（城市城镇居民比例增加、乡村居民比例减少,如图），引入分布因子的概念，根据目前的城市化进程和国家制定的城市化政策规划，我们设定了如下变化的因子：

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即是我国未来一段时间内，人口的变化趋势。国家统计局根据2005年全国人口1%抽样调查结果推算2005年底人口为130756万人。根据我们的模型预测结果，我国总人口在2010年将达到13.565亿人。预计到2028年将达到峰值，约为15.17亿，此后开始缓慢下降。如果长期保持1.8的生育水平，则到2040年达到14.843亿。

图4-2 我国未来50年（2005年起）人口变化曲线（预测图）

其预测的中国人口发展趋势为（如表4-1）：

                                                                                                                                      
2、人口年龄结构的变化：

由于生育率的快速下降和长期的低生育水平，我国未来年龄结构变化将体现出两大特征[3]：显著的的‘人口红利’（图4-2）和人口老龄化（图4-3）。‘人口红利’是指人口转变过程中出现的劳动年龄人口比例上升带来的经济增长效应。它的时间长度与程度取决于死亡率、生育率的下降的速度和时间差。他将在未来对我国的就业带来很大的压力和挑战。

    我国未来发展的另一大挑战是持续、快速的人口老龄化。无论是老年人规模还是其比例都在迅速的上升。下面的图我们可以看出。

    人口的老龄化将导致抚养比不断的提高，对社会的保障体系和公共服务体系的压力加大，并影响到社会代际关系的和谐。

    同时由于城镇化的进程加快，进城打工人数的增加，将使得农村的老龄化形式更为严峻。

五、模型的评价

1、模型的优点：

    本模型在建立的过程中，考虑了多个因素（出生率、死亡率、各个年龄男女比等）比以往的人口预测模型如指数增长模型、阻滞增长模型（Logistic模型）有了很大的改进。由于我们是根据人口密度推总人口分布的，人口与时间和年龄都有联系，所以我们可以推出将来的人口的年龄结构。

    此外，我们的模型还涉及到女性性别比、总和生育率，这对我国的计划生育政策提供依据和手段。

   对今后一段时间内的预测很好的吻合了我国最新的国家人口战略发展研究报告。

2、模型存在的问题：

从控制论的观点看，在方程（2）描述的人口系统中p(r,t)可视为状态变量，p(0,t)=f(t)视为控制变量，是分布参数系统的边界控制函数。式表明控制输入中含有状态变量，形成反馈，视为反馈增益。并且通常是正反馈，即人口密度函数p(r,t)增加，通过婴儿出生率f(t)又使p(r,t)进一步增大。方程（4）式中因子f(t-r)表明这种反馈有很大的滞后作用，所以一旦政策失误，使p(r,t)在一段时间内增长得过多过快，再想同过控制手段和h(r,t)把人口增长降下来，就很困难，可能需要几代人。显然，随着科学的发展，统计预测的加强，政策方面的大失误，我们暂不考虑。

 参考文献

[1]姜启源  谢金星等，数学模型，高等教育出版，2003

[2]苏金明  阮沈勇    matlab实用教程，2005

[3]陈伟  中国未来人口发展趋势：2005～2050   人口研究  2006.7（4）

[4]刘祯 张秋辉 我国人口发展现状研究 科技统计 河南科技 2007.5（上）

源程序附录

%pp.m

load a01-05.mat

for r=1:91

    p1(r,1)=(a2001(r,2)+a2001(r,4));

    p2(r,1)=(a2001(r,6)+a2001(r,8));

    p3(r,1)=(a2001(r,10)+a2001(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)=(a2002(r,2)+a2002(r,4));

    p2(r,2)=(a2002(r,6)+a2002(r,8));

    p3(r,2)=(a2002(r,10)+a2002(r,12));

    p(r,2)=(p1(r,2)*(z01_05(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)=(a2003(r,2)+a2003(r,4));

    p2(r,3)=(a2003(r,6)+a2003(r,8));

    p3(r,3)=(a2003(r,10)+a2003(r,12));

    p(r,3)=(p1(r,3)*(z01_05(2,3)+z01_05(3,3))+p2(r,3)*(z01_05(4,3)+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)=(a2004(r,2)+a2004(r,4));

    p2(r,4)=(a2004(r,6)+a2004(r,8));

    p3(r,4)=(a2004(r,10)+a2004(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_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)=(a2005(r,2)+a2005(r,4));

    p2(r,5)=(a2005(r,6)+a2005(r,8));

    p3(r,5)=(a2005(r,10)+a2005(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_05(4,5)+z01_05(5,5)+z01_05(6,5)+z01_05(7,5));

end;

krt=0.460012;

bata=1.8;

swl=sum(a2005(:,2).*a2005(:,3))+sum(a2005(:,4).*a2005(:,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:90

        if t<=r

             jf1=((sum(a2005((r-t+1):(r+1),2).*a2005((r-t+1):(r+1),3))+sum(a2005((r-t+1):(r+1),4).*a2005((r-t+1):(r+1),5)))*io1(t+1)+(sum(a2005((r-t+1):(r+1),6).*a2005((r-t+1):(r+1),7))+sum(a2005((r-t+1):(r+1),8).*a2005((r-t+1):(r+1),9)))*io2(t+1)+(sum(a2005((r-t+1):(r+1),10).*a2005((r-t+1):(r+1),11))+sum(a2005((r-t+1):(r+1),12).*a2005((r-t+1):(r+1),13)))*io3(t+1))/100000-lp(t+1)+sw(t+1);

            if t>4

                 p(r+1,t+1)=p(r-t+1,1)*exp(-jf1);

            end;

        end;

        if t>r

            jf2=((sum(a2005(1:(r+1),2).*a2005(1:(r+1),3))+sum(a2005(1:(r+1),4).*a2005(1:(r+1),5)))*io1(t+1)+(sum(a2005(1:(r+1),6).*a2005(1:(r+1),7))+sum(a2005(1:(r+1),8).*a2005(1:(r+1),9)))*io2(t+1)+(sum(a2005(1:(r+1),10).*a2005(1:(r+1),11))+sum(a2005(1:(r+1),12).*a2005(1:(r+1),13)))*io3(t+1))/100000-lp(t+1)+sw(t+1);

            if t>4

                p(r+1,t+1)=f(t-r+1)*exp(-jf2);

            end;

        end;

    end;

    for j=15:49

        f(t+1)=f(t+1)+(a2005(j+1,15))*krt*p(j+1,t+1)/100000;

    end;

    f(t+1)=bata*f(t+1);

    if f(t+1)<0.002

        f(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:50

    rk(i+1)=rk(i)*(1+zrzzl(i+1));

end;

t=0:50;

plot(t,rk,"*");

grid on;

hold on;