微分方程人口模型以及Matlab代码实现

马尔萨斯模型

模型假设

1.设x(t)表示t时刻的人口数,且x(t)连续可微。
2.人口的增长率r是常数(增长率=出生率-死亡率)。
3.人口数量的变化是封闭的,即没有人口迁移，且人口数量的增加与减少只取决于人口中个体的生育和死亡,且每一个体都具有同样的生育能力与死亡率。

建模
$$\begin{gathered} \frac{dx(t)}{dt}=\underset{\Delta t\to 0}{\lim }\frac{x(t+\Delta t)-x(t)}{\Delta t}=\underset{\Delta t\to 0}{\lim }\frac{x(t)\cdot (1+r{)}^{\Delta t}-x(t)}{\Delta t} \\ =x(t)\underset{\Delta t\to 0}{\lim }\frac{(1+r{)}^{\Delta t}-1}{\Delta t}=x(t)\underset{\Delta t\to 0}{\lim }\frac{e^{ln(1+r)\cdot \Delta t}-1}{\Delta t}=x(t)ln(1+r)=x(t)r \end{gathered}$$
注：因为增长率通常很小，可以近似 $\ln (1+r)\approx r$

代码演示
假设x(0)=1000,r=1%，查看50个单位时间内人口变化情况

r=0.01
syms x(t)
eqn=(diff(x,t)==r*x)
cond=(x(0)==1000)
f=dsolve(eqn,cond)
figure(1)
fplot(f,[0,50],'-*r')

阻滞增长模型（Logistic模型）

阻滞增长模型通常用于描述在资源有限的环境中，种群或某种量的增长受到约束的情况。由于资源是有限的，模型假设随着人口增加，增长率呈现线性下降，当人口达到环境容纳量K时，增长率为0，人口数量达到平衡或开始出现负增长。

建模

$$\begin{gathered} r(x)=r_0-sx \\ r(K)=r_0-sK=0\Rightarrow s=\frac{r_0}{K} \\ \frac{dx(t)}{dt}=r\cdot x=r_0\cdot (1-\frac{x}{K})\cdot x \end{gathered}$$

求出微分方程x(t),求二阶导后可以发现规律：x=K/2对应的t为拐点，在这之前增长率保持上升，之后增长率下降直至0。

代码演示
初始增长率5%，x(0)=1000,K=5000,查看200个单位时间内人口数量

%阻滞增长模型
r0=0.05
x_0=1000
K=5000
syms x(t)
eqn=(diff(x,t)==r0*(1-x/K)*x)
cond=(x(0)==x_0)
f2=dsolve(eqn,cond)
figure(1)
fplot(f2,[0,200],'-*r')

实战案例

题目

解答
①求解logistic模型一般解

%阻滞增长模型
x0=3.9
syms x(t)  K r0
eqn=(diff(x,t)==r0*(1-x/K)*x)
cond=(x(1790)==x0)
f=dsolve(eqn,cond)

$$

f =K/(exp(1790*r0 + log((10*K)/39 - 1) - r0*t) + 1)
 

②输入数据，打开cftoo拟合l工具箱

year=1790:10:2000
population=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4
cftool

选择数据->自定义方程->输入方程->打开高级选项，设置合理的初值和范围->得到满意的拟合曲线->导出代码

③得到拟合的K和r0之后，进行后续人口预测

[fitresult, gof] = createFit(year, population)
K=342.45
r0=0.0274
year_pred=2001:2031
population_pred=K./(exp(1790*r0 + log((10*K)/39 - 1) - r0.*year_pred) + 1)
figure(1)
plot(year,population,'*',year_pred,population_pred)