用偏微分方程进行人口仿真.docx
《用偏微分方程进行人口仿真.docx》由会员分享，可在线阅读，更多相关《用偏微分方程进行人口仿真.docx（24页珍藏版）》请在冰豆网上搜索。
用偏微分方程进行人口仿真
系统仿真课程设计
题目：
专业：
小组成员：
用偏微分方程进行人口仿真
摘要：
建立中国人口增长的数学模型，由建立的人口发展的偏微分方程来预测中国未来人口的数量和结构。
在预测的基础上，考虑到降低生育率与人口数量和老龄化有着直接的关系，所以在预测人口基础之上，我们进一步拓展对未来人口控制进行研究。
即在对人口数量预测的同时对其控制及其优化做出探讨。
关键词：
一、提出背景
人口增长预测是随着社会经济发展而提出来的。
由于人类社会生产力水平低，生产发展缓慢，人口变动和增长也不明显，生产自给自足或进行简单的以货易货，因而对未来人口发展变化的研究并不重要，根本不用进行人口增长预测。
人口是社会经济活动的主体,人口的发展变动趋势,对社会经济发展的影响关系极大,因此人口预测在社会经济实践中占有十分重要的地位。
现阶段，中国在享受计划生育政策带来的红利的同时，依然面连着人口结构性失调的严重性问题，而当今社会，经济发展迅速，生产力达到空前水平，这时的生产不仅为了满足个人需求，还要面向社会的需求，所以必须了解供求关系的未来趋势。
而人口增长预测是对未来进行预测的各环节中的一个重要方面。
准确地预测未来人口的发展趋势，制定合理的人口规划和人口布局方案具有重大的理论意义和实用意义。
二、问题重述与分析
1）、基本假设
1．假设本问题所使用的数据均真实有效，具有统计分析价值。
2．假设本问题所研究的是一个封闭系统，也就是说不考虑我国与其它国家的人口迁移问题。
3.无重大毁灭性自然灾害和疾病，无战争等暴烈活动，即扰动
人口发展的因素只有人。
4.在对人口进行分段处理时，假设同一年龄段的人死亡率相同，同一年龄段的育龄妇女生育率相同。
5.生育模式在预测时间内保持不变，并且假设一胎只生一个孩
子。
6．人类的生育观念不发生太大改变，如没有集体不愿生小孩的想法
7.中国各地各民族的人口政策相同。
8.人口生存环境为一般常态的自然和社会环境。
9.中短期内，总和生育率、死亡率和出生性别比不发生大的波
动。
2）、符号说明
r:
表示人的年龄；
rm：
表示人的最高年龄；
t:
表示时间；
F（r,t）：
时刻t年龄小于r的人口的分布函数；
p（r,t）dr：
表示时刻t年龄在区间[r,r+dr）内的人数；
p（r,t）=：
表示人口密度函数；初始密度函数（r）；
μ（r,t）：
表示时刻t年龄r的人的死亡率；
β（t）：
表示总和生育率或生育胎次；
f（t）：
表示婴儿出生率；
h（r,t）：
表示年龄为r女性的生育加权因子，称生育模式；
N（t）：
表示人口总数函数；R（t）：
表示平均年龄函数；
S（t）：
表示平均寿命函数；w（t）：
表示老龄化指数。
3）、数据来源
数据来自。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
三、整体思路及主要涉及因素
本次课设是一个关于人口预测的问题，与其他不同，本问题需要根据中国特殊的国情在特有的现状之下去研究，我们根据对问题的分析并结合实际情况认为对人口产生主要影响的因素有以下四个：
生育率、死亡率、年龄结构、男女比例。
在这里需要说明的是对于人口产生影响的一些因素，如经济发展状况，生态环境情况、已婚夫妇对生育所持的态度、医疗技术的发展等，我们认为它们对人口的增长是通过作用于以上四个指标而间接发挥作用的。
而对于诸如战争爆发、疾病流行等突发因素，由于其不可预测性，我们不考虑。
在考虑上述主要因素的前提下，我们根据不同的因素所对应的不同情况建立与其最接近的偏微分方程，从而在综合情况下对未来人口的总趋势做出一个预测。
以下对四种主要因素做一个简单概述：
1．生育率
生育率代表育龄妇女生育人口的能力，从一定意义上讲生育率的高低控制着人口增长率高低，通常来说生育率越高人口增长率越高，所以说生育率是人口增长的源头。
生育率的影响因素很多，首先是年龄因素，不同年龄段的育龄妇女的生育率不同，通常20岁至30岁的育龄妇女的生育率最强；
此外是地域因素，受政策因素、观念认识、周边环境等影响乡村育龄妇女的生育率高于城市育龄妇女的生育率；还有其它因素的影响，比如大规模疾病会降低育龄妇女的生育率。
2．死亡率
死亡率表示一定时期内一个人口群体中死亡的人数占该人口群体的比值，和生育率一样死亡率的高低同样控制着人口增长率高低，
如果说生育率是人口增长的源头，则死亡率是人口增长的汇点。
同样影响死亡率的因素很多，首先不同年龄段的死亡率不同，通常老年人和刚出生的婴儿的死亡率较高；从长远来看，随着医疗水平的提高，整个人口群体的死亡率将会成下降趋势；此外一些突发事件，如战争、疾病等，将会使使那一段的人口死亡率大幅度提高。
3．年龄结构
年龄结构反映了总体人口在各年龄段分布情况，年龄结构蕴涵的信息量很大，从其中我们可以实现对很多问题的分析，比如从年龄结构我们可以分析出社会的老年化程度，此外从年龄结构我们可以判断出不同时间段人口出生的情况，比如年龄结构不仅反映了总体人口在各年龄段分布情况，而且考虑到不同年龄段人口生育率、死亡率不同等情况，我们可以在年龄结构中有效反映这些差异。
4．男女比例
男女比例反映了总体人口中男性与女性人数的比较关系，男女比例值能反映出体人口中男性与女性人数是否协调，男女比例主要受男女出生比和男女死亡率的影响，男女出生比正常范围在103－107，也就是说出生100个女儿的同时会有103—107个男儿出生，但是在现实社会中，女性死亡率低于男性，所以男性与女性人数大致相等，社会维持在一个稳定状态。
但目前我国男女出生比超过110，这不仅将导致男女比例失调，还会对人口的预测产生影响，所以在人口预测时必须将男女比例问题考虑进去。
四、模型建立
1）、人口发展的偏微分方程建立
造成人口数量和结构变化的因素不外乎出生、死亡和迁移。
为简化，预测期间只考虑自然的出生与死亡，不计迁移等社会因素的影响。
为研究任意时刻不同年龄的人口数量，引入人口的分布函数和密度函数。
时刻t年龄小于r的人口称为人口分布函数，记F（r,t），其中t，（r>0），均为连续变量，设F连续、可微，时刻t的人口总数记作N（t），最高年龄记作rm，理论推导时设rm→∞，于是对非负非降函数F（r,t）有
（1）
人口密度函数定义：
（2）
其中，p（r,t）dr表示时刻t年龄在区间[r,r+dr）内的人数。
记μ（r,t）为时刻t年龄r的人的死亡率，其含义是，μ（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这段时间内死亡的人数为μ（r,t）p（r,t）drdt。
于是
（3）
变为：
dr1=dt可得：
（4）
即人口密度函数p（r,t）的一阶偏微方程，死亡率μ（r,t）为已知函数。
方程（4）有两个定解条件：
初始密度函数记为p（r,0）=p0（r），是已知函数;单位时间出生的婴儿数记为p（0,t）=f（t）,是婴儿出生率。
故所得模型为：
（4）
（5）
（6）
方程（5）求解过程比较复杂，在这儿只给出一种特殊情况求解。
社会安定且不太长的时间内，死亡率大致与时间无关，故可近似的设μ（r,t）=μ（r），可得（5）的解为：
（7）
这个解在Otr平面上有一个前线的解释：
图中对应的对角线r=t将Otr面的第一象限分为两部分，tr区域内，p（r,t）则由未来的生育状况f（t-r）和死亡率μ（s）（0≤s2）相关主要因素的方程建立和定义
（5）或（7）中p0（r）和μ（r）可从人口统计数据中得到，μ（r,t）也可从μ（r,0）粗略得到。
这样，为了预测和控制人口的发展状况，人们主要关注和可以用作控制手段的就是婴儿出生率f（t）,以下对f（t）进一步分解：
记女性性别比函数为k（r,t）,即时刻t年龄在[r,r+dr）的女性人数为k（r,t）p（r,t）dr,将这些女性单位时间内平均美人的生育数记为b（r,t）,设育龄区间为[r1,r2],则：
（8）
定义：
（9）
于是
β（t）=（10）
f（t）=β（t）（11）
h（r）=,r>r1（12）
并取θ=2，α=n/2，此时有：
rc=r1+n-2（13）
可以看出，提高r1意味着晚婚，而增加n意味着晚育。
这样，人口发展方程（5）和单位时间出生的婴儿数f（t）的表达式（12），构成了我们的连续型人口模型。
模型中死亡率函数μ（r,t）、性别比函数k（r,t）和初始密度函数p0（r）可以有人口统计资料直接得到，或在资料的基础上估计，而生育率β（t）和生育模式h（r,t）则是可以用于控制人口发展过程的两种手段，β（t）可以控制生育的多少，h（r,t）可以控制生育的早晚或疏密。
在上面的模型分布函数F（r,t）固然是人口发展过程最完整的描述，但使用起来并不方便。
下面定义一些人口指数的及它们与p（r,t）
等数量之间的关系。
1）人口总数N（t）
（15）
2）平均年龄R（t）
（15）
3）平均寿命S（t）它表示时刻t出生的人不论活到什么时候，
死亡率都是按时刻t的μ（r,t）计算，这些人的平均存活时间
dτ（16）
S（t）实际上是预估寿命。
通常说目前平均寿命已达多少岁，是指今年
出生婴儿的预估寿命，即S（0）。
根据统计资料得到当前的死亡率
μ（r,0）后就可以算出S（0）.
4）老龄化指数w（t）定义：
（17）
显然，平均年龄R（t）越大w（t）越大；对于R（t）相同的两个国家
或地区，平均寿命S（t）大的，表示健康水平高，一个人能工作的时间一生中占得比例越大，于是老龄化的指数w（t）较小。
五、数据处理与仿真
1）参数的确立
1.初始密度函数（r）的确定
初始密度选择的是2010年的人口数据。
之所以没有选择最近的数据，是因为这样更方便对预测数据的检验。
*人口单位数为千人
2.死亡率函数μ（r,t）的确定
在社会和平，不出现大的瘟疫等突发情况下，可认为人口的死亡率基本不变。
我们分别对2008--2012年和世界发达国家的各年龄段分性别的死亡率做加权平均，作为我国未来人口的死亡率。
3.性别比例函数k（r,t）的确定
在未来一段时期，我们认为各年龄段性别比例基本保持不变。
我们根据2008--2012年的人口统计数据，采取平均计算得出各年龄段的性别比例，将之作为计算所用的性别比例。
4.生育模式函数h（r,t）的确定
生育模式是对年龄为r的女性的生育加权因子，在社会稳定情况下可认为其基本不变。
根据惯例，我们认为生育年龄为15--49岁。
根据2008--2012年的统计数据，分别对各年龄段的生育率取平均，然后再做归一化处理，将得到的数据作为计算所需的生育模式函数。
2）方程求解
先取β（t）=1.5，根据方程（5）和已确定的参数，求出每隔5年的各个年龄段的数据，然后再根据方程（15），（16），
（17），（18）分别求出人口总数和老龄化指数。
然后改变β（t）的值，重复以上步骤计算。
（运行renkou1.m程序）分别调整不同时间阶段的β（t）的值，然后再计算求出各个年龄段的数据，人口总数和老龄化指数。
将取得的数据作横向比较，选出较优结果。
六、计算结果与分析
β取不同值时的各个时间段的人口总数和老龄化指数（表1）
+
β取不同值时的人口总数和老龄化指数最大值（表2）
β
总人口数最大值（千人）
老龄化最大值
1.0
1.419e+006
0.7088
1.2
1.458e+006
0.6759
1.4
1.498e+006
0.648
1.5
1.522e+006
0.6346
1.6
1.545e+006
0.6216
1.8
1.607e+006
0.6025
2
1.688e+006
0.5856
以下是β取1.5和1.8所对应的人口增长曲线和老龄化指数曲线（由于图形太多，只取其中两个具有代表性的）：
图1
图2
根据相关文献，中国目前的β值约为1.5，由表一也可知β＝1.5时的计算结果与真实的误差是最小的。
此时人口的峰值约为15.52亿，老龄化指数的峰值为0.6346β＝1.8时人口的峰值约为16.07亿。
也就是说，当每个妇女平均生育的小孩数从1.5提高到1.8时，整个国家就会多出生九千万人。
所以我们现在不应该也不能放开二胎的限制。
由表２可知，没有一组数据同时满足人口峰值小于15亿和老龄化峰值小于0.6，下面我们通过调整生育政策来调整人口结构。
七、实验拓展
由上述可知，最大老龄化指数是与最大人口数呈反相关的。
我们可以分别调整不同时间阶段的β（t）的值，以达到优化人口总数和老龄化指数的目的。
下面我们分别计算（运行youhua1.m）
2010年β（t）＝1.5，2020年β（t）＝1.8
2010年β（t）＝1.5，2030年β（t）＝1.8
2010年β（t）＝1.5，2030年β（t）＝2.0
图3
2010年β（t）=1.52020年β（t）=1.8
图4
2010年β（t）=1.52030年β（t）=1.5
图5
2010年β（t）=1.52030年β（t）=2.0
经分析可知，图3和图4均满足要求。
它们的人口峰值和老龄化指数峰值分别为１４．9５亿、0.597和14.79亿、0.596.从峰值上比较，二者并无太大区别。
但2100年时，后一幅图的值比前者多出近一个亿，达到12.5亿。
所以第一幅图代表方案更优，即2010--2020年每个妇女平均生育小孩数控制在1.5附近，2020--2100年生育政策可以适度放松，将每个妇女平均生育小孩数控制在1.7附近，这样人口将逐步优化。
八、结论概述
由以上仿真结果可以看出，在忽略一些不可控和不可预测的因素条件下，我们可以根据当时所特有的情况确定符合的β值，从而预测某年的人口总数和老龄化，对未来人口数量和结构有个了解。
而在改变β值得取值范围我们也可以对未来人口的生育率有所了解，在此基础之上，可以采取相应的措施对我国人口的发展趋势有所掌控。
以上结果表明，在中短期，仍然要坚持严格的人口政策，中长期可以适度放松，但是每个妇女平均生育小孩数仍要控制在2.0以内，否则人口的峰值会很大，超出环境所能承受的极限。
附录
程序及其说明
1renkou1.m
clc,clf,clear
b=1.5%每个妇女一生生的小孩数
r=[75533708827490899889127413101014971381180161247541055957875381312586674111332972238521337356321984]';%初始人口数据
x=[119.1250118.6560116.2410108.1700100.9540101.3160104.0020104.7820104.0190103.7790105.1400102.1210103.4730101.884099.001089.705079.369060.095049.5520];%性别比例
sh=[0005.9369.4784.0845.8418.717.514.68000
000000];%出生率
sm=[1.340.350.380.540.680.841.131.632.343.345.627.9112.6320.9736.5158.2193.61139.51200.60];%男人死亡率
sw=[1.280.230.230.260.300.360.510.731.101.602.894.237.2512.8124.1440.2470.81110.72180.20];%女人死亡率
n=length（r）;
shsum=sum（sh）;
h=b*sh/shsum;
fori=1:
n
nv（i）=100/（100+x（i））;%归一化性别比例
h（i）=h（i）*nv（i）;
end
fori=1:
n%计算死亡率
sm（i）=sm（i）*x（i）;
x（i）=x（i）+100;
sw（i）=sw（i）*100;
s（i）=（sm（i）+sw（i））/x（i）;
s（i）=s（i）/1000;
s（i）=s（i）/4;
end
A=zeros（19）;
A（1,:
）=（1-s
（1））*h;
fori=2:
n%计算以后每年的人口数
s（i）=（1-s（i））^5;
A（i,i-1）=s（i）;
end
M=zeros（19,20）;%储存未来100年的人口数
M（:
1）=r;
fori=2:
20
M（:
i）=A*M（:
i-1）;
end
zsh=sum（M）;%人口总数
[zshmax,k]=max（zsh）%求最大人口总数及其下标
NL=zeros（19）;
fori=1:
19
NL（i,i）=5*i-3;
end
N=NL*M;
zn=sum（N）;
fori=1:
20
zn（i）=zn（i）/zsh（i）;%平均年龄
shm（i）=70+i;%预期寿命
w（i）=zn（i）/shm（i）;%老龄化指数
end
[wmax,p]=max（w）%求w的最大值及其下标
fori=1:
20
nf（i）=2005+5*i;
end
subplot（1,2,1）
[nf（k）zshmax];
plot（nf,zsh,nf（k）,zshmax,'ro'）,set（gca,'fontsize',15）
text（nf（k）,zshmax,['\bullet\leftarrowTheMaxValue
is',num2str（zshmax）,'at\it{nf}=',num2str（nf（k））]）;
gridon;
xlabel（'年份'）;
ylabel（'人口总数'）;
title（'\it{人口增长曲线b=1.5}'）
subplot（1,2,2）
[nf（p）wmax];
plot（nf,w,nf（p）,wmax,'ro'）,set（gca,'fontsize',15）
text（nf（p）,wmax,['\bullet\leftarrowTheMaxValue
is',num2str（wmax）,'at\it{nf}=',num2str（nf（p））]）;
gridon;
xlabel（'年份'）;
ylabel（'老龄化指数'）;
title（'\it{老龄化指数变化曲线b=1.5}'）
zsh,w
2.youhua1.m
clc,clf,clear;
b=1.5,b1=1.8,%每个妇女一生生的小孩数
j=7;%从（2010+5*（j-1））年改变生育政策
r=[75533708827490899889127413101014971381180161247541055957875381312586674111332972238521337356321984]';%
初始人口数据
x=[119.1250118.6560116.2410108.1700100.9540101.3160104.0020104.7820104.0190103.7790105.1400102.1210103.4730101.884099.001089.705079.369060.095049.5520];%性别比例
sh=[0005.9369.4784.0845.8418.717.514.68000
000000];%出生率
sm=[1.340.350.380.540.680.841.131.632.343.345.627.9112.6320.9736.5158.2193.61139.51200.60];%男人死亡率
sw=[1.280.230.230.260.300.360.510.731.101.602.894.237.2512.8124.1440.2470.81110.72180.20];%女人死亡率
n=length（r）;
shsum=sum（sh）;
h=b*sh/shsum;
fori=1:
n
nv（i）=100/（100+x（i））;%归一化性别比例
h（i）=h（i）*nv（i）;
end
fori=1:
n%计算死亡率
sm（i）=sm（i）*x（i）;
x（i）=x（i）+100;
sw（i）=sw（i）*100;
s（i）=（sm（i）+sw（i））/x（i）;
s（i）=s（i）/1000;
s（i）=s（i）/4;
end
A=zeros（19）;
A（1,:
）=（1-s
（1））*h;
fori=2:
n%计算以后每年的人口数
s（i）=（1-s（i））^5;
A（i,i-1）=s（i）;
end
M=zeros（19,20）;%储存未来100年的人口数
M（:
1）=r;
fori=2:
20
M（:
i）=A*M（:
i-1）;
end
zsh=sum（M）;%人口总数
[zshmax,k]=max（zsh）;%求最大人口总数及其下标
NL=zeros（19）;
fori=1:
19
NL（i,i）=5*i-3;
end
N=NL*M;
zn=sum（N）;
fori=1:
2
zn（i）=zn（i）/zsh（i）;%平均年龄
shm（i）=70+i;%预期寿命
w（i）=zn（i）/shm（i）;%老龄化指数
end
[wmax,p]=max（w）;%求w的最大值及其下标
x=[119.1250118.6560116.2410108.1700100.9540101.3160104.0020104.7820104.0190103.7790105.1400102.1210103.4730101.884099.001089.705079.369060.095049.5520];%性别比例
sh=[0005.9369.4784.0845.8418.717.514.68000
000000];%出生率
sm=[1.340.350.380.540.680.841.131.632.343.345.627.9112.6320.9736.5158.2193.61139.51200.60];%男人死亡率
sw=[1.280.230.230.260.300.360.510
升级会员 