数学建模:人口增长模型

模型目标:
通过给定的一组人口增长数据,预测后续的人口增长情况.

一、指数增长模型

假设增长率不变:

若已知人口年增长率为r,今年人口为 x0x_0x0,预测k年后的人口可以用简单的公式得到:
xk=x0(1+r)kx_k = x_0(1+r)^kxk=x0(1+r)k
*以美国人口为例,数据点取下表:

用matlab输入好数据:

p = [3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4];

该公式得到的人口增长模型如下图(matlab作图):

选取的增长率0.2与实际的人口数据点贴合度差,于是进一步建立模型,找到偏差更小的增长率r

精确化增长率r

把人口看作关于时间的可微函数 x(t)x(t)x(t),记 x(0)=x0x(0) = x_0x(0)=x0
rx0rx_0rx0 即为单位时间 x(t)x(t)x(t) 的增长量
所以得到微分方程:
dxdt=rx,x(0)=x0\frac{dx}{dt} = rx , x(0)=x_0dtdx=rx,x(0)=x0
解得:
x(t)=x0ertx(t) = x_0e^{rt}x(t)=x0ert
上述公式即位指数增长模型

数据拟合:

接下来对指数增长模型的参数进行估计,即数据拟合

·法一

用人口数据和线性最小二乘法,对上式取对数:
y=rt+a,y=lnx,a=lnx0y = rt +a , y=lnx,a=lnx_0y=rt+a,y=lnx,a=lnx0
用matlab进行编程拟合:

%取1790年为t=0,数据点取到2000年
t=0:10:210;
lnp = log(p);
cftools

解得:
r=0.202,x0=e1.8=6.049r = 0.202,x_0 = e^{1.8} = 6.049r=0.202,x0=e1.8=6.049
带入得x(t)=6.049e0.0202tx(t) = 6.049 e^{0.0202t}x(t)=6.049e0.0202t

·法2
对人口数据做数值微分,计算平均值r‘,x0直接选用原始数据.
函数在各点的近似导数值为(数值微分中点公式):
x′(tk)=xk+1−xk−12Δt,(k=1,2,3,...,n−1)x'(t_k) = \frac{x_{k+1}-x_{k-1}}{2\Delta t},\\\ \\(k=1,2,3,...,n-1)x′(tk)=2Δtxk+1−xk−1, (k=1,2,3,...,n−1)
x′(0)=4x1−3x0−x22Δt,x′(n)=−4xn−1+3xn+xn−22Δtx'(0) = \frac{4x_1-3x_0-x_2}{2 \Delta t}, x'(n) = \frac{-4x_{n-1}+3x_{n}+x_{n-2}}{2 \Delta t}x′(0)=2Δt4x1−3x0−x2,x′(n)=2Δt−4xn−1+3xn+xn−2
那么增长率为: r(tk)=x′(tk)x(tk)r(t_k)=\frac{x'(t_k)}{x(t_k)}r(tk)=x(tk)x′(tk)

*公式相关推导可参考数值计算方法第六章数值积分和数值微分

增长率 rk=x′(tk)x(tk)r_k=\frac{x'(t_k)}{x(t_k)}rk=x(tk)x′(tk)再取平均值得到r = 0.0205

改进的指数增长模型

上述模型对于增长率r不变的假设导致预测曲线与实际偏差较大
所以改进模型中假设 r 为 t 的函数 r(t)r(t)r(t) , 根据上述法2的 x′(t)x'(t)x′(t) 画r—t图像:

r=[];
for i=1:22if i == 1r(i)=(4*p(i+1)-3*p(i)-p(i+2))/(20*p(i));elseif i == 22r(i)=(-4*p(i-1)+3*p(i)+p(i-2))/(20*p(i));elser(i)=(p(i+1)-p(i-1))/(20*p(i));end
end
plot(t,r,'.','MarkerSize',20);
ylim ([0.005,0.04]);
xlabel('t');
ylabel('增长率r');

根据散点图假设r(t)=r0−r1tr(t)=r_0-r_1tr(t)=r0−r1t的线性函数,用最小二乘法线性拟合得到:
r0=0.03252,r1=0.0001143r_0 = 0.03252,\space r_1 = 0.0001143r0=0.03252, r1=0.0001143
根据微分方程:
dx/dt=r(t)xdx/dt=r(t)xdx/dt=r(t)x
⇒x(t)=x0e(r0t−r1t2/2)\Rightarrow x(t) = x_0 e^{(r_0t-r_1t^2/2)}⇒x(t)=x0e(r0t−r1t2/2)
把拟合后的参数带入:

xt2 = 3.9.*exp(0.03252.*t-0.0001143.*t.^2./2);
plot(t,xt2,'LineWidth',1);

最终根据改进模型预测的2010年人口为290million,与实际数据281.4吻合度较高.显然改进后的模型优于前两者.

二、logistic模型

改进的指数增长模型中增长率线性下降,但没有体现其下降的相关影响因素,只是以时间为变量.logistic模型考虑了自然资源,环境等对人口增长的阻滞作用.

未完…