第四章 微分方程组模型

微分方程是包含连续变化的自变量、未知函数及其变化率（函数的微分或导数）的方程式，当我们的研究对象涉及某个过程或物体随时间连续变化的规律时，通常会建立微分方程模型。

虽然动态过程的变化规律一般要用微分方程建立的动态模型来描述，但是对于某些实际问题，建模的主要目的并不是要寻求动态过程每个瞬时的性态，而是研究某种意义下稳定状态的特征，特别是当时间充分长以后动态过程的变化趋势。这时常常不需要求解微分方程（并且我们将会看到，即使对于不太复杂的方程，解析解也不是总能得到的)，而可以利用微分方程稳定性理论，直接研究平衡状态的稳定性就行了。

指数增长模型

指数模型是一个早就用于描述生物群体增长的简单模型。 该模型假设在研究的时间范围内，只有生殖现象而没有死亡现象，而且生物群体可以获得无限的生长条件。

推导

美国人口统计.xlsx

已知今年人口为x_0，年增长率为r，下面的公式就能退出k年后的人口x_k：

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

这个公式的前提是年增长率k保持不变。

当考察一个国家或一个较大地区的人口随着时间延续而变化的规律时，为了利用微积分这一数学工具，可以将人口看作连续时间t的连续可微函数x(t).记初始时刻(t=0)的人口为x。假设单位时间人口增长率为常数r,rx(t)就是单位时间内x(t)的增量\frac{\mathrm dx}{\mathrm dt}，于是得到x(t)满足的微分方程和初始条件

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

解得：

x\left(t\right)=x_{0}\operatorname{e}^{rt}

称该式为指数增长模型，又称为Malthus人口模型。

那么如何对该模型的参数进行估计呢，使用最小二乘法，对上式取对数可得到：

y=rt+a,y=\ln x,a=\ln x_0

使用编程计算即可。

改进的指数增长模型

我们知道人口增长率不会一直为一个常数，根据经验，我们可以知道人口增长率会随着人口的增多而逐渐减小。

我们将r视为t的函数，假设r(t)=r_0-r_1t，并且将指数增长模型改写为：

\dfrac{\mathrm{d}x}{\mathrm{d}t}=r\left(t\right)x=\left(r_0-r_1t\right)x,x\left(0\right)=x_0

解得：

x{\left(t \right)} = x_{0} e^{t \left(r_{0} - \frac{r_{1} t}{2}\right)}

import sympy as sp

r0 = sp.Symbol('r_0')
r1 = sp.Symbol('r_1')
t = sp.Symbol('t')
x = sp.Function('x')(t)
x0 = sp.Symbol('x_0')

solve = sp.Eq(x.diff(t), (r0-r1*t)*x)
# 求解x(0)=x_0时的微分方程x(t)
sol = sp.dsolve(solve, x, ics={x.subs(t, 0): x0})

改进的指数增长模型图像

改进的指数增长模型考虑到增长率随时间的变化，计算结果有所改善，但是它没有反映增长率下降的成因，增长率函数形式也呈不确定性，给应用带来不便。为了建立更符合实际情况的人口增长模型，必须考察人口增长率下降的机理，修改原来的模型假设。

阻滞增长模型

在推出阻滞增长模型前，我们需要先推出指数增长模型。

阻滞增长模型（logistic模型）

资源和环境对人口增长的阻滞作用体现在对增长率，的影响上，使r随着人口数量x的增加而下降.将r表示为x的函数r(x),并且取既简单又便于应用的线性减函数r(x)=a+bx.为了赋予增长率函数r(x)中的系数a,b以实际含义，引入两个参数：

1. 内禀(bǐng)增长率r：r是x=0时的增长率，所以a=r。
2. 人口容量x_m：x_m是自愿和环境所能容纳的最大人口数量，当x=x_m时人口不再增长，即r(x_m)=r+bx+m=0。得到b=-\frac{r}{x_m}。

得到以下的人口增长模型：

\frac{\mathrm{d}x}{\mathrm{d}t}=rx\Big(1-\frac{x}{x_{m}}\Big),x(0)=x_{0}

求解得到模型的解为：

x{\left(t \right)} = \frac{x_{m}}{1 - \frac{\left(x_{0} - x_{m}\right) e^{- r t}}{x_{0}}}

上式称为logistic曲线。

import sympy as sp

r = sp.Symbol('r')
x_m = sp.Symbol('x_m')
t = sp.Symbol('t')
x_0 = sp.Symbol('x_0')

x = sp.Function('x')(t)

solve = sp.Eq(x.diff(t), r * x * (1 - x / x_m))
res = sp.dsolve(solve, x, ics={x.subs(t, 0): x_0})

利用所给数据进行最小二乘法拟合：

import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
# 绘图更改字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 读取数据
data = pd.read_excel('美国人口统计.xlsx', sheet_name='Sheet1')
data = data.to_dict('list')
x = [i[0] for i in data.values()]
del x[0]

t = np.arange(0, len(x), 1)
x = np.array(x)

def f(t, x_m, x_0, r):
    return x_m / (1-((x_0-x_m)*np.exp(-r*t))/x_0)

# 曲线拟合
popt, pcov = curve_fit(f, t, x)
x_m = popt[0]
x_0 = popt[1]
r = popt[2]
print(x_m, x_0, r)
# 绘制图像
yvals = f(t, x_m, x_0, r)
plot1 = plt.plot(t, x, 'o', label='原始数据')
plot2 = plt.plot(t, yvals, label='拟合曲线')
plt.xlabel('t')
plt.ylabel('x')
plt.legend(loc=4)
plt.title('美国人口统计')
plt.show()

得到三个参数分别为：

x_m = 486.1611~~ x_0 =8.1921~~ r = 0.2080

对所得数据进行相关性检验：

def get_indexes(y_predict, y_data):
    n = y_data.size

    # 计算 残差平方方差 SSE
    SSE = ((y_data - y_predict)**2).sum()

    # 均方差 MSE
    MSE = SSE / n

    # 均方根 RMSE 越接近 0，拟合效果越好
    RMSE = np.sqrt(MSE)

    # 求R方，0<=R<=1，越靠近1,拟合效果越好
    u = y_data.mean()
    SST = ((y_data - u)**2).sum()
    SSR = SST - SSE
    R_square = SSR / SST

    return SSE, MSE, RMSE, R_square

print(get_indexes(yvals, x)[3])

得到R^2=0.9975，拟合效果很好。

传染病模型

SI模型

SI模型这个模型将人群分为两类：易感染者(susceptible)和已感染者(infective)，分别简称健康人和患者.在传染病模型中取两个词的英文字头称为SI模型。设所考察地区的总人数（记作N）不变，时刻t健康人和患者在总人数中的比例分别记作s(t)和i(t),显然有s(t)+i(t)=1。

假设每个患者每天有效接触的人数是常数\lambda，称接触率，且当健康人被有效接触后立即被感染成为患者，\lambda也称感染率。

根据上述假设，每个患者每天有效接触的健康人数是入s(t),全部患者Ni每天有效接触的健康人数是Nλs(t)(t),这些健康人立即被感染成为患者，于是患者比例i(t)满足微分方程：

\frac{\mathrm{d}i}{\mathrm{d}t}=\lambda si

已知s=1-i，初始时刻患者比例记为i_0，得：

\dfrac{\mathrm{d}i}{\mathrm{d}t}=\lambda i(1-i),i(0)=i_0

这是标准的logistic模型。也就是说随着时间，感染人数的速度先快后慢，最后所有人都被感染。这不符合实际情况。究其原因，是SI模型只考虑了健康人可以被感染，而没有考虑到患者可以被治愈的情况.下面的模型将增加关于患者被治愈的假设。

SIS模型

SIS的假设条件是：患者每天被治愈的人数比例是常数\mu，称为治愈率.

于是我们得到了新方程：

\frac{\mathrm{d}i}{\mathrm{d}t}=\lambda si-\mu i

定义：\sigma=\lambda/\mu，将s=1-i代入上式可得：

\dfrac{\mathrm{d}i}{\mathrm{d}t}=\lambda i\Big[(1-1/\sigma)-i\Big]

当\sigma>1，上式还是logistic曲线。当\sigma\le1，曲线将单调下降。

由此可知，\sigma决定了患者比例是持续增加还是持续减少。

SIR模型

SIR模型许多传染病如天花、流感、麻疹等，治愈后免疫力很强，于是患者愈后不会成为易感染的健康人，当然也不再是患者.传染病模型中将病愈免疫后和因病死亡的人称为移除者(removed),由此将这个模型称为SIR模型.仍设总人数N不变，时刻t移除者在总人数中的比例记作r(t),有：

s(t)+i(t)+r(t)=1

得到一下方程组

\frac{\mathrm{d}i}{\mathrm{d}t}=\lambda si-\mu i

\frac{\mathrm{d}s}{\mathrm{d}t}=-\lambda si

\frac{\mathrm{d}r}{\mathrm{d}t}=\mu i

这个线性方程组没有解析解，只有数值解。利用编程计算可以得到图像如下：

\lambda = 0.6,\mu=0.3,s(0)=0.99,i(0)=0.01

\lambda = 0.5,\mu=0.4,s(0)=0.99,i(0)=0.01

# SIR模型绘图
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# 绘图更改字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

t = np.arange(0, 50, 0.1)
lambda_ = 0.5  # 感染率
mu = 0.4  # 移除率
s0 = 0.99  # 初始易感人群比例
i0 = 0.01  # 初始感染人群比例
r0 = 0  # 初始移除人群比例

def f(y, t):
    """SIR模型微分方程

    Args:
        y (list): [s, i, r], s:易感人群比例, i:感染人群比例, r:移除人群比例
        t (np.arange): 时间

    Returns:
        list: [ds/dt, di/dt, dr/dt]
    """
    s, i, r = y
    return [-lambda_ * s * i, lambda_ * s * i - mu * i, mu * i]

result = odeint(f, [s0, i0, r0], t)  # 求解微分方程
plt.plot(t, result[:, 0], label='s(t)')
plt.plot(t, result[:, 1], label='i(t)')
plt.plot(t, result[:, 2], label='r(t)')
plt.grid()
plt.title('SIR模型')
plt.legend()
plt.show()

参数时变的SIR模型

假设\lambda(t)和\mu(t)表示第t天的感染率和移除率，SIR门票销售可以表示为以下方程组：

\begin{aligned} \frac{\mathrm{d}s}{\mathrm{d}t}=-\lambda\left(t\right)s\left(t\right)i\left(t\right) \\ \frac{\mathrm{d}i}{\mathrm{d}t} =\lambda\left(t\right)s\left(t\right)i\left(t\right)-\mu\left(t\right)i(t) \\ \frac{\mathrm{d}r}{\mathrm{d}t}=\mu\left(t\right)i\left(t\right) \end{aligned}

引入不可控带菌者和疑似已感染者的模型

在以上的模型中再增加两个人群：不可控带菌者和疑似已感染者，记为c(t)和e(t)。以下符号分别代表：

5个人群之间的转移关系如图所示：

根据上图就可以建立起以下的模型：

\begin{gathered} s\left(t\right)+c\left(t\right)+e\left(t\right)+i\left(t\right)+r\left(t\right)=1 \\ \begin{aligned}\frac{\mathrm{d}s}{\mathrm{d}t}=\beta e-\lambda cs\end{aligned} \\ \frac{\mathrm{d}c}{\mathrm{d}t}=\lambda(1-\alpha)cs-\boldsymbol{\varepsilon}c \\ \frac{\mathrm{d}e}{\mathrm{d}t}=-\left(\beta+\gamma\right)e+\lambda\alpha cs \\ \begin{aligned}\frac{\mathrm{d}i}{\mathrm{d}t}&=\varepsilon c-\mu i+\gamma e\end{aligned} \\ \frac{\mathrm{d}r}{\mathrm{d}t}=\mu i \end{gathered}

引入尚未隔离和已经隔离已感染者的模型

结合下图来说明各个人群的构造的转移关系：

我们可以通过一下的方程组来建立模型：

\begin{gathered} \frac{\mathrm{d}E}{\mathrm{d}t}= \beta\left(t\right)\left(kE+I\right)-\left(\varepsilon+\lambda\right)E \\ \frac{\mathrm{d}I}{\mathrm{d}t}=\boldsymbol{\varepsilon}E-\left(\theta+\delta\right)I \\ \frac{\mathrm{d}Q}{\mathrm{d}t}=\lambda E-\sigma Q \\ \frac{\mathrm{d}J}{\mathrm{d}t}=\theta I+\sigma Q-\left(\gamma+\delta\right)J \\ \frac{\mathrm{d}R}{\mathrm{d}t}=\gamma J \end{gathered}

差分方程模型

例题 管住嘴迈开腿

**问题提出：**如何科学的制订减肥计划？

1.体重增加正比于吸收的热量，平均每8000kcal(kcal为非国际单位制单位，(1kcal=4.2kJ)增加体重1kg;
2.身体正常代谢引起的体重减少正比于体重，每周每千克体重消耗热量一般在200kcal至320kcal之间，且因人而异，这相当于体重70kg的人每天消耗2000kcal至3200 kcal;
3.运动引起的体重减少正比于体重，且与运动形式和运动时间有关；
4.为了安全与健康，每周吸收热量不要小于10000kcal,且每周减少量不要超过1000kcal,每周体重减少不要超过1.5kg。

基本模型

基本模型记第k周（初）体重为w(k)(kg),第k周吸收热量为c(k)(kcl),k=1,2,…。设热量转换（体重的）系数为α，身体代谢消耗系数为\beta，根据模型假设，正常情况下（不考虑运动）体重变化的基本方程为：；

w(k+1)=w(k)+\alpha c(k)-\beta w(k),~~k=1,2,...

当增加运动时，只需将\beta改为\beta+\beta_1，\beta_1由运动的形式和时间决定。

具体的计划

假设某人身高1.70m,体重100kg,BMI高达34.6.自述目前每周吸收20000kcal热量，体重长期未变.试为他按照以下方式制订减肥计划，使其体重减至75kg(即BMI=26)并维持下去：

1. 在正常代谢情况下安排一个两阶段计划，第一阶段：吸收热量由20000kcal每周减少1000kcal,直至达到安全下限（每周10000kcal)；第二阶段：每周吸收热量保持下限，直至达到减肥目标。

2. 为加快进程而增加运动，重新安排两阶段计划。

3. 给出达到目标后维持体重不变的方案。

首先我们应该确定某人的代谢消耗系数\beta，由c=20000kcal热量，体重保持不变，由式：w(k+1)=w(k)=w，c(k)=c得：

w=w+\alpha c-\beta w

可得\beta=0.025。假设\alpha=\frac 1 {8000}kg/kcal。

由第一阶段的计划：

c(k)=20000-1 000k, k=1,2,...,10

将其代入差分方程：

\begin{aligned} w(k+1)&=(1-\beta)w(k)+\alpha (20000-1000k)\\&=0.975w(k)+2.5-0.125k,~k=1,2,...,10 \end{aligned}

利用Python，求得第11周（初）的体重w(11)。

# 差分方程
def w_k(alpha, beta, times):
    if times == 1:
        return 100  # w(1)=100
    times -= 1
    return (1-beta)*w_k(alpha, beta, times) + alpha * (20000-1000*times)

if __name__ == '__main__':
    print(w_k(1/8000, 0.025, 11))

解得：w(11)=93.6157。

第二阶段要求吸收热量保持下限c_{\min}，类推可得：

w(k+1)=(1-\beta)w(k)+\alpha c_{\min}

编程计算：

def w_k_in_2(alpha, beta):
    res = 93.6157
    k = 11
    while res > 75:
        res = (1-beta)*res + alpha*10000
        k += 1
    return {k: res}

if __name__ == '__main__':
    print(w_k_in_2(1/8000, 0.025))

解得w(33)=74.9887。所以在不运动的情况下能够在32周减到74kg。

为了加快进程而增加运动，设每周运动时间为th，热量消耗为\gamma，将之前的式子中\beta改为\beta+\alpha\gamma t。试取\alpha\gamma t=0.005，则继续用编程法可得到：w(11)=89.3318，w(23)=74.7388。

此外，我们可以将两次减肥的情况做成曲线图，可得：

# 引入绘图
import matplotlib.pyplot as plt
import numpy as np

# 绘图更改字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 差分方程

def w_k(alpha, beta, times, res):
    if times == 1:
        return 100  # w(1)=100
    times -= 1
    weight = (1-beta)*w_k(alpha, beta, times, res) + alpha * (20000-1000*times)
    res[times+1] = weight
    return weight

def w_k_in_2(alpha, beta, weight):
    res = {}
    k = 11
    while weight > 75:
        weight = (1-beta)*weight + alpha*10000
        k += 1
        res[k] = weight
    return res

if __name__ == '__main__':
    res1 = {1: 100}
    w_11 = w_k(1/8000, 0.025, 11, res1)
    res1.update(w_k_in_2(1/8000, 0.025, w_11))
    print(res1)
    x1 = [i[0] for i in res1.items()]
    y1 = [i[1] for i in res1.items()]
    res2 = {1: 100}
    w_11 = w_k(1/8000, 0.03, 11, res2)
    res2.update(w_k_in_2(1/8000, 0.03, w_11))
    print(res2)
    x2 = [i[0] for i in res2.items()]
    y2 = [i[1] for i in res2.items()]
    # 绘制曲线，并标注曲线代表的内容
    plt.plot(x1, y1, label='正常代谢')
    plt.plot(x2, y2, label='增加运动')
    plt.annotate('正常代谢', xy=(x1[15], y1[15]), xytext=(x1[15], y1[15]+1))
    plt.annotate('增加运动', xy=(x2[15], y2[15]), xytext=(x2[15]-5, y2[15]))
    # 添加网格
    plt.grid(True)
    plt.xlabel('周数')
    plt.ylabel('体重')
    plt.show()

我们还需要检查每周的下降体重是不是超过了1.5kg。

def max_weight_loss(res):
    m_weight_loss = 0
    l = [i for i in res.items()]
    for i in range(1, len(l)):
        weight_loss = l[i][1] - l[i-1][1]
        if weight_loss < m_weight_loss:
            m_weight_loss = weight_loss
    return m_weight_loss

if __name__ == '__main__':
		res1 = {1: 100}
    w_11 = w_k(1/8000, 0.025, 11, res1)
    res1.update(w_k_in_2(1/8000, 0.025, w_11))
    print(max_weight_loss(res1))
    res2 = {1: 100}
    w_11 = w_k(1/8000, 0.03, 11, res2)
    res2.update(w_k_in_2(1/8000, 0.03, w_11))
    print(max_weight_loss(res2))

第一种方法失去的体重最大为1.1183，第二种方法失去的体重最大为1.4742，均未超过1.5。

达到目标后维持体重不变的方案，最简单的是寻求每周吸收热量保持某常数值c，使体重w=75kg不变。由式可得：

w=w+\alpha c-(\beta+\alpha\gamma t)w

求解得到：

c=\frac{(\beta+\alpha\gamma t)w}{\alpha}

由该式，在正常代谢下（\gamma=0）可以算出c=15000kcal，若增加\gamma t=40的运动，则c=18000kcal。

import sympy as sp

w = sp.symbols("w")
alpha = sp.symbols("alpha")
c = sp.symbols("c")
beta = sp.symbols("beta")
gamma = sp.symbols("gamma")
t = sp.symbols("t")

exp = alpha*c-(beta+alpha*gamma*t
               )*w
# 求解c
c = sp.solve(exp, c)[0]
c
# 正常代谢
c_1 = sp.solve(
    exp.subs({w: 75, alpha: 1/8000, beta: 0.025, gamma: 1, t: 0}), dict=True)
c_1
# 增加运动
c_2 = sp.solve(
    exp.subs({w: 75, alpha: 1/8000, beta: 0.025, gamma: 1, t: 40}), dict=True)
c_2

如果某人的摄入量一开始就一直保持c=15000kcal，那么的体重会随着时间推移逐渐接近到75kg。

体重随时间变化曲线：

import sympy as sp

beta, alpha, c = sp.symbols('beta, alpha, c', positive=True)
w = sp.Function('w')
w_star = sp.symbols('w^*')
n = sp.symbols('n', positive=True, integer=True)

s = sp.summation((1-beta)**k, (k, 0, n - 1))

exp = (1 - beta) ** n * w(1) + alpha * c * s
sp.simplify(exp)
# 构造方程
solve = sp.Eq(exp, w_star)
res_n = sp.solve(solve, n)[0]  # n的结果, n表示需要第几周
# 绘制不同摄入的曲线
res_w_star = sp.solve(solve, w_star)[0].subs(
    {alpha: 1/8000, beta: 0.025, w(1): 100})

p = sp.plot((res_w_star.subs({c: 15000}), (n, 0, 80)),
            legend=True, label='c=15000', show=False, axis_center=(0, 65))
p.extend(sp.plot((res_w_star.subs({c: 14000}),
         (n, 0, 80)), legend=True, label='c=14000'))
p.extend(sp.plot((res_w_star.subs({c: 13000}),
         (n, 0, 80)), legend=True, label='c=13000'))
p.extend(sp.plot((res_w_star.subs({c: 12000}),
         (n, 0, 80)), legend=True, label='c=12000'))
p.show()