数学建模方法之微分方程

数学建模方法之微分方程

图与网络还没有完全分析完，后续会继续有关内容的补充。因为组会要求，所以先提前完成这篇有关微分方程的blog！！！冲！！！

---

肝！

---

前言

微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步:
（1）根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。
（2）找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。
（3）运用这些规律列出方程和定解条件。

---

列方程常见的方法有以下几种：
（1）按规律直接列方程
在数学、力学、物理、化学等学科中许多自然现象所满足的规律已为人们所熟悉，并直接由微分方程所描述。如牛顿第二定律、放射性物质的放射性规律等。我们常利用这些规律对某些实际问题列出微分方程。

（2）微元分析法与任意区域上取积分的方法
自然界中也有许多现象所满足的规律是通过变量的微元之间的关系式来表达的。对于这类问题，我们不能直接列出自变量和未知函数及其变化率之间的关系式，而是通过微元分析法，利用已知的规律建立一些变量（自变量与未知函数）的微元之间的关系式，然后再通过取极限的方法得到微分方程，或等价地通过任意区域上取积分的方法来建立微分方程。

（3）模拟近似法
在生物、经济等学科中，许多现象所满足的规律并不很清楚而且相当复杂，因而需要根据实际资料或大量的实验数据，提出各种假设。在一定的假设下，给出实际现象所满足的规律，然后利用适当的数学方法列出微分方程。

---

在实际的微分方程建模过程中，也往往是上述方法的综合应用。不论应用哪种方法，通常要根据实际情况，做出一定的假设与简化，并要把模型的理论或计算结果与实际情况进行对照验证，以修改模型使之更准确地描述实际问题并进而达到预测预报的目的。

---

一、常微分方程问题的数学模型

在数学、力学、物理、化学等学科中已有许多经过实践检验的规律和定律，如牛顿运动定律、基尔霍夫电流和电压定律、物质的放射性规律、曲线的切线的性质等，这些都涉及某些函数的变化率。我们就可以根据相应的规律，列出常微分方程。

【例1】建立物体冷却过程的数学模型
将某物体放置于空气中，在时刻$t=0$时，测量得它的温度为$u_0=100℃$，20分钟后测量得它的温度为$u_1=60℃$。要求建立此物体的温度$u$和时间$t$的关系，并计算经过多长此物体的温度将达到30℃，其中我们假设空气的温度保持为20℃。
【解】
Newton冷却定律是温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律：当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比，比例系数称为热传递系数，记为$k$。

假设该物体在时刻t时的温度为$u=u(t)$，则由Newton冷却定律，得到
$$\frac{du}{dt}=-k(u-20),(6.1)$$

其中，$k>0$，方程(6.1)就是物体冷却过程的数学模型。

可将方程(6.1)改写为
$$\frac{d(u-20)}{u-20}=-kdt,$$
两边积分得到
$${\int }_{100}^u\frac{d(u-20)}{u-20}={\int }_0^t-kdt,$$
化简得
$$u=20+80e^{-kt}.(6.2)$$

把条件$t=20$时，$u=u_1=60$代入式(6.2)，得$k=\frac{ln2}{20}$，所以此物体的温度u和时间t的关系为$u=20+80e^{-\frac{ln2}{20}t}$。令$30=20+80e^{-\frac{ln2}{20}t}$ ，解之得$t=60$，即60分钟后物体的温度为30℃。

【例2】 目标跟踪问题
设位于坐标原点的甲舰向位于$x$轴上点$Q_0(1,0)$处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度$v_0$（$v_0$是常数）沿平行于$y$轴的直线行驶，导弹的速度是$5v_0$，求导弹运行的曲线。又乙舰行驶多远时，导弹将它击中？
【解】
设导弹的轨迹曲线为$y=y(x)$，并设经过时间 ，导弹位于点$P(x,y)$,乙舰位于点$Q(1,v_{0}t)$。由于导弹头始终对准乙舰，故此时直线$PQ$就是导弹的轨迹曲线弧$OP$在点$P$处的切线，如图6.1。则有
$$\frac{dy}{dx}=\frac{v_{0}t-y}{1-x}.(6.2)$$

由于模型中含有参变量$t$，故要求解该模型应增加附件条件。解决这个问题可从问题描述中寻求办法。

方法一：任意两个匀速运动的物体，相同时间内所经过的距离与其速度成正比。
由已知，导弹的速度5倍于乙舰，即在同一时间段$t$内，导弹运行轨迹的总长亦应5倍于乙舰，即 $\widehat{OP}$的弧长5倍于线段$\overline{Q_{0}Q}$的长度。由弧长计算公式可得
$${\int }_0^x\sqrt{1+{\left(\frac{dy}{dx}\right)}^2}dx=5v_{0}t，(6.4)$$
方程两边关于$x$求导，得
$$\sqrt{1+{\left(\frac{dy}{dx}\right)}^2}=5v_0\frac{dt}{dx}.(6.5)$$

方法二：利用速度分量合成的概念。
由于在点$P(x,y)$导弹的速度
恒为$5v_0$，而该点的速度大小等于该点在$x$轴和$y$轴上的速度分量$x^{\prime }(t),y^{\prime }(t)$的合成，故有
$$\sqrt{{\left(\frac{dx}{dt}\right)}^2+{\left(\frac{dy}{dt}\right)}^2}=5v_0，(6.6)$$
或改写成
$$\frac{dx}{dt}\sqrt{1+{\left(\frac{dy}{dx}\right)}^2}=5v_0.(6.7)$$
两边同除以$\frac{dx}{dt}$即得(6.5)。由此可见，利用弧长的概念或速度的概念得到的结果是一致的。

为了消去中间变量$t$，把方程(6.3)改写为
$$(1-x)\frac{dy}{dx}=v_{0}t-y.(6.8)$$
然后两边关于$x$求导，得
$$(1-x)\frac{d^{2}y}{d{x}^2}-\frac{dy}{dx}=v_0\frac{dt}{dx}-\frac{dy}{dx}，$$
整理后，得
$$(1-x)\frac{d^{2}y}{d{x}^2}=v_0\frac{dt}{dx}，(6.9)$$

联立(6.5)和(6.9)得
$$\{\begin{array}{l}\sqrt{1+(\frac{dy}{d2}{)}^2}=5v_0\frac{dt}{dx}\\ (1-x)\frac{d^{2}y}{d{x}^2}=v_0\frac{dt}{dx}\end{array}$$

消去中间变量$\frac{dt}{dx}$，得关于轨迹曲线的二阶非线性常微分方程：
$$(1-x)\frac{d^{2}y}{d{x}^2}=\frac{1}{5}\sqrt{1+(\frac{dy}{d2}{)}^2},0<x<1.$$

要求此问题的定解，还需要给出两个初始条件。事实上，初始时刻轨迹曲线通过坐标原点，即$x=0$时，$y(0)=0$；此外在该点的切线平行于$x$轴，因此有$y^{\prime }(0)=0$。归纳可得导弹轨迹问题的数学模型为
$$\{\begin{array}{l}(1-x)\frac{d^{2}y}{d{x}^2}=\frac{1}{5}\sqrt{1+(\frac{dy}{d2}{)}^2},0<x<1,\\ y(0)=0,y^{\prime }(0)=0.\end{array}(6.10)$$
此模型为二阶常微分方程初值问题。求解此类问题，通常采用降价法。

令$p=y^{\prime }，$则KaTeX parse error: Undefined control sequence: \pprime at position 3: y^\̲p̲p̲r̲i̲m̲e̲=\frac{dp}{dx}，则(6.10)式变为关于p的常微分方程初值问题：
$$\{\begin{array}{l}(1-x)\frac{dp}{dx}=\frac{1}{5}\sqrt{1+p^2},0<x<1,\\ p(0)=0.\end{array}(6.11)$$
利用分离变量法，求解并代入初始条件得
$$ln\left(p+\sqrt{1+p^2}\right)=-\frac{1}{5}ln(1-x)，$$
化简得
$$p+\sqrt{1+p^2}=(1-x{)}^{-1/5}.$$
为求得p的显式表达式，利用上式作如下等式变换：
$$-p+\sqrt{1+p^2}=\frac{1}{p+\sqrt{1+p^2}}=(1-x{)}^{1/5}.$$
以上两式相减，得关于$p$的表达式，从而得到关于$y$的一阶微分方程的初值问题：
$$\{\begin{array}{l}\frac{dy}{dx}=p=\frac{1}{2}[(1-x{)}^{-\frac{1}{5}}-(1-x{)}^{\frac{1}{5}}],\\ p(0)=0.\end{array}$$
求解此微分方程，即得导弹运行的轨迹曲线方程为：
$$y=-\frac{5}{8}(1-x{)}^{4/5}+\frac{5}{12}(1-x{)}^{6/5}+\frac{5}{24}.(6.12)$$
如何击中乙舰？在解（6.12）中，令$x=1$，得$y=\frac{5}{24}$，即在（1，$\frac{5}{24}$）处击中乙舰。
何时击中乙舰“？击中乙舰时，乙舰航行的距离$y=\frac{5}{24}$，由$y=vt$，得$t=\frac{5}{24v_0}$时击中乙舰。

---

二、传染病预测问题

· 问题提出
世界上存在着各种各样的疾病，许多疾病是传染的，如SARS、艾滋病、禽流感等，每种病的发病机理与传播途径都各有特点。如何根据其传播机理预测疾病的传染范围及染病人数等，对传染病的控制意义十分重大。

1. 指数传播模型

基本假设
（1）所研究的区域是一封闭区域，在一个时期内人口总量相对稳定，不考虑人口的迁移（迁入或迁出）。
（2）$t$时刻染病人数$N(t)$是随时间连续变化的、可微的函数。
（3）每个病人在单位时间内的有效接触（足以使人致病）或传染的人数为$\lambda$（$\lambda >0$为常数）。

模型建立与求解
记N(t)为t时刻染病人数，则 时刻的染病人数为$N(t+\Delta t)$。从$t\to t+\Delta t$时间内，净增加的染病人数为$N(t+\Delta t)-N(t)$，根据假设（3），有
$$N(t+\Delta t)-N(t)=\lambda N(t)\Delta t.$$
若记$t=0$时刻，染病人数为$N_0$，则由假设（2），在上式两端同时除以$\Delta t$，并令$\Delta t\to 0$，得传染病染病人数的微分方程预测模型：
$$\{\begin{array}{ll}\frac{dN(t)}{dt}=\lambda N(t),& t>0\\ N(0)=N_0.& \end{array}(6,13)$$

利用分离变量法可很容易地得到该模型的解析解为
$$N(t)=N_0{e}^{\lambda t}.$$

结果分析与评价
模型结果显示传染病的传播是按指数函数增加的。一般而言在传染病发病初期，对传染源和传播路径未知，以及没有任何预防控制措施的情况下，这一结果是正确的。此外，我们注意到，当$t\to \infty$时，$N(t)\to \infty$，这显然不符合实际情况。

事实上，在封闭系统的假设下，区域内人群总量是有限的。预测结果出现明显失误。为了与实际情况吻合，有必要在原有基础上修改模型假设，以进一步完善模型。

2. SI模型

基本假设
（1）在传播期内，所考察地区的人口总数为$N$，短期内保持不变，既不考虑生死，也不考虑迁移。
（2）人群分为易感染者（susceptible）和已感染者（infective），即健康人群和病人两类。
（3）设$t$时刻两类人群在总人口中所占的比例分别为$s(t)$和$i(t)$，则$s(t)+i(t)=1$。
（4）每个病人在单位时间（每天）内接触的平均人数为常数$\lambda$，$\lambda$称为日感染率，当病人与健康者有效接触时，可使健康者受感染成为病人。
（5）每个病人得病后，经久不愈，且在传染期内不会死亡。

模型的建立与求解
根据假设（4），每个病人每天可使$\lambda s(t)$个健康者变为病人，而$t$时刻病人总数为$Ni(t)$，故在$t\to t+\Delta t$时段内，共有$\lambda Ns(t)i(t)\Delta t$个健康者被感染。
于是有
$$\frac{Ni(t+\Delta t)-Ni(t)}{\Delta t}=\lambda Ns(t)i(t).$$
令$\Delta t\to 0$，得微分方程
$$\frac{di(t)}{dt}=\lambda s(t)i(t).$$

又由假设（3）知，$s(t)=1-i(t)$，代入上式得
$$\frac{di(t)}{dt}=\lambda i(t)(1-i(t)).$$
假定起始时（$t=0$），病人占总人口的比例为$i(0)=i_0$。于是SI模型可描述为
$$\{\begin{array}{ll}\frac{di(t)}{dt}=\lambda i(t)(1-i(t)),& t>0\\ i(0)=i_0.& \end{array}(6.14)$$
用分离变量法求解此微分方程初值问题，得解析解为
$$i(t)=\frac{1}{1+\left(\frac{1}{i_0}-1\right)e^{-\lambda t}}.(6.15)$$

结果分析与评价
模型(6.14)事实上就是Logistic模型。病人占总人口的最大比例为1，即当$t\to \infty$时，区域内所有人都被传染。
医学上称$\frac{di}{dt}$~$t$为传染病曲线，它表示传染病人增加率与时间的关系，见图6.2（A）所示。预测结果曲线见图6.2（B）所示。

由模型(6.14)易知，当病人总量占总人口比值达到$i=\frac{1}{2}$时， 达到最大值，即$\frac{d^{2}i}{d{t}^2}=0$，也就是说，此时达到传染病传染高峰期。利用(6.15)易得传染病高峰到来的时刻为
$$t_m=\frac{1}{\lambda }ln\left(\frac{1}{i_0}-1\right).$$

医学上，这一结果具有重要的意义。由于$t_m$与$\lambda$成反比，故当$\lambda$（反应医疗水平或传染控制措施的有效性）增大时，$t_m$将变小，预示着传染病高峰期来得越早。若已知日接触率$\lambda$（由统计数据得出），即可预报传染病高峰到来的时间$t_m$，这对于防治传染病是有益处的。

当$t\to \infty$时，由(6.15)式可知，$i(t)\to 1$，即最后人人都要生病。这显然是不符合实际情况的。其原因是假设中未考虑病人得病后可以治愈，人群中的健康者只能变为病人，而病人不会变为健康者。而事实上对某些传染病，如伤风、痢疾等病人治愈后免疫力低下，可假定无免疫性。于是病人被治愈后成为健康者，健康者还可以被感染再变成病人。

3. SIS模型

SIS模型在SI模型假设的基础上，进一步假设：
（6）每天被治愈的病人人数占病人总数的比例为$\mu$。
（7）病人被治愈后成为仍可被感染的健康者。

于是SI模型可被修正为SIS模型：
$$\{\begin{array}{ll}\frac{di(t)}{dt}=\lambda i(t)(1-i(t))-\mu i(t),& t>0\\ i(0)=i_0.& \end{array}(6.16)$$

模型(6.16)的解析解可表示为
$$\{\begin{array}{ll}[\frac{\lambda }{\lambda -\mu }+(\frac{1}{i_0}-\frac{\lambda }{\lambda -\mu })e^{-(\lambda -\mu )t}{]}^{-1},& \lambda \ne \mu \\ (\lambda t+\frac{1}{i_0}{)}^{-1},& \lambda =\mu \end{array}(6.17)$$

若令
$$\sigma =\frac{\lambda }{\mu }，$$
称$\sigma$为传染强度。

利用$\sigma$的定义，方程(6.16)可改写为
$$\{\begin{array}{ll}\frac{di}{dt}=-\lambda i[i-(1-\frac{1}{\sigma })],& t>0,\\ i(0)=i_0.& \end{array}(6.18)$$
相应地，模型的解析解可表示为
$$i(t)=\{\begin{array}{ll}[\frac{1}{1-\frac{1}{\sigma }}+(\frac{1}{i_0}-\frac{1}{1-\frac{1}{\sigma }})e^{-\lambda (1-\frac{1}{\sigma })t}{]}^{-1},& \sigma \ne 1，\\ (\lambda t+\frac{1}{i_0}{)}^{-1}& \sigma =1\end{array}(6.19)$$

结果分析与评价
由(6.19)得，当$t\to \infty$时，有
$$i(t)=\{\begin{array}{ll}1-\frac{1}{\sigma },& \sigma >1，\\ 0,& \sigma \le 1\end{array}(6.20)$$
由上式可知，$\sigma =1$是一个阈值。

若$\sigma \le 1$，随着时间的推移，$i(t)$逐渐变小，当$t\to \infty$时趋于零。这是由于治愈率大于有效感染率，最终所有病人都会被治愈。

若$\sigma >1$，则当$t\to \infty$时，$i(t)$趋于极限$1-\frac{1}{\sigma }$，这说明当治愈率小于传染率时，总人口中总有一定比例的人口会被传染而成为病人。
大多数传染病，如天花、麻疹、流感、肝炎等疾病经治愈后均有很强的免疫力。病愈后的人因已具有免疫力，既非健康者（易感染者）也非病人（已感染者），即这部分人已退出感染系统。

4. SIR模型

基本假设
（1）人群分健康者、病人和病愈后因具有免疫力而退出系统的移出者三类。设任意时刻$t$，这三类人群占总人口的比例分别为：$s(t)$,$i(t)$和$r(t)$。
（2）病人的日接触率为 ，日治愈率为$\mu$，传染强度$\sigma =\lambda /\mu$。
（3）人口总数$N$为固定常数。

模型建立
类似于前述问题的建模过程，依据假设，有
对所有人群， $s(t)+i(t)+r(t)=1$； (6.21)
对系统移出者， $N\frac{dr}{dt}=\mu Ni$； (6.22)
对病人， $N\frac{di}{dt}=\lambda Nsi-\mu Ni$； (6.23)
对健康者， $N\frac{ds}{dt}=-\lambda Nsi$. (6.24)

联立(6.22)- (6.24)，可得SIR模型：
$$\{\begin{array}{l}\frac{di}{dt}=\lambda si-\mu i,\\ \frac{ds}{dt}=\lambda si,\\ \frac{dr}{dt}=\mu i,\\ i(0)=i_0,s(0)=s_0,r(0)=0.\end{array}(6.25)$$
SIR模型是一个较典型的系统动力学模型，其突出特点是模型形式为关于多个相互关联的系统变量之间的常微分方程组。类似的建模问题有很多，如河流中水体各类污染物质的耗氧、复氧、反应、迁移、吸附、沉降等，食物在人体中的分解、吸收、排泄，污水处理过程中的污染物降解，微生物、细菌增长或衰减等。这些问题很难求得解析解，可以使用软件求数值解。

---

三、常微分方程的求解

对于常微分方程，只有一小部分可以求得解析解，大部分常微分方程是无法求得解析解，只能求数值解。常微分方程数值解的算法我们就不介绍了，有兴趣的读者可以参看数值分析等相关书籍。下面介绍使用Matlab软件求微分方程的符号解和数值解。

1. 常微分方程的符号

Matlab符号运算工具箱提供了功能强大的求常微分方程符号解函数dsolve。

【例3】 求解微分方程$y^{\prime }=-2y+2x^2+2x$，$y(0)=1$。
【解】
求得的符号解为$y=e^{-2x}+x^2.$

 clc, clear, syms y(x)
y=dsolve(diff(y)==-2*y+2*x^2+2*x, y(0)==1)

【例4】 求解二阶线性微分方程$y^{\prime \prime }-2y^{\prime }+y=e^x$，$y(0)=1$，$y^{\prime }(0)=-1$。
【解】
求得二阶微分方程的解为$y=e^x+\frac{x^2{e}^x}{2}-2x{e}^x$。

clc, clear, syms y(x)  %定义符号函数ｙ，自变量为x
dy=diff(y);  %定义y的一阶导数，目的是下面赋初值
y=dsolve(diff(y,2)-2*diff(y)+y==exp(x), y(0)==1, dy(0)==-1)

【例5】 已知输入信号为$u(t)=e^{-t}cos(t)$，试求下面微分方程的解。
$y^{(4)}(t)+10y^{(3)}(t)+35y^{\prime \prime }(t)+50y^{\prime }(t)+24y(t)=u^{\prime \prime }(t).$
$y(0)=0，y^{\prime }(0)=-1，y^{\prime \prime }(0)=1，y^{\prime \prime \prime }(0)=1.$
【解】
求得的解为$y=-\frac{7}{3}{e}^{-t}+\frac{9}{2}{e}^{-2t}-\frac{14}{5}{e}^{-3t}-\frac{e^{-t}sint}{5}+\frac{19}{30}{e}^{-4t}$.

clc, clear, syms y(t)
dy=diff(y); d2y=diff(y,2); d3y=diff(y,3); %定义y的前3阶导数，是为了赋初值
u=exp(-t)*cos(t);
y=dsolve(diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==diff(u,2),...
    y(0)==0,dy(0)==-1,d2y(0)==1,d3y(0)==1)

【例6】 试求解下列Cauchy问题
$$\{\begin{array}{l}\frac{dx}{dt}=Ax,\\ x(0)=[1,1,1{]}^T.\end{array}$$
的解，其中$\mathbf{x}(t)=[x_1(t),x_2(t),x_3(t){]}^T$，$\mathbf{A}=\left[\begin{array}{ccc}3& -1& 1\\ 2& 0& -1\\ 1& -1& 2\end{array}\right].$
【解】
求得的解为
$$\mathbf{x}(t)=\left[\begin{array}{c}\frac{4}{3}{e}^{3t}-\frac{1}{2}{e}^{2t}+\frac{1}{6}\\ \frac{2}{3}{e}^{3t}-\frac{1}{2}{e}^{2t}+\frac{5}{6}\\ \frac{2}{3}{e}^{3t}+\frac{1}{3}\end{array}\right].$$

clc, clear
syms x(t)  [3, 1]  %定义符号向量函数，x(t)后面有空格
A=[3 -1 1;2 0 -1;1 -1 2];
[s1,s2,s3]=dsolve(diff(x)==A*x,x(0)==ones(3,1))

代码运行结果：

【例7】 试解初值问题

$\left[\begin{array}{c}{x_1}^{\prime }(t)\\ {x_2}^{\prime }(t)\\ {x_3}^{\prime }(t)\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& -2\\ 3& 2& 1\end{array}\right]\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ e^{t}cos2t\end{array}\right]，\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ 1\end{array}\right].$

【解】
求得的解为
$$\mathbf{x}(t)=\left[\begin{array}{c}\frac{4}{3}{e}^{3t}-\frac{1}{2}{e}^{2t}+\frac{1}{6}\\ \frac{2}{3}{e}^{3t}-\frac{1}{2}{e}^{2t}+\frac{5}{6}\\ \frac{2}{3}{e}^{3t}+\frac{1}{3}\end{array}\right].$$
$$\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right]=\left[\begin{array}{c}0\\ e^t\left[cos(2t)-sin(2t)-\frac{tsin(2t)}{2}\right]\\ e^t\left[cos(2t)+\frac{5sin(2t)}{4}+\frac{tcos(2t)}{2}\right]\end{array}\right]$$

clc,clear
syms x(t)  [3,1]  %定义符号向量函数，x(t)后面有空格
A=[1,0,0;2,1,-2;3,2,1]; B=[0;0;exp(t)*cos(2*t)];
x0=[0;1;1]; %初值条件
[s1,s2,s3]=dsolve(diff(x)==A*x+B,x(0)==x0) %求符号解

【例8】 是求常微分方程
$$\{\begin{array}{l}{f}^{{}^{\prime \prime }}+g=3\\ g^{{}^{\prime }}+f^{{}^{\prime }}=1\end{array}$$
在初边值条件$f^{{}^{\prime }}(1)=0,f(0)=0,g(0)=0$时的解。
【解】
$$f(x)=x-\frac{e-3}{e^2+1}{e}^x+\frac{e(3e+1)}{e^2+1}{e}^{-x}-3，$$
$$g(x)=\frac{e-3}{e^2+1}{e}^x-\frac{3e^2+e}{e^2+1}{e}^{-x}+3.$$

clc,clear, syms f(x) g(x)  %定义符号函数
df=diff(f);  %定义f的一阶导数
[s1,s2]=dsolve(diff(f,2)+g==3, diff(g)+diff(f)==1,...
        df(1)==0, f(0)==0, g(0)==0)

2. 初值问题的Matlab数值解

Matlab的工具箱提供了几个解常微分方程数值解的函数，如ode45，ode23，ode113，其中ode45采用四五阶龙格库塔方法（以下简称RK方法），是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。

在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组
$$\frac{d\mathbf{y}}{dx}=\mathbf{A}\mathbf{y}+\mathbf{\Phi }(x)，(6.26)$$
其中$\mathbf{y},\mathbf{\Phi }\in R^m,\mathbf{A}$为$m$阶方阵。

若矩阵$\mathbf{A}$的特征值${\lambda }_i(i=1,2,\cdots ,m)$满足关系

则称方程组(6.26)式为刚性方程组或Stiff方程组，称数

为刚性比。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长$h$不作任何限制。

Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s，ode23t，ode23tb。
这里只简单介绍ode45求数值解的用法。

对一阶微分方程或方程组的初值问题
$$\{\begin{array}{l}{y}^{\prime }=f(t,y),\\ y(t_0)=y_0,\end{array}$$
其中$y$和$f$可以为向量.

函数ode45有如下两种调用格式：

[t,y]=ode45(fun,tspan,y0) 或 s=ode45(fun,tspan,y0)

其中$fun$是用$M$函数或匿名函数定义$f(t,y)$的函数文件名或匿名函数返回值，
$tspan=[t_0,tfinal]$(这里t_0必须是初值条件中自变量的取值，$tfinal$可以比$t_0$小)
是求解区间，$y_0$是初值。返回值$t$是MATLAB自动离散化的区间$[t_0,tfinal]$上的点，$y$的列是对应于$t$的函数值；如果只有一个返回值$s$，则$s$是一个结构数组。

利用结构数组$s$和MATLAB函数$deval$，我们可以计算区间$tspan$中任意点$x$的函数值，调用格式为
$$y=deval(s,x)，$$
其中$x$为标量或向量，返回值$y$的行是对应于$x$的数值解。

【例9】 是求常微分方程
【解】 （续例3）求微分方程$y^{\prime }=-2y+2x^2+2x$，$y(0)=1$ ，$0\le x\le 0.5$的数值解；并在同一个图形界面上画出数值解和符号解的曲线。

clc, clear, close all, syms y(x)
y=dsolve(diff(y)==-2*y+2*x^2+2*x, y(0)==1)
dy=@(x,y)-2*y+2*x^2+2*x;
[sx, sy]=ode45(dy, [0,0.5], 1)
fplot(y,[0,0.5]), hold on
plot(sx, sy, '*'); legend({'符号解','数值解'})
xlabel('$x$','Interpreter','Latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)

代码运行结果：

Matlab无法直接求解高阶微分方程或方程组的数值解，必须化成一阶微分方程组，才能求数值解。

【例10】（续例2）求二阶常微分方程(6.10)式的数值解。
【解】
求数值解时，需要把二阶微分方程转化为一阶微分方程组，引进$y_1=y$，$y_2=y^{\prime }$，则方程(6.10)式可以转化为如下的一阶微分方程组：
$$\{\begin{array}{ll}{y}_1^{{}^{\prime }}=y_2,& y_1(0)=0,\\ y_2^{{}^{\prime }}=\frac{1}{5(1-x)}\sqrt{1+y_2^2},& y_2(0)=0.\end{array}$$
最后得到的导弹轨迹曲线如图6.3所示。

clc, clear, close all
dy=@(x,y)[y(2); sqrt(1+y(2)^2)/5/(1-x)];
[x,y]=ode45(dy,[0,0.999999],[0;0])
plot(x, y(:,1)), xlabel('$x$','Interpreter','Latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)

【例11】 Lorenz模型的混沌效应。
【解】
Lorenz模型是由美国气象学家Lorenz在研究大气运动时，通过简化对流模型，只保留3个变量提出的一个完全确定性的一阶自治常微分方程组（不显含时间变量），其方程为
$$\{\begin{array}{l}x=\sigma (y-x),\\ y=\rho x-y-xz,\\ z=xy-\beta z.\end{array}$$
其中，参数$\sigma$为Prandtl数，$\rho$为Rayleigh数，$\beta$为方向比。

Lorenz模型如今已经成为混沌领域的经典模型，第一个混沌吸引子—Lorenz吸引子也是在这个系统中被发现的。系统中三个参数的选择对系统会不会进入混沌状态起着重要的作用。图6.4（A）给出了Lorenz模型在$\sigma =10$，$\rho =28$，$\beta =\frac{8}{3}$时系统的三维演化轨迹。由图6.4（A）可见，经过长时间运行后，系统只在三维空间的一个有限区域内运动，即在三维相空间里的测度为零。图6.4（A）显示出我们经常听到的“蝴蝶效应”。

图6.4（B）给出了系统从两个靠的很近的初值出发（相差仅0.00001）后，解的偏差演化曲线。随着时间的增大，可以看到两个解的差异越来越大，这正是动力学系统对初值敏感性的直观表现，由此可断定此系统的这种状态为混沌态。混沌运动是确定性系统中存在随机性，它的运动轨道对初始条件极端敏感。

clc, clear, close all, rng(2)  %为了进行一致性比较，每次运行取相同随机数
sigma=10; rho=28; beta=8/3; T=80;
g=@(t,f)[sigma*(f(2)-f(1)); rho*f(1)-f(2)-f(1)*f(3); f(1)*f(2)-beta*f(3)]; 
xyz0=rand(3,1); %初始值
[t,xyz]=ode45(g,[0,T],xyz0); %求数值解
subplot(121), plot3(xyz(:,1),xyz(:,2),xyz(:,3)) %轨线图
xlabel('$x(t)$','Interpreter','latex')
ylabel('$y(t)$','Interpreter','latex') 
zlabel('$z(t)$','Interpreter','latex'), box on %加盒子线以突出立体感
so=ode45(g,[0,T],xyz0+0.00001) %初值变化后，再求数值解
xyz2=deval(so,t); %返回值为3行的矩阵，计算对应的x,y,z的值
subplot(122), plot(t,xyz(:,1)-xyz2(1,:)','.-') 
xlabel('$x(t)$','Interpreter','latex')
ylabel('$x_1(t)-x_2(t)$','Interpreter','latex')

所画出的图形如图6.4所示。

【例12】 一个慢跑者在平面上按如下规律跑步
$X=10+20cost$，$Y=20+15sint.$
突然有一只狗攻击他，这只狗从原点出发，以恒定速率 跑向慢跑者，狗运动方向始终指向慢跑者。分别求出$w=20$，$w=5$时狗的运动轨迹。

【解】
设时刻$t$人的坐标为(X(t),Y(t))，狗的坐标$(x(t),y(t))$。狗的速度大小恒为$w$，则
${\left(\frac{dx}{dt}\right)}^2+{\left(\frac{dy}{dt}\right)}^2=w^2，$

由于狗始终对准人，故狗的速度方向平行于狗与人位置的差向量，
$\left[\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right]=\lambda \left[\begin{array}{c}X-x\\ Y-y\end{array}\right],\lambda >0$

消去$\lambda$得
$$\{\begin{array}{l}\frac{dx}{dt}=\frac{w}{\sqrt{(X-x{)}^2+(Y-y{)}^2}}(X-x),\\ \frac{dy}{dt}=\frac{w}{\sqrt{(X-x{)}^2+(Y-y{)}^2}}(Y-y).\end{array}$$

因而建立狗的运动轨迹的如下方程
$$\{\begin{array}{l}\frac{dx}{dt}=\frac{w}{\sqrt{(10+20cost-x{)}^2+(20+15sint-y{)}^2}}(10+20cost-x),\\ \frac{dy}{dt}=\frac{w}{\sqrt{(10+20cost-x{)}^2+(20+15sint-y{)}^2}}(20+15sint-y),\\ x(0)=0,y(0)=0.\end{array}$$

利用Matlab软件求得当$w=20$时，在$t=4.1888$时，狗追上人。人和狗之间的距离见图6.5（A）。
当$w=5$时，狗永远追不上人。人和狗之间的距离见图6.5（B）。

clc, clear, close all, w=20;
df=@(t,f,w)[w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(10+20*cos(t)-f(1))
    w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(20+15*sin(t)-f(2))];
t0=0; tf=5; %时间的终值是适当估计的
s1=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0]); %求微分方程的数值解
d1=@(t)sqrt((deval(s1,t,1)-10-20*cos(t)).^2+...
    (deval(s1,t,2)-20-15*sin(t)).^2); %t时刻人和狗的距离
[tmin,fmin]=fminbnd(d1,0,tf) %求两者距离的最小值及对应的时间
t=0:0.1:tf; subplot(121), plot(t,d1(t))  %画出两者之间的距离
xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')
w=5; tf=60; 
[t,s]=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0]); %求微分方程的数值解
d2=sqrt((s(:,1)-10-20*cos(t)).^2+(s(:,2)-20-15*sin(t)).^2); %计算两者距离
subplot(122), plot(t,d2)  %画出两者之间的距离
xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')

【例13】 求解二阶线性微分方程$y^{\prime \prime }-2y^{\prime }+y=e^x$，$y(0)=1$ ，$y^{\prime }(0)=-1$，在区间$[-1,1]$上的数值解。

【解】
求高价微分方程数值解时，首先做变量替换化成一阶微分方程组。设y_1=y，y_2=y’，则把上述二阶微分方程化成如下的一阶微分方程组：
$$\{\begin{array}{ll}{y}_1^{{}^{\prime }}=y_2,& y_1(0)=1,\\ y_2{}^{\prime }=2y_2-y_1+e^x,& y_2(0)=-1.\end{array}$$

所求得的数值解如图6.6所示。

clc, clear, close all
df=@(x,f)[f(2); 2*f(2)-f(1)+exp(x)]; hold on
s1=ode45(df,[0,1],[1;-1]); %求[0,1]区间上的数值解
s2=ode45(df,[0,-1],[1;-1]); %求[-1,0]区间上的数值解
fplot(@(x)deval(s1,x,1),[0,1],'-ok') %画第一个分量y在[0,1]区间上的数值解
fplot(@(x)deval(s2,x,1),[-1,0],'-*k') %画第一个分量y在[-1,0]区间上的数值解
xlabel('$x$', 'Interpreter', 'latex')
ylabel('$y$', 'Interpreter', 'latex', 'Rotation', 0)

【例14】 已知阿波罗卫星的运动轨迹$(x,y)$满足下面方程：
$$\{\begin{array}{l}\frac{d^{2}x}{d{t}^2}=2\frac{dy}{dt}+x-\frac{\lambda (x+\mu )}{r_1^3}-\frac{\mu (x-\lambda )}{r_2^3},\\ \frac{d^{2}y}{d{t}^2}=-2\frac{dx}{dt}+y-\frac{\lambda y}{r_1^3}-\frac{\mu y}{r_2^3}.\end{array}$$
其中$\mu =1/82.45$，$\lambda =1-\mu$，$r_1=\sqrt{(x+\mu {)}^2+y^2}$，$r_2=\sqrt{(x+\lambda {)}^2+y^2}$，试在初值$x(0)=1.2$，$x^{\prime }(0)=0$，$y(0)=0$，$y^{\prime }(0)=-1.0494$下求解，并绘制阿波罗卫星轨迹图。

【解】
做变量替换，令$z_1=x$，$z_2=\frac{dx}{dt}$，$z_3=y$，$z_4=\frac{dy}{dt}$，则原二阶微分方程组可以化为如下的一阶方程组
$$\{\begin{array}{ll}\frac{d{z}_1}{dt}=z_2,& z_1(0)=1.2,\\ \frac{d{z}_2}{dt}=2z_4+z_1-\frac{\lambda (z_1+\mu )}{((z_1+\mu {)}^2+z_3^2{)}^{3/2}}-\frac{\mu (z_1-\lambda )}{((z_1+\lambda {)}^2+z_3^2{)}^{3/2}},& z_2(0)=0，\\ \frac{d{z}_3}{dt}=z_4,& z_3(0)=0,\\ \frac{d{z}_4}{dt}=-2z_2+z_3-\frac{\lambda z_3}{((z_1+\mu {)}^2+z_3^2{)}^{3/2}}-\frac{\mu z_3}{((z_1+\lambda {)}^2+z_3^2{)}^{3/2}},z_4(0)=-1.0494.& \end{array}$$

clc, clear, close all
mu=1/82.45; lamda=1-mu;
dz=@(t,z)[z(2); 2*z(4)+z(1)-lamda*(z(1)+mu)/((z(1)+mu)^2+z(3)^2)^(3/2)-...
          mu*(z(1)-lamda)/((z(1)+lamda)^2+z(3)^2)^(3/2)
          z(4); -2*z(2)+z(3)-lamda*z(3)/((z(1)+mu)^2+z(3)^2)^(3/2)-...
          mu*z(3)/((z(1)+lamda)^2+z(3)^2)^(3/2)];
[t,z]=ode45(dz,[0,100],[1.2 0 0 -1.0494])
plot(z(:,1), z(:,3), 'k')  %画轨迹图
xlabel('$x$', 'Interpreter', 'latex')
ylabel('$y$', 'Interpreter', 'latex', 'Rotation', 0)

3. 边值问题的Matlab数值解

Matlab中用bvp4c和bvpinit命令求解常微分方程的两点边值问题，微分方程的标准形式为
$y^{\prime }=f(x,y)$，$bc(y(a),y(b))=0$，
或者是
$y^{\prime }=f(x,y,p)$，$bc(y(a),y(b),p)=0$，
式中，$p$是有关的参数。这里$y,f$可以为向量函数，求解的区间为$[a,b]$，$bc$为边界条件。

【例15】 考察描述在水平面上一个小水滴横截面形状的标量方程
$\frac{d^2}{d{x}^2}h(x)+(1-h(x))(1+(\frac{d}{dx}h(x){)}^2{)}^{3/2}=0，h(-1)=h(1)=0，$
这里$h(x)$表示$x$处水滴的高度。
【解】
设$y_1(x)=h(x)，y_2(x)=\frac{dh(x)}{dx}$，把上述微分方程写成两个一阶微分方程组
$$\{\begin{array}{l}\frac{d}{dx}{y}_1(x)=y_2(x),\\ \frac{d}{dx}{y}_2(x)=(y_1(x)-1)(1+y_2^2(x){)}^{3/2}.\end{array}$$

上述微分方程组可以由如下函数表示
function yprime=drop(x,y);
yprime=[y(2);(y(1)-1)*(1+y(2)^ 2)^(3/2)];
边界条件通过残差函数指定，边界条件通过如下函数表示
function res=dropbc(ya,yb);
res=[ya(1); yb(1)];

使用$y_1(x)=\sqrt{1-x^2}$和$y_2(x)=-\frac{x}{\left(0.1+\sqrt{1-x^2}\right)}$（这里分母加0.1是为了避免奇性）作为初始猜测解（初始解可以是任意取的，如取$y_1(x)=x^2$和$y_2(x)=2x$），由如下函数定义

function yinit=dropinit(x);
yinit=[sqrt(1-x.^2); -x./(0.1+sqrt(1-x.^2))];

利用如下的程序就可以求微分方程的边值问题并画出图6.8。

solinit=bvpinit(linspace(-1,1,20), @dropinit);
sol=bvp4c(@drop, @dropbc,solinit);
fill(sol.x, sol.y(1,:), [0.7,0.7,0.7]), axis([-1,1,0,1])
xlabel('$x$','Interpreter', 'latex', 'FontSize',12)
ylabel('$h$','Interpreter', 'latex', 'Rotation', 0, 'FontSize',12)

这里调用函数bvpinit，计算区间$[-1,1]$上等间距的20个点的数据，然后调用函数bvp4c，得到数值解的结构sol，填充命令fill填充$x-y_1$平面上的解曲线。

一般地，bvp4c的调用格式如下：

sol=bvp4c(@odefun,@bcfun,solinit,options,p1,p2,…);

函数odefun的格式为

yprime=odefun(x,y,p1,p2, …);

函数bcfun的格式为

res=bcfun(ya,yb, p1,p2, …);

初始猜测结构solinit有两个域，solinit.x提供初始猜测的$x$值，排列顺序从左到右排列，其中solinit.x(1)和solinit.x(end)分别为$a$和$b$。对应地，solinit.y(:,i)给出点solinit.x(i)处初始猜测解。

输出参数sol是包含数值解的一个结构，其中sol.x给出了计算数值解的 点，sol.x(i)处的数值解由sol.y(:,i)给出，类似地，sol.x(i)处数值解的一阶导数值由sol.yp(:,i)给出。
可以把上面的所有函数都放在一个文件中，程序如下：

clc, clear, close all
solinit=bvpinit(linspace(-1,1,20), @dropinit);
sol=bvp4c(@drop, @dropbc, solinit)
fill(sol.x, sol.y(1,:), [0.7,0.7,0.7]), axis([-1,1,0,1])
xlabel('$x$','Interpreter', 'latex', 'FontSize',12)
ylabel('$h$','Interpreter', 'latex', 'Rotation', 0, 'FontSize',12)
 

function yprime=drop(x,y);
yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)]; end
 
function res=dropbc(ya,yb);
res=[ya(1);yb(1)]; end
 
function yinit=dropinit(x);
yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))]; end

【例16】 描述x=0处固定，x=1处有弹性支持，沿着x轴平衡位置以均匀角速度旋转的绳的位移方程
$\frac{d^2}{d{x}^2}y(x)+\mu y(x)=0，$
具有边界条件
$y(0)=0，{(\frac{d}{dx}y(x))|}_{x=0}=1，{(y(x)+\frac{d}{dx}y(x))|}_{x=1}=0.$

【解】
这个边值问题是一个特征问题，必须找到参数\mu的值使得方程的解存在。如果提供了参数\mu的猜测值和对应解的猜测值，也可以利用函数bvp4c求解特征问题。上述微分方程可以写成下面的微分方程组
$$\{\begin{array}{l}\frac{d}{dx}{y}_1(x)=y_2(x),\\ \frac{d}{dx}{y}_2(x)=-\mu y_1(x).\end{array}$$
使用$\mu =5，y_1(x)=sinx$和$y_2(x)=cosx$作为初始猜测解。编写程序如下：

clc, clear
eq=@(x,y,mu)[y(2);-mu*y(1)]; %定义一阶方程组的匿名函数
bd=@(ya,yb,mu)[ya(1);ya(2)-1;yb(1)+yb(2)];  %定义边值条件的匿名函数
guess=@(x)[sin(2*x);2*cos(2*x)]; %定义初始猜测解的匿名函数
guess_structure=bvpinit(linspace(0,1,10),guess,5); %给出初始猜测解的结构,mu=5
sol=bvp4c(eq,bd,guess_structure); %计算数值解
plot(sol.x, sol.y(1,:), '-', sol.x, sol.yp(1,:), '--', 'LineWidth',2)
xlabel('$x$','FontSize',12,'Interpreter','Latex')
ylabel('$y$','FontSize',12,'Interpreter','Latex','Rotation',0)
legend('$y_1$','$y_2$','Interpreter','Latex')

【例17】 微分方程组为
$$\{\begin{array}{l}{u}^{{}^{\prime }}=0.5u(w-u)/v,\\ v^{{}^{\prime }}=-0.5(w-u),\\ w^{\prime }=(0.9-1000(w-y)-0.5w(w-u))/z,\\ z^{{}^{\prime }}=0.5(w-u),\\ y^{{}^{\prime }}=-100(y-w),\end{array}$$
边界条件为$u(0)=v(0)=w(0)=1$，$z(0)=-10$；$w(1)=y(1)$。

【解】
使用如下猜测解：
$$\{\begin{array}{l}u(x)=1,\\ v(x)=1,\\ w(x)=-4.5x^2+8.91x+1,\\ z(x)=10,\\ y(x)=-4.5x^2+9x+0.91.\end{array}$$
编写如下程序：

clc, clear
eq=@(x,y)[0.5*y(1)*(y(3)-y(1))/y(2)
          -0.5*(y(3)-y(1))
          (0.9-1000*(y(3)-y(5))-0.5*y(3)*(y(3)-y(1)))/y(4)
          0.5*(y(3)-y(1))
          100*(y(3)-y(5))]; %定义一阶方程组的匿名函数
bd=@(ya,yb)[ya(1)-1;ya(2)-1;ya(3)-1;ya(4)+10;yb(3)-yb(5)]; %定义边值条件
guess=@(x)[1;1;-4.5*x.^2+8.91*x+1;-10;-4.5*x.^2+9*x+0.91]; %定义初始猜测解
guess_structure=bvpinit(linspace(0,1,5),guess); %给出初始猜测解的结构
sol=bvp4c(eq,bd,guess_structure); %计算数值解
plot(sol.x,sol.y(1,:),'-*',sol.x,sol.y(2,:),'-D',sol.x,sol.y(3,:),':S',...
    sol.x,sol.y(4,:),'-.O',sol.x,sol.y(5,:),'--P') %画出5条解曲线
legend({'$u$','$v$','$w$','$z$','$y$'}, 'Interpreter', 'latex', ...
       'Location','southwest')  %图例标注在左下角

---

四、微分方程建模实例

1. Malthus模型（马尔萨斯模型）

1789年，英国神父Malthus在分析了一百多年人口统计资料之后，提出了Malthus模型。

1. 模型假设

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

2. 建模与求解

由假设，$t$时刻到$t+\Delta t$时刻人口的增量为$x(t+\Delta t)-x(t)=rx(t)\Delta t$，于是得
$$\{\begin{array}{l}\frac{dx}{dt}=rx,\\ x(0)=x_0,\end{array}(6.27)$$
其解为
$$x(t)=x_0{e}^{rt}.(6.28)$$

3. 模型评价

考虑二百多年来人口增长的实际情况，1961年世界人口总数为$3.06\times 10^9$，在1961～1970年这段时间内，每年平均的人口自然增长率为2%，则(6.28)式可写为
$$x(t)=3.06\times 10^9\cdot e^{0.02(t-1961)}.(6.29)$$
根据1700～1961年间世界人口统计数据，发现这些数据与(6.29)式的计算结果相当符合。因为在这期间地球上人口大约每 35年增加 1倍，而(6.29)式算出每 34.6年增加1倍。

但是，利用(6.29)式对世界人口进行预测，也会得出惊异的结论，当$t=2670$年时，$x(t)=4.4\times 10^{15}$，即4400万亿，这相当于地球上每平方米要容纳至少 20人。
显然，用这一模型进行预测的结果远高于实际人口增长，误差的原因是对增长率$r$的估计过高。由此，可以对$r$是常数的假设提出疑问。

2. Logistic模型

如何对增长率r进行修正呢？我们知道，地球上的资源是有限的，它只能提供一定数量的生命生存所需的条件。随着人口数量的增加，自然资源、环境条件等对人口再增长的限制作用将越来越显著。如果在人口较少时，可以把增长率$r$看成常数，那么当人口增加到一定数量之后，就应当视r为一个随着人口的增加而减小的量，即将增长率r表示为人口$x(t)$的函数$r(x)$，且$r(x)$为$x$的减函数。

（1）设$r(x)$为$x$的线性函数，$r(x)=r-sx$（工程师原则，首先用线性）。
（2）自然资源与环境条件所能容纳的最大人口数为$x_m$，即当$x=x_m$时，增长率$r(x_m)=0$。

由假设（1），（2）可得$r(x)=r(1-\frac{x}{x_m})$，则有
$$\{\begin{array}{l}\frac{dx}{dt}=r(1-\frac{x}{x_m})x,\\ x(t_0)=x_0.\end{array}(6.30)$$
(6.30)式是一个可分离变量的方程，其解为
$$x(t)=\frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r(t-t_0)}}(6.31)$$

由(6.30)式，计算可得
$$\frac{d^{2}x}{d{t}^2}=r^2(1-\frac{x}{x_m})(1-\frac{2x}{x_m})x.(6.32)$$
人口总数x(t)有如下规律：
（1）${\lim }_{x\to \infty }x(t)=x_m$ ，即无论人口初值$x_0$如何，人口总数以$x_m$为极限。
（2）当$0<x_0<x_m$时，$\frac{dx}{dt}=r(1-\frac{x}{x_m})x>0$，这说明$x(t)$是单调增加的，又由(6.32)式知，当$x<\frac{x_m}{2}$时，$\frac{d^{2}x}{d{t}^2}>0$，$x=x(t)$为凹函数，当$x>\frac{x_m}{2}$时，$\frac{d^{2}x}{d{t}^2}<0$，$x=x(t)$为凸函数。
（3）人口变化率$\frac{dx}{dt}$在$x=\frac{x_m}{2}$时取到最大值，即人口总数达到极限值一半以前是加速生长时期，经过这一点之后，增长速率会逐渐变小，最终达到零。

3. 美国人口的预报模型

【例18】
利用表6.1给出的近两个世纪的美国人口统计数据（以百万为单位），建立人口预测模型，最后用它预报2010年美国的人口。

建模与求解

记$x(t)$为第$t$年的人口数量，设人口年增长率$r(x)$为$x$的线性函数，$r(x)=r-sx$。自然资源与环境条件所能容纳的最大人口数为$x_m$，即当$x=x_m$时，增长率$r(x_m)=0$，可得$r(x)=r(1-\frac{x}{x_m})$，建立Logistic人口模型
$$\{\begin{array}{l}\frac{dx}{dt}=r(1-\frac{x}{x_m})x,\\ x(t_0)=x_0,\end{array}$$

其解为
$$x(t)=\frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r(t-t_0)}}.(6.33)$$

参数估计
把表6.1中的全部数据保存到文本文件data6_18.txt中。
（1）非线性最小二乘估计
把表6.1中的第1个数据作为初始条件，利用余下的数据拟合式(6.33)中的参数$x_m$和$r$，求得$r=0.0274$，$x_m=342.4441$，2010年人口的预测值为282.6806百万。
（2）线性最小二乘估计
为了利用简单的线性最小二乘法估计这个模型的参数$r$和$x_m$，把Logistic方程表示为
$$\frac{1}{x}\cdot \frac{dx}{dt}=r-sx,s=\frac{r}{x_m}(6.34)$$

分别用$k=1,2,\cdots ,22$表示$1790,1800,\cdots ,2000$年。
1）利用向前差分，得到差分方程
$$\frac{1}{x(k)}\frac{x(k+1)-x(k)}{\Delta t}=r-sx(k),k=1,2,\cdots ,21，(6.35)$$
其中步长$\Delta t=10$。

利用Matlab软件求得r=0.0325，$x_m=294.3860$。2010年人口的预测值为$277.9634$百万。
2）利用向后差分，得到差分方程
$$\frac{1}{x(k)}\frac{x(k)-x(k-1)}{\Delta t}=r-sx(k),k=2,3,\cdots ,22，(6.36)$$
其中步长$\Delta t=10$。

利用Matlab软件求得$r=0.0247$，$x_m=373.5135$。2010年人口的预测值为$264.9119$百万。
从上面的三种拟合方法可以看出，拟合同样的参数，方法不同可能结果相差很大。

clc, clear, a=readmatrix('data6_18.txt');
x=a([2:2:6],:)'; x=x(~isnan(x));
fn=@(r,xm,t)xm./(1+(xm/3.9-1)*exp(-r*(t-1790))); %定义匿名函数
ft=fittype(fn,'independent','t');  t=[1800:10:2000]'; 
[f, st]=fit(t,x(2:end), ft, 'StartPoint',rand(1,2),...
    'Lower',[0,280],'Upper',[0.1,1000])  %由先验知识主观确定参数界限
xs=coeffvalues(f)  %显示拟合的参数值
xh=f(2010)  %求2010年的预测值
 
a=[ones(21,1), -x(1:end-1)];  %向前差分
b=diff(x)./x(1:end-1)/10;
cs=a\b; r=cs(1), xm=r/cs(2)
xh=fn(r,xm,2010)  %求2010年的预测值
 
a=[ones(21,1), -x(2:end)];  %向后差分
b=diff(x)./x(2:end)/10;
cs=a\b; r=cs(1), xm=r/cs(2)
xh=fn(r,xm,2010)  %求2010年的预测值

4. 两个种群的相互作用模型

基本假设
设在 时刻两个种群群体的数量分别为$x_1(t)$和$x_2(t)$。假定：
（1）初始时两个种群的数量均较小，分别为$x_1^0,x_2^0$，此后在一定时期内各自以自然增长率增长。
（2）每一种群群体的增长都受到Logistic规律的制约。设其自然增长率分别为$r_1$和$r_2$，在同一资源、环境制约下只维持第一个或第二个种群群体的生存极限数分别为$K_1$和$K_2$。
（3）两个种群依靠同一种资源生存，这两个种群的数量越多，可获得的资源就越少，从而物种的种群增长率就会降低。而且随着种群群体数量的增加，各自种群的增量都会对对方种群的变化产生一定的限制影响。

模型建立
对第一个种群而言，若第二个种群的每个个体消耗资源量相当于第一个种群每个个体消耗资源量的${\alpha }_1$倍，则第一个种群群体数量的增长率为
$r_1\left(1-\frac{x_1+{\alpha }_1{x}_2}{K_1}\right)x_1.$
同理，设第一个种群的每个个体消耗的资源量为第二个种群每个个体消耗资源量的${\alpha }_2$倍，则第二个种群群体数量的增长率为
$r_2\left(1-\frac{x_2+{\alpha }_2{x}_1}{K_2}\right)x_2.$

综上所述，两个种群竞争系统的群体总数$x_1(t)$和$x_2(t)$应满足微分方程组：
$$\{\begin{array}{l}\frac{d{x}_1}{dt}=r1(1-\frac{x_1+\alpha 1x2}{K1})x_1,\\ \frac{d{x}_2}{dt}=r2(1-\frac{x_2+\alpha 2x1}{K2})x_2,\\ x_1(0)=x_1^0,x_2(0)=x_2^0.\end{array}(6.37)$$

弱肉强食模型
下面讨论另一类生物链问题的数学建模问题，即弱肉强食问题。此类问题广泛存在于自然界中，如大鱼吃小鱼、狼群与羊群等。
设$t$时刻第一个种群的数量和第二个种群的数量分别为$x_1(t)$和$x_2(t)$。初始时种群数量分别为$x_1^0,x_2^0$。

基本假设
（1）第一个种群的生物捕食第二个种群的生物，其种群数量的变化除了自身受Logistic规律的制约外，还受到被捕食的第二个种群的数量影响。
（2）第二个种群的数量变化除了自身受自限规律影响外，还受其天敌第一个种群的数量影响。第二个种群的种群数量越多，被捕杀的机会越多，从而第一个种群的繁殖越快。
（3）设两个种群的自然增长率分别为$r_1$和$r_2$，各自独自生存的生存极限数分别为$K_1$和$A{K}_2$。

模型建立
对强食型的第一个种群，其种群数量的增长率受自身Logistic规律限制，同时还受第二个种群供应水平影响：供应能力（即第二个种群数量）越强，对第一个种群的数量增长刺激越明显。不妨假定单位时间内第一个种群的单位个体与第二个种群的有效接触数为$b_{12}{x}_2(t)$，其中$b_{12}>0$为比例系数。
$$\frac{d{x}_1}{dt}=r_1\left(1-\frac{x_1}{K_1}\right)x_1+b_{12}{x}_1{x}_2.$$

另一方面，第一个种群数量越多，第二个种群被捕杀的数量也就越多，从而种群数量减少得越快，考虑到自限规律的因素，第二个种群的增长率应为
$$\frac{d{x}_2}{dt}=r_2\left(1-\frac{x_2}{K_2}\right)x_2-b_{21}{x}_1{x}_2，$$
其中，$b_{21}$为正实数，为两个种群之间的接触系数。

因而$x_1(t)$,$x_2(t)$应满足微分方程组：
$$\{\begin{array}{l}\frac{d{x}_1}{dt}=r_1(1-\frac{x_1}{K_1})x_1+b_{12}{x}_1{x}_2\\ \frac{d{x}_2}{dt}=r_2(1-\frac{x_2}{K_2})x_2-b_{21}{x}_1{x}_2\\ x_1(0)=x_1^0,x_2(0)=x_2^0.\end{array}$$

该模型是较经典的捕食者-被捕食者模型的一种形式。更多的讨论可基于食物链中各种群的增长规律以及种群之间的相互依存关系，建立基于各类具体问题的数学模型。
把类似的讨论应用到其他研究领域，也可以得到相似的模型，如市场中同类商品价格的相互竞争问题等。
在不同的假设下，可以得到不同的捕食者-被捕食者模型。

【例19】 捕食者-被捕食者方程组
$$\{\begin{array}{l}\frac{dx}{dt}=ax-bxy,x(0)=60,\\ \frac{dy}{dt}=-cy+dxy,y(0)=30.\end{array}(6.38)$$
其中$x(t)$表示$t$个月之后兔子的总体数量，$y(t)$表示狐狸的总体数量，参数$a,b,c,d$未知。利用表6.2的13个观测值，拟合式(6.38)中的参数$a,b,c,d$。

【解】 微分方程对应的差分方程为
$$\{\begin{array}{l}(x_{k+1}-x_k)/(t_{k+1}-t_k)=a{x}_k-b{x}_k{y}_k,\\ (y_{k+1}-y_k)/(t_k+1-t_k)=-c{y}_k+d{x}_k{y}_k,k=0,1,\cdots ,11,\end{array}(6.39)$$
其中$x_k,y_k,t_k$分别为$x,y,t$的第k个观测值，$k=0,1,\cdots ,12$。

可以将(6.39)式改写成下列格式
$$[\begin{array}{cccc}{x}_k& -x_k{y}_k& 0& 0\\ 0& 0& -y_k& x_k{y}_k)\end{array}][\begin{array}{c}a\\ b\\ c\\ d\end{array}]=[\begin{array}{c}(x_{k+1}-x_k)/(t_{k+1}-t_k)\\ (y_{k+1}-y_k)/(t_{k+1}-t_k)\end{array}],k=0,1,\cdots ,11.(6.40)$$

上述所有的差分方程可以写成矩阵格式
$$\left[\begin{array}{cccc}{x}_0& -x_0{y}_0& 0& 0\\ ⋮& ⋮& ⋮& ⋮\\ x_{11}& -x_{11}{y}_{11}& 0& 0\\ 0& 0& -y_0& x_0{y}_0\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& -y_{11}& x_{11}{y}_{11}\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\\ d\end{array}\right]=\left[\begin{array}{c}(x_1-x_0)/(t_1-t_0)\\ ⋮\\ (x_{12}-x_{11})/(t_{12}-t_{11})\\ (y_1-y_0)/(t_1-t_0)\\ ⋮\\ (y_{12}-y_{11})/(t_{12}-t_{11})\end{array}\right].$$
利用线性最小二乘法，求得$a=0.1907，b=0.0048，c=0.4829，d=0.0095$。

clc, clear
t0=[0	1	2	3	4	5	6	8	10	12	14	16	18]';
x0=[60	63	64	63	61	58	53	44	39	38	41	46	53]';
y0=[30	34	38	44	50	55	58	56	47	38	30	27	26]';
dt=diff(t0); dx=diff(x0); dy=diff(y0);
temp=x0(1:end-1).*y0(1:end-1);
mat=[x0(1:end-1), -temp,zeros(12,2)
    zeros(12,2), -y0(1:end-1), temp]; %构造线性方程组的系数矩阵
const=[dx./dt; dy./dt];  %构造线性方程组的常数项列
abcd=mat\const  %拟合参数a,b,c,d

【例20】 捕食者-被捕食者方程组
捕食者-被捕食者方程组
$$\{\begin{array}{l}\frac{dx}{dt}=0.2x-0.005xy,x(0)=70,\\ \frac{dy}{dt}=-0.5y+0.01xy,y(0)=40.\end{array}(6.41)$$

其中$x(t)$表示$t$个月之后兔子的总体数量，$y(t)$表示狐狸的总体数量。
研究如下问题：
（1）$x(t)$和$y(t)$的总体数量变化的周期；
（2）$x(t)$的总体数量的最大值和最小值；
（3）$y(t)$的总体数量的最大值和最小值。

【解】
$$\{\begin{array}{l}0.2x-0.005xy=0,\\ -0.5y+0.01xy=0,\end{array}$$
得临界点$(50,40)$ ，它是一个稳定的中心，表示50只兔子和40只狐狸的平衡数量。方程组(6.41)的轨线和方向场见图6.10。方向场说明点$(x(t),y(t))$以逆时针方向沿其轨道运行，并且兔子和狐狸的总体数量分别在它们的最大值和最小值之间周期性振动。相位平面图像的一个不足之处是它没有给出每条轨道变化的速度。

通过做每个物种总体数量的（关于时间$t$的）函数的图像，可以弥补这个“缺少时间意义”的不足。在图6.11中，做出了$x(t)$和$y(t)$数值解的图像。通过两个相邻极小点的时间间隔就可以得到每个物种总体数量变化的周期$T$为20.2201个月，计算得$x(t)$的最大值为69.8810，x(t)的最小值为34.2513。$y(t)$的最大值为66.7664，$y(t)$的最小值为21.5644。

clc, clear, close all, L=100; %求解时间长度
x=20:5:80; y=10:5:80; [x,y]=meshgrid(x,y);
u=0.2*x-0.005*x.*y; v=-0.5*y+0.01*x.*y;
quiver(x,y,u,v), hold on
dxy=@(t,z)[0.2*z(1)-0.005*z(1)*z(2)
          -0.5*z(2)+0.01*z(1)*z(2)]; %定义微分方程组的右端项
sol=ode45(dxy,[0,100],[70;40]);
xt=@(t)deval(sol,t,1); yt=@(t)deval(sol,t,2);
fplot(xt,yt,[0,L]), xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex','Rotation',0)
figure(2), hold on
fplot(xt,[0,100],'--','LineWidth',1.5)  %画x(t)的解曲线
fplot(yt,[0,100],'LineWidth',1.5) %画y(t)的解曲线
xlabel('$t$','Interpreter','Latex')
legend({'$x(t)$','$y(t)$'},'Interpreter','latex')
fprintf('请在虚线上选取两个相邻的极小点！\n')
[t1,xx]=ginput(2) %用鼠标取x(t)上的相邻2个极小点
fprintf('请在实线上选取两个相邻的极小点！\n')
[t2,yy]=ginput(2) %用鼠标取y(t)上的相邻2个极小点
T1=diff(t1), T2=diff(t2)  %求近似的周期
[xt1,fx1]=fminbnd(xt,0,L)   %求x的最小点及最小值
[xt2,fx2]=fminbnd(@(t)-xt(t),0,L) %求x的最大点及最大值
[yt1,fy1]=fminbnd(yt,0,20)   %求y在[0,20]上的最小点及最小值
[yt2,fy2]=fminbnd(yt,20,40)   %求y在[20,40]上的最小点及最小值
T=yt2-yt1  %求精确的周期
[yt3,fy3]=fminbnd(@(t)-yt(t),0,L) %求y的最大点及最大值

5. 放射性废料的处理

【例20】
美国原子能委员会以往处理浓缩的放射性废料的方法，一直是把它们装入密封的圆桶里，然后扔到水深为90多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辩说这是不可能的。为此工程师们进行了碰撞实验。发现当圆桶下沉速度超过12.2 m/s 与海底相撞时，圆桶就可能发生碰裂。
这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少? 已知圆桶质量$m=239.46kg$，体积$V=0.2058m^3$，海水密度$\rho =1035.71kg/m^3$，若圆桶速度小于12.2 m/s就说明这种方法是安全可靠的，否则就要禁止使用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数$k=0.6$。现要求建立合理的数学模型，解决如下实际问题：

（1）判断这种处理废料的方法是否合理?
（2）一般情况下，$v$大，$k$也大；$v$小，$k$也小。当$v$很大时，常用$kv$来代替$k$，那么这时速度与时间关系如何? 并求出当速度不超过12.2 m/s，圆桶的运动时间$t$和位移$s$应不超过多少? ($k$的值仍设为0.6)

1. 问题一的模型

又因为$F=ma=m\frac{dv}{dt}=m\frac{d^{2}s}{d{t}^2}$，$G=mg$，$H=\rho gV$以及$f=kv=k\frac{ds}{dt}$，可得到圆桶的位移和速度分别满足下面的微分方程
$$m\frac{d^{2}s}{d{t}^2}=mg-gV-k\frac{ds}{dt}，(6.43)$$
$$m\frac{dv}{dt}=mg-\rho gV-kv.(6.44)$$

根据方程(6.43)，加上初始条件${\frac{ds}{dt}|}_{t=0}=s|t=0=0$，求得位移函数为
$$s(t)=-171511.0+429.7444t+171511.0e^{-0.002505638t}.(6.45)$$
由方程(6.44)，加上初始条件${v|}_{t=0}=0$，求得速度函数为
$$v(t)=429.7444-429.7444e^{-0.002505638t}.(6.46)$$

由$s(t)=90m$，求得圆筒到达水深$90m$的海底需要时间$t=12.9994s$，再把它带入(6.46)式，求出圆桶到达海底的速度为$v=13.7720\mathrm{m}/\mathrm{s}$。
显然此圆桶的速度已超过$12.2m/s$，可以得出这种处理废料的方法不合理。因此，美国原子能委员会已经禁止用这种方法来处理放射性废料。

2. 问题二的模型

由题设条件，圆桶受到的阻力应改为$f=k{v}^2=k(\frac{ds}{dt}{)}^2$， 类似问题一的模型，可得到圆桶的速度应满足如下的微分方程
KaTeX parse error: Undefined control sequence: \vthicksp at position 19: …frac{dv}{dt}=mg\̲v̲t̲h̲i̲c̲k̲s̲p̲-\rho gV-kv^2. …
根据方程(6.47)，加上初始条件KaTeX parse error: Undefined control sequence: \vthicksp at position 20: …ft|t=0\right.=0\̲v̲t̲h̲i̲c̲k̲s̲p̲，求出圆桶的速度
$v(t)=20.7303tanh(0.0519426t)$，
这时若速度要小于$12.2m/s$，那么经计算可得圆桶的运动时间就不能超过$T=13.0025s$，利用位移$s(T)={\int }_0^{T}v(t)dt$，计算得位移不能超过$84.8439m$。
通过这个模型，也可以得到原来处理核废料的方法是不合理的。

clc,clear
syms m V rho g k s(t) v(t) %定义符号常数和符号变量
ds=diff(s); %定义s的一阶导数，为了初值条件赋值
s1=dsolve(m*diff(s,2)-m*g+rho*g*V+k*diff(s),s(0)==0,ds(0)==0);
s1=subs(s1,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});  %常数赋值
ss1=vpa(s1,7)  %显示小数形式的位移函数
v1=dsolve(m*diff(v)-m*g+rho*g*V+k*v,v(0)==0);
v1=subs(v1,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
vv1=vpa(v1,7)  %显示小数形式的速度函数
t1=solve(s1-90); t1=double(t1)   %求到达海底90米处的时间
vvv1=subs(v1,t,t1); vvv1=double(vvv1)  %求到底海底90米处的速度
v2=dsolve(m*diff(v)-m*g+rho*g*V+k*v^2,v(0)==0);
v2=subs(v2,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
v2=simplify(v2); vv2=vpa(v2,6)  %显示小数形式的速度函数
t2=solve(v2-12.2); t2=double(t2)  %求时间的临界值
s2=int(v2,0,t2); s2=double(s2)  %求位移的临界值

3. 结果分析

由于在实际中$k$与$v$的关系很难确定， 所以上面的模型有它的局限性，而且对不同的介质，比如在水中与在空气中$k$与$v$的关系也不同。如果假设$k$为常数的话，那么水中的这个$k$就比在空气中对应的$k$要大一些。在一般情况下，$k$应是$v$的函数，即$k=k(v)$，至于是什么样的函数，这个问题至今还没有解决。

这个模型还可以推广到其它方面，比如说一个物体从高空落向地面的道理也是一样的。尽管物体越高，落到地面的速度越大，但决不会无限大。

---

总结

终于两天顶住了！