3. 微分方程模型

数学建模

张瑞
中国科学技术大学数学科学学院

微分方程模型

常微分方程

此类模型的建模思想是:

  • 模型的因素随时间、空间而变化,可看作是关于随时间、空间的连续可微函数;
  • 因素之间的联系可用微分方程来表示;
  • 微分方程的求解是解决实际问题的核心。

常见方法：

  1. 按规律直接列方程。许多现象所满足的规律已知知道，如：牛顿第二定律
  2. 微元分析，或区域积分
  3. 模拟近似

一阶微分方程

(1)
\[\frac{dy}{dt}=f(t,y) \]

一阶微分方程组

(2)
\[\left\{\begin{array}{lcl} y'_1 &=& f_1(t,y_1,y_2,\cdots,y_m)\\ y'_2 &=& f_2(t,y_1,y_2,\cdots,y_m)\\ &\vdots& \\ y'_m &=& f_m(t,y_1,y_2,\cdots,y_m) \end{array}\right. \]

高阶微分方程

(3)
\[y^{(n)}=f(t,y,y',y'',\cdots,y^{(n-1)}) \]

假设给定高阶微分方程

(4)
\[y^{(n)}=f(t,y,y',y'',\cdots,y^{(n-1)}) \nonumber \]

我们令

\[y_1=y, y_2=y', y_3=y'', \cdots, y_n=y^{(n-1)} \]

则引进的新变量$y_1,y_2,\cdots,y_n$ 满足如下一阶微分方程组

(5)
\[\left\{\begin{array}{lcl} y'_1 &=& y_2\\ y'_2 &=& y_3\\ &\vdots& \\ y'_n &=& f(t,y_1,y_2,\cdots,y_n) \end{array}\right. \nonumber \]

考虑以下初值问题

(6)
\[\left\{\begin{array}{l} \dfrac{dy}{dt}=f(t,y) \\ y(t_0)=y_0 \end{array}\right. \]

这里$y$是$t$的未知函数，我们希望从(6)中给出的信息去构造它。

  • 第一个等式给出$y(t)$在任意$t$处的斜率$f(y,t)$，其中$f$是指定的；
  • 第二个等式给出$y(t)$在$t_0$处的一个特定值。

例如

(7)
\[\left\{\begin{array}{l} y'=(t+\sin y)^2 \\ y(0)=3 \end{array}\right. \nonumber \]

传染病的传播

SARS(Severe Acute Respiratory Syndrome,严重急性呼吸道综合症,俗称:非典型肺炎)是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响,我们从中得到了许多重要的经验和教训,认识到定量地研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。请你们对SARS的传播建立数学模型。

困难

  • 传染病传播涉及的因素很多，如已感染病人的多少，易受传染者的多少，传染率的大小，潜伏期的长短，传染途径如何。
  • 如果还要考虑人群的流动情况、疾病预防与卫生医疗条件等影响，那么传染病的传播就变得非常复杂。

简化

  • 一开始就把众多的因素考虑在内将会使建模者陷入困境而无所适从，倒不如抓住主要因素，简化问题，建立初步的数学模型。
  • 将初步模型所得结果与实际情况分析比较，找出不足之处，修改原有假设，逐步建立与实际相吻合的模型。

假定

  1. 每个病人在单位时间内接触$\alpha_0$个人，而且立即被感染
  2. 每个病人一直是病人（不考虑被治愈的情况）

模型一 设时刻$t$时的病人数为$i(t)$，最初的染病者为$i(0)=i_0$，则在$\Delta t$时间内增加的病人数为

\[i(t+\Delta t)-i(t)=\alpha_0 i(t) \Delta t, \]

于是得到微分方程

\[\left\{\begin{array}{l} \dfrac{d i(t)}{dt}=\alpha_0 i(t) \\ i(0)=i_0 \end{array}\right. \label{infection-model-one} \]

其解为 $i(t)=i_0 e^{\alpha_0 t}$.

  • 模型结果表明，传染病的传播是按指数函数发展的。
  • 这个结果可以推出当$t\rightarrow\infty$时$i(t)\rightarrow\infty$。
  • 如果发生了这种传染病，则要使$\alpha_0=0$，也就是完全隔离。

这个模型可能与传染病传播的初期比较吻合,但长期看显然是不合实际情况的。

问题在于假设不合理，特别是假设1（每个病人单位时间内传染的人数是常数）。

  • 在传播初期，传染病人少未被传染者众多，
  • 而在传播中后期，传染病人增多未被传染者逐渐减少，因而在不同时期的传染率是不同的。

为了吻合实际情况，我们要修改原有假设并建立新的模型。

假设

  1. 每个病人在单位时间内传染的人数与此时未被传染的人数成正比，
  2. 每个病人一直是病人（不考虑被治愈的情况）

模型二 $i(t),s(t)$分别表示已感染病人数和未被传染人数，总人数为$N=i(t)+s(t)$; 基于以上假设，每个病人单位时间内可传染人数为$\alpha s(t)$，则可建立微分方程

\[\left\{\begin{array}{l} \dfrac{d i(t)}{dt}=\alpha s(t) i(t) \\ s(t)+i(t)=N \\ i(0)=i_0 \end{array}\right. \label{infection-model-two} \]

由分离变量法解得 $i(t)=\dfrac{N}{1+(N/i_0-1)e^{-\alpha N t}}$.

mm3-ex-spread mm3-ex-spread-dt

由

\[\dfrac{d i(t)}{dt}=\dfrac{\alpha N^2(N/i_0-1)e^{-\alpha Nt}}{[1+(N/i_0-1)e^{-\alpha N t}]^2} \]

可知，当$t_1=\dfrac{\ln(N/i_0-1)}{\alpha N}$时，$\dfrac{di(t)}{dt}$的达到极大值。

  • 当传染病的传染强度$\alpha$或总人数$N$增加时，$t_1$都将变小，即传染病的高峰来得快，这与实际情况相吻合。
  • 同时，如果知道了传染病强度（$\alpha$可由统计数据得出），即可预报传染病高峰到来的时间$t_1$,这对于传染病的防治工作是有益处的。

不足

当$t\rightarrow\infty$时，$i(t)\rightarrow N$,即最后人人都要得病。这显然不合实际情况，造成不合理的原因是假设了“人得病后经久不愈”。

假设

  1. 假设患过传染病而痊愈的人具有较长期免疫力，不考虑反复受传染的情形。
  2. 另外设传染病的潜伏期很短，可忽略不计，也就是说一个人患了病后立即成为传染者。

于是我们把居民分为三类：

  • [第一类] 由能够把疾病传染给别人的那些染病者组成，用$i(t)$表示$t$时刻的第一类人数。
  • [第二类] 由未被传染但可能得病成为传染者的那些人组成，用$s(t)$表示$t$时刻的第二类人数。
  • [第三类] 包括患病死去的人，病愈后具有长期免疫力的人，以及被隔离起来的人，用$r(t)$表示$t$时刻的第三类人数。

模型三（SIR）

  1. 易受传染者人数$s(t)$的变化率正比于第一类人数$i(t)$与第二类人数$s(t)$的乘积；
  2. 由第一类$i(t)$向第三类$r(t)$转变的速率与第一类人数成正比。
\[\left\{\begin{array}{l} \dfrac{d s(t)}{dt}=-\alpha s(t) i(t) \\ \dfrac{d i(t)}{dt}=\alpha s(t) i(t)-\beta i(t) \\ \dfrac{d r(t)}{dt}=\beta i(t) \end{array}\right. \]

其中$\alpha$是传染率，$\beta$是排除率(治愈+死亡)，两者皆为比率常数。

mm3-spread-sir

$\alpha=0.6$, $\beta=0.1$

三个方程相加，可以得到

\[\dfrac{d}{dt}\Big[s(t)+i(t)+r(t)\Big]=0, \]

即总人数$s(t)+i(t)+r(t)$保持不变

由方程组的前两个方程

\[\left\{\begin{array}{l} \dfrac{d s(t)}{dt}=-\alpha s(t) i(t) \\ \dfrac{d i(t)}{dt}=\alpha s(t) i(t)-\beta i(t) \end{array}\right. \label{first-two-eqns} \]

可得

\[\dfrac{d i}{d s}=\dfrac{di(t)}{dt}/\dfrac{ds(t)}{dt}= \dfrac{\alpha s i - \beta i}{-\alpha s i}=-1+\dfrac{\beta}{\alpha s}, \]

解出$i$关于$s$的函数关系

\[i(s)=-s+\dfrac{\beta}{\alpha}\ln s + C. \]

若$t=0$时，$i(0)=i_0$, $s(0)=s_0$，则有

\[i(s)=(s_0+t_0)-s+\dfrac{\beta}{\alpha}\ln \frac{s}{s_0} \]

记$\rho=\dfrac{\beta}{\alpha}$, 则有

\[i'(s)=-1+\dfrac{\rho}{s} \ \left\{\begin{array}{l} <0,\ \mathrm{when} \ s>\rho \\ =0,\ \mathrm{when} \ s=\rho \\ >0,\ \mathrm{when} \ s<\rho \end{array}\right. \nonumber \]
  • 如果$s_0<\rho$, $i'(s)>0$，随着$s(t)$的减少，$i(t)$也减小。当$i(t)$单调减小到$0$, $\frac{di}{dt}=0$, $\frac{ds}{dt}=0$
  • 如果$s_0>\rho$, 则随着$s(t)$减小到$\rho$,$i(t)$增加，且当$s=\rho$时$i(t)$达到最大值，而当$s(t)<\rho$时$i(t)$开始减小。此时，$i(t)$会先增加，然后减少。

由以上分析可得出结论：只当居民中易受传染者的人数超过阀值时，即$s>\rho$,传染病才会蔓延。

我们可以用常识来检验上述结论是否符合实际。

  • 当人口密度高，卫生医疗条件不好，隔离不良导致排除率低时，阀值$\rho$较小，传染病就会很快蔓延；
  • 反之，人口密度低，社会条件好，有良好的公共卫生医疗条件而排除率高时，阀值$\rho$较大，传染病在有限范围内出现便会很快被扑灭。

预防手段

  1. 接种疫苗，减少易感人数$s$
  2. 提高医疗技术，增加治愈率$\beta$。提高卫生条件水平，降低感染率$\alpha$

mm3-spread-sir mm3-spread-sir-2

左：$\alpha=0.6$, $\beta=0.1$ 右：$\alpha=0.6$, $\beta=0.5$

常微分方程求解

变量分离法：

只适用于能将自变量与因变量分离的情形，这一类微分方程可以写成如下形式

(8)
\[q(y)dy=p(t)dt. \]

由于方程的左边是仅含$y$的函数，而右边是仅含$t$的函数，因此通过对两边直接积分就得到解。 这种方法称为变量分离法，是解析求解微分方程的一种最基本的方法。

在实际中很少能直接得到微分方程的公式解，即给出$y(t)$关于$t$的函数的显式表达。

代之，通常构造下列形式的函数值表：

\[\begin{tabular}{c|c|c|c|c|c} $t_0$ & $t_1$ & $t_2$ & $t_3$ & $\cdots$ & $t_n$\\ \hline $y_0$ & $y_1$ & $y_2$ & $y_3$ & $\cdots$ & $y_n$ \end{tabular} \]

即所谓的数值解。

Taylor级数： 举个例子

(9)
\[\left\{\begin{array}{l} y'=t^2 + \cos t - \sin y \\ y(-1)=3 \end{array}\right. \label{ode-example-1} \]

考虑$y$的Taylor展开级数，我们记之为

\[y(t+h)=y(t)+h y'(t)+\frac{h^2}{2!}y''(t)+\frac{h^3}{3!}y'''(t)+\cdots \]

上式中的导数可从微分方程中依据链式法则得到，它们是

\[\begin{array}{lcl} y'' &=& 2t - \sin t - y'\cos y \\ y''' &=& 2 - \cos t - y''\cos y + (y')^2\sin y \\ &\vdots& \end{array} \]

下面是求解初值问题的一个算法，它从$t_0=-1$开始，步长取$h=0.01$. 我们求$y(t)$在区间$[-1,1]$上的数值解，要执行 $N=200$ 步。

  • 输入 $t_0,y_0,h$
  • for $k=1$ to $N$ do
    1. $y'\leftarrow t^2+\cos t-\sin y$
    2. $y''\leftarrow 2t-\sin t-y'\cos y$
    3. $y'''\leftarrow 2-\cos t-y''\cos y+(y')^2\sin y$
    4. $y\leftarrow y+h(y'+\frac{h}{2}(y''+\frac{h}{3}(y''')))$
    5. $t\leftarrow t+h$
    6. 输出 $t(k),y(k)$
  • end do

Euler方法

上述Taylor级数法中，特殊地只取到一阶近似即称为Euler方法，它具有如下形式

(10)
\[y(t+h)\approx y(t)+h y'(t)=y(t)+h f(t,y). \]

mm3-ode-euler

Runge-Kutta法

从$y(t+h)$的Taylor级数入手

\[y(t+h)=y(t)+h y'(t)+\frac{h^2}{2!}y''(t)+\frac{h^3}{3!}y'''(t)+\cdots \]

根据微分方程我们有

(11)
\[\begin{array}{lcl} y' &=& f \\ y'' &=& f_t + f_y y' = f_t + f_y f \\ y''' &=& f_{tt} + f_{ty} f + (f_t+f_y f)f_y + f(f_{yt}+f_{yy}f) \\ &\vdots& \end{array} \]

取展开式中前三项可写成如下形式

(12)
\[\begin{array}{ll} y(t+h) &= y + h f + \frac{h^2}{2}(f_t + f_y f) + O(h^3)\\ &= y + \frac12 h f + \frac12 h(f + h f_t + h f f_y) + O(h^3) \end{array} \]

又由于

(13)
\[f(t+h,y+hf) = f + h f_t + h f f_y + O(h^2), \]

$y(t+h)$的二阶展开式(11)可改写成

(14)
\[y(t+h) = y + \frac12 h f + \frac12 h f(t+h,y+hf) + O(h^3) \]

因此，步进求解的迭代公式是

(15)
\[y(t+h) = y(t) + \frac12 h f(t,y) + \frac12 h f(t+h,y+hf(t,y)) \]

等价地可写成

(16)
\[\left\{\begin{array}{l} F_1 = h f(t,y) \\ F_2 = h f(t+h,y+F_1) \\ y(t+h) = y(t) + \frac12(F_1+F_2) \end{array}\right. \label{runge-kutta-2nd-equal} \]

反复地使用这个公式，每次就前进一步求解。 它被称为二阶 Runge-Kutta法。

下面是经典的四阶Runge-Kutta方法：

(17)
\[y(t+h) = y(t) + \frac16(F_1+2F_2+2F_3+F_4), \label{runge-kutta-4th-5ev} \]

其中

(18)
\[\left\{\begin{array}{l} F_1 = h f(t,y), \\ F_2 = h f(t+\frac12 h,y+\frac12 F_1), \\ F_3 = h f(t+\frac12 h,y+\frac12 F_2), \\ F_4 = h f(t+h,y+F_3). \end{array}\right. \nonumber \]

Python库

首先是用sympy的dsolve解常微分方程，直接贴代码 

import sympy as sy
 
def differential_equation(x,f):
    return sy.diff(f(x),x,2)+f(x)#f(x)''+f(x)=0 二阶常系数齐次微分方程
x=sy.symbols('x')#约定变量
f=sy.Function('f')#约定函数
print(sy.dsolve(differential_equation(x,f),f(x)))#打印
sy.pprint(sy.dsolve(differential_equation(x,f),f(x)))#漂亮的打印

数值： scipy.integrate.odeint

这个用起来就比较麻烦了，不过用于画图还是很棒的。先看一个一阶微分方程的例子。 

import numpy as np
from scipy.integrate import odeint
#一阶微分方程的例子
def diff_equation(y,x):
    #dy/dx=y，其实手工解的话，很容易知道，y=Cexp(x)
    return np.array(y)#微分方程格式，左边一定是dy/dx,返回右边
x=np.linspace(0,1,num=100)#初始点是0
result=odeint(diff_equation,1,x)#中间那个是y0初值，即x=0时y=1
plt.plot(x,result[:,0])#result整个矩阵的第一列
plt.grid()#网格
plt.show()#这是y=exp(x)的图像

二阶的话，就有点麻烦了。 

import numpy as np
from scipy.integrate import odeint
#二阶微分方程
def diff_equation(y_list,x):
    #y''+y=0 二阶的话，要换成两个一阶微分方程的方程组
    #设y'=z
    #那么z'=y''=-y
    #手工解也很容易知道是 y=C1sin(x)+C2cos(x)
    y,z=y_list
    return np.array([z,-y])
x=np.linspace(0,np.pi*2,num=100)
y0=[1,1]#y(0)=1,y'(0)=1
result=odeint(diff_equation,y0,x)
plt.plot(x,result[:,0],label='y')#y的图像，y=cos(x)+sin(x)
plt.plot(x,result[:,1],label='z')#z的图像，也就是y'的图像，z=-sin(x)+cos(x)
plt.legend()
plt.grid()
plt.show()

Mathematica

  • 在Mathematica中，可以用 NDSolve 来数值求解微分方程
  • 在Mathematica中，可以用 Solve 来求微分方程的解析解

万有引力定理的推导

万有引力定理的推导。

背景：

  1. 哥白尼(1473-1543)日心说（1543）
  2. 伽利略(1564-1642)惯性原理，落体定理
  3. 开普勒(1571-1630)行星运动的三大定理

行星运动的受力情况是什么？胡克（英国），惠更斯（荷兰）都取得了一些成果。牛顿(1642-1703)以微积分为工具，在开普勒三定律和牛顿第二定律的基础上，用演绎法得到了万有引力定律（1687,《自然科学之数学原理》）

模型假设: 已知行星的运行轨道是椭圆，太阳在一个焦点上；行星在每个时间内行星-太阳向径扫过的面积是常数；行星运行周期与平均半径的1.5次方成正比；受力为质量乘加速度。

  1. 以太阳为极点，建立极坐标$(r,\theta)$，则椭圆方程
    \[r=\frac{p}{1+e\cos\theta}, p=\frac{b^2}a, b^2=a^2(1-e^2) \]
    $a$, $b$为长、短半轴，$e$为离心率。用$\vec r(t)=(r(t),\theta(t))$表示行星的位置。
  1. 单位时间内向径$r$扫过的面积是常数$A$，即
    \[\frac12r^2\frac{d\theta}{dt}=A \]
  2. 运行周期$T$满足
    \[T^2=\lambda a^3 \]
  3. 行星运行时受力
    \[\vec f=m\frac{d^2\vec r}{dt^2} \]

推导： 先引入2个基向量

\[\begin{cases} \vec u_r=(\cos\theta, \sin\theta) \\ \vec u_\theta=(-\sin\theta,\cos\theta) \end{cases} \]

即，$\vec u_r$指向向径方向，$u_\theta$与$u_r$垂直。这样

\[\begin{cases} \frac{d\vec u_r}{dt}=(-\sin\theta\frac{d\theta}{dt}, \cos\theta\frac{d\theta}{dt})=\vec u_\theta\frac{d\theta}{dt} \\ \vec u_\theta=(-\cos\theta,-\sin\theta)\frac{d\theta}{dt}=-\frac{d\theta}{dt}\vec u_r \end{cases} \]

由

\[\vec r=r\vec u_r \]

因此

\[\frac{d\vec r}{dt}=\frac{dr}{dt}\vec u_r+r\frac{d\theta}{dt}\vec u_\theta \]
\[\frac{d^2\vec r}{dt^2}=(\frac{d^2r}{dt^2}-r\frac{d\theta}{dt}^2)\vec u_r +(r\frac{d^2\theta}{dt^2}+2\frac{dr}{dt}\frac{d\theta}{dt})\vec u_\theta \]

由

\[\frac12r^2\frac{d\theta}{dt}=A \]

有

\[\frac{d\theta}{dt}=\frac{2A}{r^2}, \frac{d^2\theta}{dt^2}=\frac{-4A}{r^3}\frac{dr}{dt} \]

可以看到$r\frac{d^2\theta}{dt^2}+2\frac{dr}{dt}\frac{d\theta}{dt}=0$，因此

\[\frac{d^2\vec r}{dt^2}=(\frac{d^2r}{dt^2}-r\frac{d\theta}{dt}^2)\vec u_r \]

由椭圆方程，可得

\[\frac{dr}{dt}=p\frac{-1}{(1+e\cos\theta)^2}(-e\sin\theta)\frac{d\theta}{dt} =\frac1p r^2e\sin\theta\frac{d\theta}{dt}=\frac{2Ae}{p}\sin\theta \]

和

\[\frac{d^2r}{dt}=\frac{2Ae}{p}\cos\theta\frac{d\theta}{dt} =\frac{4A^2e}{pr^2}\cos\theta=\frac{4A^2(p-r)}{p r^3} \]

则有，

\[\frac{d^2\vec r}{dt^2}=\frac{-4A^2}{p r^2}\vec u_r \]

所以，行星运动时，受到指向太阳的力

\[\vec f=\frac{-4A^2m}{p r^2}\vec u_r \]

最后，还需要说明$\frac{A^2}{p}$与行星无关。

由行星运行一周，正好扫过整个椭圆，则

\[TA=\pi a b \]

可得

\[\frac{A^2}{a^2b^2}=\frac{\pi^2}{T^2} \]

又

\[\frac{A^2}{a^2b^2}=\frac{A^2}{a^2 pa}=\frac{A^2}{p\frac{T^2}{\lambda}} \]

因此

\[\frac{A^2}{p}=\frac{A^2}{a^2b^2}\frac{T^2}{\lambda}=\frac{\pi^2}{\lambda} \]

因此，太阳对行星的力为

\[\vec f=-\frac{4\pi^2}{\lambda}\frac{m}{r^2}\vec u_r \]

其中$\frac{4\pi^2}{\lambda}$与行星无关，是跟太阳相关的量

人口模型

中国是一个人口大国,人口问题始终是制约我国发展的关键因素之一。 根据已有数据,运用数学建模的方法,对中国人口做出分析和预测是一个重 要问题。 近年来中国的人口发展出现了一些新的特点,例如,老龄化进程加速、 出生人口性别比持续升高,以及乡村人口城镇化等因素,这些都影响着中国 人口的增长。

试从中国的实际情况和人口增长的上述特点出发,参考相关数 据(也可以搜索相关文献和补充新的数据),建立中国人口增长的数学模型 ,并由此对中国人口增长的中短期和长期趋势做出预测;特别要指出你们模 型中的优点与不足之处。

假设

  1. 单位时间人口增长与当前人口数成比例，这个比例（人口增长率）为常数

模型一(Malthus, 1798)

设$x(t)$是时间$t$时刻的人口总数，人口增长率为$r$，则有微分方程

\[\begin{cases} \frac{dx(t)}{dt}=r x(t) \\ x(0)=x_0 \end{cases} \]

可以得到

\[x(t)=x_0e^t \]

为指数增长模型

改进的指数增长模型 人口增长率为时间的函数$r(t)$，则

\[\begin{cases} \frac{dx(t)}{dt}=r(t) x(t) \\ x(0)=x_0 \end{cases} \]

利用实测的数据，可以近似得到$r(t)$的表达式。

适用情形

  1. 19世纪前欧洲一些国家
  2. 加拿大初期欧洲移民后代
  3. 满清入关后的中国（所谓大乱后必有大治）

假设

人口增长到一定程度后，增长率下降的主要原因是，自然资源、环境条件等因素对人口的增长起阻滞作用。 而且随着人口的增长，阻滞作用也越大。

模型二(logistic模型) 人口的增长率$r$也是人口数$x$的函数，$r=r(x)$。简单地取$r=a+bx$。

  1. 理论上$x=0$时，增长率最大为$r$。即$r(0)=r$，则$a=r$
  2. $x_m$是环境所能容纳的最大人口，则当$x=x_m$时，增长率为$0$，即$r(x_m)=0$

可以得到$r(x)=r(1-\frac{x}{x_m})$，人口模型方程

\[\begin{cases} \frac{dx(t)}{dt}=r\left(1-\frac{x(t)}{x_m}\right) x(t) \\ x(0)=x_0 \end{cases} \]

分离变量后，得到解

\[x(t)=\frac{x_m}{1+\left(\frac{x_m}{x_0}-1\right)e^{-rt}} \]

计算参数$r$, $x_0$, $x_m$。

  1. 由方程

    \[\frac{\frac{dx}{dt}}{x}=r-\frac{r}{x_m}x \]

    若有历年的人口数据$\left\{(t_n,x_n)\right\}$，则取$x_1$为第一个数据，用数据$q_n=\frac{x_{n+1}-x_n}{x_n}$与$x_n$做线性拟合$q=a+bx$，就可以得到$r$与$x_m$。

  2. 直接用数据做非线性最小二乘法

血液中的酒精浓度

  1. 据报载,2003年全国道路交通事故死亡人数为10.4372万,其中因饮酒驾车造成的占有相当的比 例。针对这种严重的道路交通情况,国家质量监督检验检疫局2004年5月31日发布了新的《车辆驾驶 人员血液、呼气酒精含量阈值与检验》国家标准,新标准规定,车辆驾驶人员血液中的酒精含量大 于或等于20毫克/百毫升,小于80毫克/百毫升为饮酒驾车,血液中的酒精含量大于或等于80毫克/百毫升为醉酒驾车。
  2. 大李在中午12点喝了一瓶啤酒,下午6点检查时符合新的驾车标准,紧接着他在吃晚饭时又喝 了一瓶啤酒,为了保险起见他呆到凌晨2点才驾车回家,又一次遭遇检查时却被定为饮酒驾车,这让 他既懊恼又困惑,为什么喝同样多的酒,两次检查结果会不一样呢?

请你参考相关的数据(或自己收集资料)建立饮酒后血液中酒精含量的数学模型

假设

  1. 酒到胃液后，被全部吸收
  2. 近似地将胃以外的人体看作是会吸收酒精的容器，其中的酒精浓度就是检测出来的血液中的酒精浓度
  3. 吸收的酒精会被分解排出（如，肝脏可以分解酒精）
  4. 整个过程，胃和人体的血液体积保持不变
  5. 胃中的酒精吸收入人体的速度与胃中酒精的含量成比例
  6. 人体中的酒精排出酒精的速度与人体中的酒精的含量成比例

模型

  1. 设$x_1(t)$是胃中的酒精含量，$x_2(t)$是身体血液中的酒精含量。它们的精度浓度分别为$c_1(t)$和$c_2(t)$
  2. $V_1$, $V_2$为胃和身体的体积，为常数
  3. 设酒精的吸收和排出速率为$k_1$, $k_2$
  4. 摄入的酒精总量为$D$

可以得到模型

\[\begin{cases} \frac{d x_1(t)}{dt}=-k_1 x_1(t) \\ \frac{d x_2(t)}{dt}=-k_2 x_2(t)+k_1 x_1(t) \\ x_1(0)=D \\ x_2(0)=0 \end{cases} \]

注意到$c_1(t)V_1=x_1(t)$, $c_2(t)V_2=x_2(t)$，则有

\[\begin{cases} \frac{d c_1(t)}{dt}=-k_1 c_1(t) \\ \frac{d c_2(t)}{dt}=-k_2 c_2(t)+k_1 c_1(t)\frac{V_1}{V_2} \\ x_1(0)=D \\ x_2(0)=0 \end{cases} \]

可以解得

\[\begin{cases} c_1(t)=D e^{-k_1 t} \\ c_2(t)=D\frac{k_1 V_1}{V_2(k_2-k_1)}(e^{-k_1t}-e^{-k_2t}) \end{cases} \]

式子中，$V_1$, $V_2$, $k_1$, $k_2$都是因人而异的。

当成人摄入2瓶啤酒时，可以设置$D=2$。同时测量经过多个时间点后的身体血液中的精度浓度的数据，来估计$\frac{k_1 V_1}{V_2(k_2-k_1)}$, $k_1$, $k_2$的值。

数据

体重约70kg的某人在短时间内喝下2瓶啤酒后,隔一定时间测量他的血液中酒精含量(毫克/百毫升),得到数据如下:

时间(小时) 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5
酒精含量 30 68 75 82 82 77 68 68 58 51 50 41
时间(小时) 6 7 8 9 10 11 12 13 14 15 16
酒精含量 38 35 28 25 18 15 12 10 7 7 4

利用Matlab(函数lsqcurvefit)可以得到

\[c_2(t)=61.955(e^{-0.089t}-e^{-6.603t}) \]

大李中午喝酒后，到下午6点的酒精浓度为

\[c_2(6)=\frac12*61.955*(e^{-0.089*6}-e^{-6.603*6})=18.16<20 \]

不算酒驾

假定他再次喝酒是晚上7点，则凌晨2点时的酒精浓度为（需要加上中午的酒）

\[c_2(14)+c_2(7)=8.910+16.614>20 \]

是饮酒驾驶

捕食者-食饵模型

自然界中，不同种群之间存在着一种既有依存、又有制约的生存方式：种群甲靠自然资源生长，种群乙靠捕食甲为生。如：食用鱼和鲨鱼，落叶松和蚜虫。称种群甲为食饵(prey)，种群乙为捕食者(predator)。

问：农田中有一种害虫吃农作物，有一种益虫以害虫为食。现在有一种农药，会同时杀死这两种虫。是否应用使用这种农药？

模型

设食饵（食用鱼）和捕食者（鲨鱼）在时刻$t$的数量为$x(t)$, $y(t)$。假定，在单独生存时，由于资源丰富，食饵以指数规律成长，即$\frac{dx}{dt}=rx$。而捕食者的存在，使食饵的增长率减少，设减少率与捕食者数量成正比，则有

\[\frac{dx}{dt}=x(r-ay)=rx-ayx \]

比例系数$a$反映捕食者的掠食能力。

捕食者无法离开食饵生存，设它独自存在时死亡率为$d$，即$\frac{dy}{dt}=-dy$。而食饵的存在，减小了捕食者的死亡率，促使其增长，且增长率与食饵的数量成正比。则有

\[\frac{dy}{dt}=y(-d+bx)=-dy+bxy \]

其中系数$b$反映食饵对捕食都的供养能力

取$r=1$, $d=0.5$, $a=0.1$, $b=0.02$, $x_0=25$, $y_0=2$，利用数值方法，可以得到

mm3-ex-volterra

呈周期性变化

由

\[\frac{dy}{dt}=y(-d+bx)=-dy+bxy \]

则

\[x=\frac1b\left({\frac{dy/dt}{y}+d}\right) \]

在一个周期$T$内积分，可以得到食饵的平均值，

\[\begin{aligned} \frac1T\int_0^Tx(t)dt=&\frac1T\int_0^T\frac1b\left[{\frac{dy/dt}{y}+d}\right]dt \\ =&\frac{d}{b}T\frac1T+\frac1T\int_0^T\frac1b\frac{dy}{dt}\frac1ydt \\ =&\frac{d}{b}+\frac1{Tb}(\ln y(T)-\ln y(0))=\frac{d}{b} \end{aligned} \]

类似地，捕食者的平均值

\[\bar y=\frac{r}{a} \]

若由于捕捞，导致食饵的增长率由$r$下降到$r-e$，而捕食者的死亡率由$d$上升到$d+e$。此时，有

\[\bar y=\frac{r-e}{a} , \bar x=\frac{d+e}{b} \]

若控制捕捞能力$e$（如由于战争导致捕捞率下降），则$\bar y$会上升，$\bar x$会下降。也就说，捕捞能力下降，会利于捕食者的种群增长。

类似地，对于农作物-害虫-益虫系统，若有某种药物在杀死害虫的同时，也会杀死益虫，则会导致害虫增加

类似地，对于农作物-害虫-益虫系统，若有某种药物在杀死害虫的同时，也会杀死益虫，则会导致害虫增加

网络消息的传播

网络消息的传播模型

与传染病模型是类似的

模型一：假设，某地区的总人口为$N$，在短期内不变。$x(t)$表示知道消息的人数所占的比例，初始时刻比例为$x_0<1$，传播率为$h$，则可以建立数学模型

\[\begin{cases} \dfrac{dx}{dt} N=hN(1-x) \\ x(0)=x_0 \end{cases} \]

可以得到解为$x(t)=(x_0-1)e^{-ht}+1$，且$\displaystyle\lim_{t\to+\infty}x(t)=1$，与实际情况不合。

模型二： 实际情况下，传播率为$h$，但会有一部分人得到消息后，并不去传播。设不传播率为$r$，则有

\[\begin{cases} \frac{dx}{dt}N=N[h-(h+r)x] \\ x(0)=x_0 \end{cases} \]

可以得到解为

\[x(t)=(x_0-\frac{h}{h+r})e^{-(h+r)t}+\frac{h}{h+r} \]

有

\[\lim_{t\to+\infty}x(t)=\frac{h}{h+r}<1 \]

比较符合实际

火箭发射

为什么发射火箭通常用三级发射

  1. 卫星进入600km高空轨道，需要的最低速度是多少？

假定:

  • 卫星以地球中心为圆心的某个平面的圆周上，并绕地球做匀速加速运动
  • 地球是固定的一外均匀球体，质量集中于球心
  • 其它卫星的引力不计

模型： 卫星的身心力应该与地球的引力相等

\[\frac{GMm}{r^2}=\frac{mv^2}r \]

其中$M$为地球质量，$m$为卫星质量，$r$为卫星到球心的距离，$v$为卫星的速度，$G$为万有引力常数。 注意到，在地面上

\[mg=\frac{GMm}{R^2} \]

其中$R$是地球的半径，$g$为地面的重力加速度。联立后，有

\[v=R\sqrt{\frac{g}{r}} \]

取$R=6400km$, $g=9.81m/s^2$，$r-R=600km$，可得

\[v\approx 7.6 km/s \]
  1. 火箭推力及升空速度

假定:

  • 火箭在喷气推动下做直线运动，不计火箭的重力和空气阻力
  • $t$时刻火箭质量为$m(t)$，速度$v(t)$，且为时间的连续可微函数
  • 火箭末端喷出气体的速度相对火箭自身为常数$u$
  • 需要把卫星的最后速度为$7.6km/s$

模型: 假定$dt$时间内，质量变化为$dm$，速度变化为$dv$，则喷出气体的总动量为$dm(v-u)$，而火箭的动量变化为$(m+dm)(v+dv)-mv$，由动量守恒，有

\[(m+dm)(v+dv)-mv=- dm(v-u) \]

即

\[(m+dm)dv=-u dm \]

令$dt\to0$，则有微分方程

\[m\frac{dv}{dt}=-u\frac{dm}{dt} \]

利用分离变量，加上初值条件

\[v(0)=v_0 , m(0)=m_0 \]

可得

\[v(t)=v_0+u\ln\frac{m_0}{m(t)} \]

因此，最后的升空速度$v(t)$与火箭推力$u$和最终的质量比$\frac{m_0}{m(t)}$有关。

对于一级运载火箭，$m_p$为有效负载（如卫星质量），$m_F$为燃料质量，$m_s$为结构质量（如外壳、推进器）。

  • 目前技术条件，相对火箭的喷气速度$u=3km/s$，及
    \[\frac{m_s}{m_s+m_F}\geq \frac19 \]
  • 初速度$v_0=0$

则有最终速度（其中$m_0=m_p+m_F+m_s$）

\[v=u\ln\frac{m_0}{m_p+m_s} \]

令$m_s=\lambda(m_F+m_s)=\lambda(m_0-m_p)$，代入上式

\[v=u\ln\frac{m_0}{\lambda m_0+(1-\lambda)m_p} \]

即使卫星质量为$0$，火箭的速度上限为

\[v=u\ln\frac1\lambda \]

取$u=3km/s$, $\lambda=1/9$，有$v=3\ln 9\approx 6.6 km/s$。可以看到，目前的技术条件下，一级火箭不足以送卫星上需要的轨道。

从前面的结果中可以看出，最终速度达不到要求，与火箭始终在推动整个结构有关。如果在整个过程中，可以把不断地把无用的结构抛弃的话，效率就会大大提高了。

理想火箭 理想状态下，火箭不断地把不需要的结构抛弃

模型 若时间$dt$内，质量变化$dm$，其中燃料占$1-\alpha$，其余的不用的结构占$\alpha$，则由动量守恒

\[(m+dm)(v+dv)-mv=- (1-\alpha)dm(v-u)+\alpha \cdot dm \cdot v \]

即有微分方程

\[-m(t)\frac{dv(t)}{dt}=(1-\alpha)\frac{dm(t)}{dt} u \]

及初始条件$v(0)=0$, $m(0)=m_0$。

可得

\[v(t)=(1-\alpha)u\ln{m_0}{m(t)} \]

最终，只剩下卫星后的速度

\[v=(1-\alpha)u\ln{m_0}{m_p} \]

因此，只要$m_0$足够大，总能使卫星达到需要的任意速度。

考虑到空气阻力和重力的影响，估计需要$v=10.5km/s$，取$u=3km/s$, $\alpha=0.1$，则可以得到$\frac{m_0}{m_p}=50$，即发射1吨重的卫星，需要50吨的理想火箭。

多级火箭 多级火箭自末级开始，逐级燃烧，当第$i$级燃料烧尽后，第$i+1$级自动点火，并抛弃第$i$级。用$m_i$表示第$i$级火箭质量，$m_p$表示卫星质量。假设：

  • 各级火箭具有相同的$\lambda$，$(1-\lambda)m_i$表示第$i$级的燃料质量。
  • 各级火箭的喷气相对于火箭的速度$u$相同
  • 各级火箭的初始质量与其负载质量之比为常数$k$，
    \[k=\frac{m_i}{m_{i+1}+m_{i+2}+\cdots+m_n+m_p} , i=1,2,\cdots,n \]

对于二级火箭，当第一级火箭燃烧完后，其速度为

\[v_1=u\ln\frac{m_1+m_2+m_p}{\lambda m_1+m_2+m_p}=u\ln\frac{k+1}{\lambda k+1} \]

第二级火箭燃烧完后，速度为

\[v_2=v_1+u\ln\frac{m_2+m_p}{\lambda m_2+m_p}=2u\ln{k+1}{\lambda k+1} \]

取$v_2=10.5km/s$，$\lambda=0.1$，$u=3km/s$有

\[6\ln\frac{k+1}{0.1k+1}=10.5 \]

可以得$k=11.2$，从而

\[\frac{m_0}{m_p}=\frac{m_1+m_2+m_p}{m_p}=(k+1)^2\approx 149 \]

同样，对于三级火箭，$\frac{m_0}{m_p}\approx 77$。

级数 1 2 3 4 5 .. 无穷
$m_0/m_p$ 无 149 77 65 60 … 50

综合考虑经济效益，使用三级火箭是最好的方案。

目录

本节读完

例 1. 谢

1.

例 2. 传染病模型