数学建模

 

模型知识点

变化模型-微分方程-人口增长模型

 

Wang Haihua

🚅 🚋😜 🚑 🚔

人口增长模型

认识人口数量的变化规律，建立人口模型，作出较准确的预报，是有效控制人口增长的前提。

下表给出了近两个世纪的美国人口统计数据，我们以此研究人口增长模型。

import pandas as pd
df = pd.DataFrame()
df["年份"] = [i for i in  range(1790,2010,10)]
df["人口"] = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,
50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204,226.5,251.4,281.4]
df.T

	0	1	2	3	4	5	6	7	8	9	...	12	13	14	15	16	17	18	19	20	21
年份	1790.0	1800.0	1810.0	1820.0	1830.0	1840.0	1850.0	1860.0	1870.0	1880.0	...	1910.0	1920.0	1930.0	1940.0	1950.0	1960.0	1970.0	1980.0	1990.0	2000.0
人口	3.9	5.3	7.2	9.6	12.9	17.1	23.2	31.4	38.6	50.2	...	92.0	106.5	123.2	131.7	150.7	179.3	204.0	226.5	251.4	281.4

2 rows × 22 columns

先简单通过散点图看一下人口的增长模式。

import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')  #设置绘图风格为ggplot
matplotlib.rcParams['font.size'] = 15  #设置全文绘图字体大小为15号字

year = [i for i in range(1790,2010,10)]
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,
50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204,226.5,251.4,281.4]
plt.scatter(year,population)   #绘制散点图
plt.xlabel("Year")
plt.ylabel("Population(Million)")

Text(0, 0.5, 'Population(Million)')

指数增长模型

最简单的人口增长模型是公认的：记今年人口为$x_0$,$k$年后人口为$x_k$， 年增长率为$r$，则

$$ x_{k}=x_{0}(1+r)^{k} $$

显然，这个公式的基本条件是年增长率$r$保持不变,这个假设显然过于简单了。

200多年前英国人口学家T. Malthus (1766-1834)调查了英国100多年的人口统计资料，得出了人口增长率不变的假设，并据此建立了著名的人口指数增长模型。 记时刻$t$的人口为$x(t)$，当考察一个国家或一个较大地区的人口时，$x(t)$是一个很大的整数。为了利用微积分这一数学工具，将$x(t)$视为连续、可微函数， 记初始时刻$(t=0)$的人口为$x_0$，假设人口增长率为常数$r$，即单位时间内$x(t)$的增量为等于$r$乘以$x(t)$，于是得到$x(t)$满足微分方程

$$ \frac{\mathrm{d} x}{\mathrm{d} t}=r x, x(0)=x_{0} $$

由这个方程,通过分离变量法很容易解出其解析解:

分离变量法求解微分方程

先分离变量

$$ \frac{\mathrm{d} x}{x}=r \mathrm{d} t $$

两边同时积分

$$ \int \frac{\mathrm{d} x}{x} = \int r \mathrm{d} t $$

$$ \ln(x) = rt+C_1 $$

$$ x = e^{rt+C_1} $$

$$ x = Ce^{rt} $$

当$x=0$时，

$$ x(0) = C = x_0 $$

因此$C = x_0$

$$ x(t)=x_{0} e^{r t} $$

$r>0$时，上式表示人口将按指数规律随时间无限增长，称为指数增长模型．

# 指数增长模型的绘图
import numpy as np
x0 = 1  # 初始人口数
r = 0.02    # 增长率
t = [i for i in range(1000)] # 时间列表
x_t = [x0 * np.exp(r * time) for time in t] # 人口增长记录
plt.plot(t, x_t) # 绘图
plt.xlabel("Time") # 加横坐标label
plt.ylabel("Population") # 加纵坐标label

Text(0, 0.5, 'Population')

模型的参数估计、检验和预报

我们来尝试用这个模型，带入19世纪的数据，预测20世纪的人口。

为了估计指数增长模型中的参数$r$和$x_0$，需将原式取对数，得

$$ y=r t+a $$

其中，

$$ y=\ln x, a=\ln x_{0} $$

from sklearn.linear_model import LinearRegression # 导入LinearRegression方法
ln_population = np.log(population)  # 对人口取对数，得到y
plt.scatter(year[0:12],ln_population[0:12])  # 绘制1790 - 1900的y-t图，并进行线性回归

## 使用LinearRegression 进行线性回归
lrModel = LinearRegression()
lrModel.fit(np.array(year[0:12]).reshape(-1,1),ln_population[0:12])
ln_population_fit = [ lrModel.intercept_  +  lrModel.coef_  *i for i in range(1780,1920,10)  ]
plt.plot(range(1780,1920,10),ln_population_fit) # 将线性回归以后的直线绘制在散点图上
plt.xlabel("Year")
plt.ylabel(r'$\ln(x)$')

Text(0, 0.5, '$\\ln(x)$')

## 我们把以上的代码合起来，放在一起。
## 对比一个世纪后的预测值和真实值

#输入原始数据
year = [i for i in range(1790,2010,10)]
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,
50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204,226.5,251.4,281.4]  
plt.scatter(year,population,label = 'Real Data')

#进行参数的线性拟合
from sklearn.linear_model import LinearRegression
ln_population = np.log(population)
lrModel = LinearRegression()
lrModel.fit(np.array(year[0:12]).reshape(-1,1),ln_population[0:12])
ln_population_fit = [ lrModel.intercept_  +  lrModel.coef_  *i for i in range(1780,1920,10)  ]

#将预测值和真实值进行对比
r = lrModel.coef_
x0 = np.exp(lrModel.intercept_)
pop_predicted = [x0 * np.exp(lrModel.coef_ * time) for time in year]
plt.plot(year,pop_predicted,'b',label = 'Predicted')
plt.legend()

<matplotlib.legend.Legend at 0x18b45762d90>

显然，使用1790 - 1900的数据得到的指数模型预测效果不好，不能够准确地预测20世纪的人口增长情况。

• 历史上，指数增长模型与19世纪以前欧洲一些地区人口统计数据可以很好地吻合，迁往加拿大的欧洲移民后代入口也大致符合这个模型。另外，用它作短期人口预测可以得到较好的结果。显然这是因为在这些情况下，人口增长率是常数这个基本假设大致成立。

• 但是长期来看，任何地区的人口都不可能无限增长，即指数模型不能描述、也不能预测较长时期的人口演变过程。这是因为，人口增长率事实上是在不断地变化着排除灾难、战争等特殊时期，一般说来，当人口较少时 ，增长较快，即增长率较大；人口增加到一定数量以后，增长就会慢下来，即增长率变小。

如果根据上面给出的数据计算一下美国人口的年增长率

## 增长率可视化
year = [i for i in range(1790,2010,10)]
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,
50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204,226.5,251.4,281.4]
## 计算增长率
rate = []
for i in range(len(population)-1):
    rate.append((population[i+1] - population[i])/ population[i])
## 可视化
plt.scatter(year[1:],rate)
plt.xlabel("Year")
plt.ylabel(r'$r$')

Text(0, 0.5, '$r$')

可以看到增长率从19世纪开始就基本上在缓慢下降。如果用一个平均的年增长率作为$r$，用指数增长模型描述美国人口的变化，会发现结果与实际数据相差很大。

看来，为了使人口预报特别是长期预报更好地符合实际情况，必须修改指数增长模型关于人口增长率是常数这个基本假设。

阻滞增长模型一logistic模型

分析人口增长到一定数量后增长率下降的主要原因，人们注意到，自然资源、环境条件等因素对人口的增长起着阻滞作用，并且随着人口的增加，阻滞作用越来越大。所谓阻滞增长模型就是考虑到这个因素，对指数增长模型的基本假设进行修改后得到的。

阻滞作用体现在对人口增长率$r$的影响上，使得$r$随着人口数$x$的增加而下降。若将$r$表示为$x$的函数$r(x)$，则它应是减函数。于是

$$ \frac{d x}{d t}=r(x) x, x(0)=x_{0} $$

对$r(x)$的一个最简单的假定是，设$r(x)$为$x$的线性函数，即

$$ r(x)=r-s x(r, s>0) $$

这里$r$称固有增长率，表示人口很少时（理论上是$x = 0$)的增长率．为了确定系数$s$的意义，引人自然资源和环境条件所能容纳的最大人口数量$x_m$，称人口容量．当 $x=x_m$时人口不再增长，即增长率$r(x_m)=0$

于是

$$ r(x)=r\left(1-x / x_{m}\right) $$

代回原方程，得到

$$ \frac{d x}{d t}=r x\left(1-\frac{x}{x_{m}}\right), \quad x(0)=x_{0} $$

(logistic_origin)

方程右端的因子$rx$体现人口自身的增长趋势，因子$(1-\dfrac{x}{x_m})$则体现了环境和资源对人口增长的阻滞作用．

显然．$x$越大，前一因子越大，后一因子越小，人口增长是两个因子共同作用的结果，上式称为阻滞增长模型,同样地，我们使用分离变量法求解这个方程

分离变量法求解logistic模型

原方程为{eq}logistic_origin

分离变量后 $$ \frac{\mathrm{d} x}{x\left(1-\dfrac{x}{x_{m}}\right)}=r \mathrm{d} t $$

等式左边变形 $$ \left( \frac{1}{x} + \frac{1}{x_m - x} \right) \mathrm{d} x =r \mathrm{d} t $$

两边同时积分 $$ \int \left( \frac{1}{x} + \frac{1}{x_m - x} \right) \mathrm{d} x =\int r \mathrm{d} t $$

$$ \ln(x) - \ln(x_m - x) = rt +C $$

$$ \ln\left(\frac{x}{x_m -x}\right) = rt +C $$

两边同时取指数

$$ \frac{x}{x_m -x} = e^{rt +C} $$

得到 $$ x = \frac{x_m}{1+C_1e^{-rt}} $$

带入初始条件 $$ x(0) = \frac{x_m}{1+C_1} = x_0 $$ 得到 $$ C_1 = \frac{x_{m}}{x_{0}}-1 $$

最终得到 $$ x(t)=\frac{x_{m}}{1+\left(\frac{x_{m}}{x_{0}}-1\right) e^{-r t}} $$

我们可以看一下 logistic模型的增长曲线

## 阻滞增长模型的绘图

## 输入初始参数
x0 = 1
xm = 30
r = 0.02

# 获得x 和 f(x)
t = [i for i in range(400)]
x_t = [xm * (1 + (xm/x0 -1)*np.exp(- r * time))**(-1) for time in t] 

# 绘图1:增长曲线
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,x_t)
plt.xlabel("Time")
plt.ylabel("Population")

# 绘图2:增长速度
plt.subplot(1,2,2)
x = [xm*i/30 for i in range(31) ]
deri_x = [r * xx * (1 - xx/xm) for xx in x]
plt.plot(x,deri_x)
plt.xlabel(r'$x$')
plt.ylabel(r'$\frac{\mathrm{d}x}{\mathrm{d}t}$')

Text(0, 0.5, '$\\frac{\\mathrm{d}x}{\\mathrm{d}t}$')

上面的阻滞增长模型，是荷兰生物数学家Verhulst 19世纪中叶提出的．它不仅能够大体上描述人口及许多物种数量（如森林中的树木、鱼塘中的鱼群等）的变化规律，而且在社会经济领域也有广泛的应用，例如耐用消费品的销售就可以用它来描述。基于这个模型能够描述一些事物符合逻辑的客观规律，人们常称它为logistic模型。

模型的参数估计、检验和预报

用阻滞增长模型进行人口预报，先要作参数估计．除了初始人口$x_0$外，还要估计$r$和$x_m$。它们可以用人口统计数据拟合得到，也可以辅之以专家的估计。

原方程可以写为 $$ \frac{\frac{\mathrm{d} x}{\mathrm{d} t}}{x}=r-s x, s=\frac{r}{x_{m}} $$

上式左端可以从实际人口数据用数值微分算出，右端对参数$r$,$s$是线性的,可借助最小二乘法获得。（什么是数值微分:差分）

## 对比一个世纪后的预测值和真实值
import seaborn as sns
from sklearn.linear_model import LinearRegression

# 输入初始数据
year = [i for i in range(1790,2010,10)]
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,
50.2,62.9,76,92,106.5,123.2,131.7,150.7,179.3,204,226.5,251.4,281.4]
y = []
for i in range(len(population)-1):
    y.append( (population[i+1] - population[i])/10 / population[i]  )
plt.figure(figsize = (15,5))

# 可视化线性回归结果
plt.subplot(1,2,1)
sns.regplot(population[1:],y,ci = 0)
plt.ylabel(r'$\frac{\mathrm{d}x}{\mathrm{d}t}/{x}$')
plt.xlabel(r'$x$')
lrModel = LinearRegression()
lrModel.fit(np.array(population[7:-1]).reshape(-1,1),y[6:-1])
r = lrModel.intercept_
xm = r/(- lrModel.coef_)

# 可视化拟合结果
plt.subplot(1,2,2)
x0 = population[0]
plt.plot(year,population,label = 'Real Data')
pop_predicted = [xm * (1 + (xm/x0 -1)*np.exp(- r * (time - 1790)))**(-1) for time in year]
plt.plot(year,pop_predicted,label = 'Predicted')
plt.xlabel('Year')
plt.xlabel('Population')
plt.legend()

D:\software_install\anaconda\lib\site-packages\seaborn\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(

<matplotlib.legend.Legend at 0x18b5e298e20>

可以看出，这个模型整体拟合效果不错。我们可以用模型计算2000年的人口，与已知的实际数据比较，来检验模型是否合适。

除此之外，我们还可以考虑使用之前讲过的非线性最小二乘拟合进行人口的预测。

## 复习前面学到的预测模型中的非线性拟合
plt.scatter(year,population,label = 'Real Data')
def logistic(t,xm,x0,r):
    return xm * (1 + (xm/x0 -1)*np.exp(- r * (t-1790) ))**(-1)
from scipy.optimize import curve_fit 
#popt1, pcov1 = curve_fit(logistic, year, population,p0 = [1,0,0]) 
popt1, pcov1 = curve_fit(logistic, year, population,p0 = [1,1,1]) 

#popt数组中，三个值分别是待求参数a,b,c  
y2 = [logistic(i, popt1[0],popt1[1],popt1[2]) for i in year]
plt.plot(year,y2,'b--',label = 'Fitting Curve')  
plt.legend()

<ipython-input-10-4343891452f1>:4: RuntimeWarning: overflow encountered in exp
  return xm * (1 + (xm/x0 -1)*np.exp(- r * (t-1790) ))**(-1)

<matplotlib.legend.Legend at 0x18b5e4cd820>

可以看到，非线性拟合的误差更小，请解释为什么？

Questions???

参考文献

• ThomsonRen 的 github https://github.com/ThomsonRen/mathmodels
• youcans Python小白的数学建模课 https://www.zhihu.com/column/c_1381900867424075776
• 司守奎 的 《Python数学实验与建模》