题目来源姜启源《数学模型 第四版》。

意在练习python的matplotlib库和使用latex进行数学公式编辑。

数学模型(一) 如何施救药物中毒

1 | 题目描述

已知条件:

• 孩子一口气服下$1100 \mathrm{mg}$药物
• 孩子有2000ml的血液
• 血药浓度达到$100\mu \mathrm{g}/\mathrm{ml}$时将出现严重中毒, 达到$200\mu \mathrm{g}/\mathrm{ml}$将致命
• 药物吸收的半衰期为5h(胃肠道浓度), 排除的半衰期为6h(血液浓度)

2 | 模型建立

记:

$$\mathrm{x}\left( \mathrm{t} \right) : \text{胃肠道中药量} \\ \mathrm{y}\left( \mathrm{t} \right) : \text{血液系统中药量}$$

做出如下假设:

1. 胃肠道中药物向血液系统的转移率与$\mathrm{x}\left( \mathrm{t} \right) $成正比, 比例系数记为$\lambda $, 总剂量1100mg的药物在t=0瞬间进入胃肠道
2. 血液系统中药物排除率与药量$\mathrm{y}\left( \mathrm{t} \right) $成正比, 比例系数记为$\mu $, t=0时血液中无药物

根据假设1, $\mathrm{x}\left( \mathrm{t} \right) $满足微分方程:

$$\frac{\mathrm{dx}}{\mathrm{dt}}=-\lambda \mathrm{x},\mathrm{x}\left( 0 \right) =1100 \tag{1}$$

根据假设2, $\mathrm{y}\left( \mathrm{t} \right) $满足微分方程:

$$\frac{\mathrm{dy}}{\mathrm{dt}}=\lambda \mathrm{x}-\mu \mathrm{y},\mathrm{y}\left( 0 \right) =0 \tag{2}$$

其中参数$\mu $和参数$\lambda $可以由半衰期确定

3 | 模型求解

微分方程(1)是可分离变量方程, 易得

$$\mathrm{x}\left( \mathrm{t} \right) =1100\mathrm{e}^{-\lambda \mathrm{t}} \tag{3}$$

利用药物吸收的半衰期为5h的条件可知:

$$\mathrm{x}\left( 5 \right) =1100\mathrm{e}^{-5\lambda}=\frac{\mathrm{x}\left( 0 \right)}{2}=550 \\ \Rightarrow \lambda =\frac{\ln 2}{5}=0.138629 \tag{4}$$

将方程(3)带入方程(2)得到一阶线性微分方程, 求解得:

$$\mathrm{y}\left( \mathrm{t} \right) =\frac{1100\lambda}{\lambda -\mu}\left( \mathrm{e}^{-\mu \mathrm{t}}-\mathrm{e}^{-\lambda \mathrm{t}} \right) \tag{5}$$

表明血液中的药量$\mathrm{y}\left( \mathrm{t} \right) $随着t先增后减并趋于0

为了根据药物排除的半衰期6h来确定$\mu $, 考虑某时刻$\tau $有$\mathrm{y}\left( \tau \right) =\mathrm{a}$, 则

$$\mathrm{y}\left( \mathrm{t} \right) =\mathrm{ae}^{-\mu \left( \mathrm{t}-\tau \right)},\mathrm{t}\geqslant \tau$$

由$\mathrm{y}\left( \tau +6 \right) =\frac{\mathrm{a}}{2}$得$\mu =\frac{\ln 2}{6}=0.115525$

将$\lambda $和$\mu $的值带入(3) (5)得

$$\mathrm{x}\left( \mathrm{t} \right) =1100\mathrm{e}^{-0.1386\mathrm{t}} \tag{6}$$ $$\mathrm{y}\left( \mathrm{t} \right) =6600\left( \mathrm{e}^{-0.1155\mathrm{t}}-\mathrm{e}^{-0.1386\mathrm{t}} \right) \tag{7}$$

4 | 作图与分析

这里使用python作图

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 

x = np.arange(0, 20, 0.1)
y1 = 1100 * np.exp(-0.1386 * x)
y2 = 6600 * (np.exp(-0.1155 * x) - np.exp(-0.1386 * x))

%matplotlib inline
plt.plot(x, y1, linestyle="-", label='x(t)')
plt.plot(x, y2, linestyle="-", label='y(t)')
plt.xlabel("t/h")
plt.ylabel("x, y/mg")
plt.legend(loc='upper right')
plt.title("胃肠道中药量x(t)和血液系统中药量y(t)")

由图可以看出若不及时抢救, 在5h后血液中药物浓度将达到400mg

由公式(6)易得y(2)=236.5mg, 计算药量达到400mg的时间:

通过解非线性方程:

$$6600\left( \mathrm{e}^{-0.1155\mathrm{t}}-\mathrm{e}^{-0.1386\mathrm{t}} \right)=400 \tag{7}$$

from scipy.optimize import fsolve
def f1(x):
    return [6600 * (np.exp(-0.1155 * x[0]) - np.exp(-0.1386 * x[0])) - 400]
print(fsolve(f1, [0]))

解得t = 4.87h

同时可以求出当t = 7.89h时血液中药量达到最大, 为442.1mg

def f2(x):
    a = 0.1386
    b = 0.1155
    return [1100 * a * np.exp(-1 * a * x[0]) - 1100 * a * b * (np.exp(-1 * b * x[0]) - np.exp(-1 * a * x[0])) / (a - b)]
print(fsolve(f2, [0]))
print()

def f3(x):
    return 6600 * (np.exp(-0.1155 * x) - np.exp(-0.1386 * x))
print(f3(7.89))