微分方程建模实例（1）

更新于
人口预报模型。

例：利用下表给出的近两个世纪的美国人口统计数据（以百万为单位），建立人口预测模型，用它估计2010年美国的人口。

年份 1790 1800 1810 1820 1830 1840 1850 1860
人口 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4
年份 1870 1880 1890 1900 1910 1920 1930 1940
人口 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7
年份 1950 1960 1970 1980 1990 2000
人口 150.7 179.3 204.0 226.5 251.4 281.4

解：记x(t)为第t年的人口数量，设人口年增长率r(x)为x的线性函数，r(x)=r-sx。自然资源和环境条件能容纳的最大人口数为xm，即当x=xm时，增长率r(xm)=0，可得
$$r(x) = r(1-\frac{x}{x_{m}})$$
建立Logistic人口模型
$$
\begin{cases}
\frac{dx}{dt} = r(1-\frac{x}{x_{m}})x \\
x(t_{0})=x_{0}
\end{cases}
$$
其解为
$$
x(t)=\frac{x_{m}}{1+(\frac{x_{m}}{x_{0}}-1)e^{-r(t-t_{0})}}
$$
接下来进行参数估计，利用提供的数据拟合xm和r

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import numpy as np
from scipy.optimize import curve_fit
a=[]
b=[]
with open("data.txt") as f:    # 打开文件并绑定对象f
    s=f.read().splitlines()  # 返回每一行的数据
    '''
    1790	1800	1810	1820	1830	1840	1850	1860
    3.9	5.3	7.2	9.6	12.9	17.1	23.2	31.4
    1870	1880	1890	1900	1910	1920	1930	1940
    38.6	50.2	62.9	76.0	92.0	106.5	123.2	131.7
    1950	1960	1970	1980	1990	2000		
    150.7	179.3	204.0	226.5	251.4	281.4
    '''
for i in range(0, len(s),2):  # 读入奇数行数据
    d1=s[i].split("\t")
    for j in range(len(d1)):
        if d1[j]!="": a.append(eval(d1[j]))  # 把非空的字符串转换为年代数据
for i in range(1, len(s), 2):  # 读入偶数行数据
    d2=s[i].split("\t")
    for j in range(len(d2)):
        if d2[j] != "": b.append(eval(d2[j])) # 把非空的字符串转换为人口数据
# c=np.vstack((a,b))
# np.savetxt("result.txt", c)  # 数据保存
x=lambda t, r, xm: xm/(1+(xm/3.9-1)*np.exp(-r*(t-1790)))
bd=((0, 200), (0.1,1000))  # 约束两个参数的下界和上界
popt, pcov=curve_fit(x, a[1:], b[1:], bounds=bd)
print(popt)
print("2010年的预测值为：", x(2010, *popt))
0%