微分方程、动力系统与混沌导论 第1章 一阶方程[书摘]
微分方程模型在模型分析中的主要问题之一是稳定性分析。
微分方程模型的稳定性及其实际意义
用微分方程方法建立的动态模型问题中的一个重要问题是：当时间充分长后 ，动态过程的变化趋势是什么？微分方程模型中，方程 ( 组 ) + 初始条件 → 解。初始条件的作用在于确定解, 它的微小变化会产生不同的解，换言之，对解的发展性态变化，往往具有影响作用。问题是这种对解的发展性态的影响作用是长期存在的，还是当时间充分大以后 , 影响作用会 “消逝”?有时候，初始条件的微小变化会导致解的性态随时间变大后 , 产生显著的差异，这时称系统是不稳定的；有时候，初始条件变化导致解的性态差异会随时间变大后而消失，这时称该系统是稳定的。
在实际问题中，初始状态不能精确地而只能近似地确定，所以稳定性问题的研究对于用微分方程方法建立的模型具有十分重要的实际意义。也就是说，在具有稳定性特征的微分方程模型中，长远来看, 最终发展结果与精确的初始状态究竟如何，两者之间没有多大关系，初始状态刻画得精确不精确是无关紧要的。
微分方程稳定性理论可以使我们在很多情况下不求解方程便可直接得到微分方程模型描绘的系统是稳定或不稳定的结论。研究者对于微分方程稳定性理论的研究兴趣往往大于该方程解有无解析表达式的研究兴趣。
在数学建模竞赛活动中，很多问题中涉及到的微分方程是一类称为自治系统的方程。
自治方程：是指方程中不显含自变量$t$的微分方程，例如一阶方程
\[x'(t) = f(x)\]
二阶方程
\[\left\{ \begin{array}{l}x'(t) = f(x,y)\\y'(t) = g(x,y)\end{array} \right.\]
自治方程中的解随时间不断变大，如有稳定变化趋势，则这个解的最终趋势值只能是该方程的平衡点。一阶微分方程$x'(t) = f(x)$的平衡点是指代数方程$f(x)=0$的根(可能不止一个根) 。二阶微分方程组的平衡点是指代数方程组
\[\left\{ {\begin{array}{*{20}{c}}{f(x\;,\;y) = 0}\\{g(x\;,\;y) = 0}\end{array}} \right.\]
的解$x_0,y_0$(可能不止一对)。
1.1 最简单的例子
微分方程
\[\frac{{dx}}{{dt}} = ax,\]
方程告诉我们，对于\(t\)的任一取值，下面的关系式成立：
\[x'(t)=ax(t).\]
设\(k\)为任一给定的实数，则函数\(x(t)=ke^{at}\)就是一个解，并且这个方程没有其它解。我们把一个微分方程的所有解的全体称为这个方程的通解。
当解在一个点\(t_0\)处的取值\(u_0\)确定后，出现在解中的常数\(k\)就完全确定。从而这个方程在满足特定初始条件\(x(t_0)=u_0\)时有唯一解。
通常将上面的问题写成初值问题的形式：
\[x' = ax,\;\;\;x(0) = {u_0}.\]
注意这个微分方程在\(k=0\)时有一个特别的解，即常值解\(x(t) \equiv 0\)。像这样的常值解称为该方程的平衡解或平衡点。平衡解通常都是微分方程最重要的解。
在方程\(x'=ax\)中，常数\(a\)可以看作一个参数，当\(a\)变化时，方程就变了，当然解也要随之变化。我们能否定性地描述解的变化呢？这里\(a\)的符号将是关键：
(1)如果\(a>0\)，当\(k>0\)时，极限\({\lim _{t \to \infty }}k{e^{at}}\)等于\(\infty\)，当\(k<0\)时，该极限等于-\(\infty\)；
(2)如果\(a=0\)，\(ke^{at}\)是常数；
(3)如果\(a<0\)，\({\lim _{t \to \infty }}k{e^{at}}=0\)。
通过作解的简图可以形象地描述解的定性行为。当平衡点附近的解都远离它时，我们称平衡点是源点；而当附近的解都趋于它时，该平衡点称为一个汇点。
我们还可以将解画在相线(phase line)上来描述解。因为解\(x(t)\)为时间的函数，我们可以将\(x(t)\)看成一个沿实直线运动的质点。在平衡点处，质点保持不动(用一个实心圆点表示，\(a=0\)时算什么点呢？)，而其它解则沿着\(x\)轴上下运动，在图1.1中用箭头表示。
在某种意义下，当\(a \ne 0\)时，方程\(x'=ax\)是稳定的。精确地说，当\(a\)用一个与之同号的\(b\)替换时，解的定性行为不发生改变。但当\(a=0\)时，\(a\)的微小改变都将根本地改变解的行为。于是，我们说单参数微分方程的族\(x'=ax\)在\(a=0\)处(而不是针对任意常数\(k\)而言的)出现了一个分岔。
1.2 Logistic物种模型
前面的微分方程\(x'=ax(a>0)\)可以看成一个简单的物群增长模型。实际的问题不可能无限制地增长。为了在物种总量模型中将这种限制考虑进来，我们进一步假设：
(1)当物种较少时，增长率几乎与总量成正比；
(2)但当物种增长到很多时，增长率就变成负的。
Logistic物种增长模型就是满足这些假设的一个微分方程，其方程为：
\[x'=ax(1-\frac{x}{N}).\]
不失一般性，我们假设\(N=1\)，也即选取单位使得承载量正好是1单位的总量，而\(x(t)\)则代表在\(t\)时刻的总量占理想总量的比例。这样合理方程就简化成
\[x'=f_a(x)=ax(1-x).\]
这是一个一阶、自治、非线性微分方程的例子。其通解为
\[x(t)=\frac{Ke^{at}}{1+Ke^{at}},\]
其中\(K\)是积分时产生的任意常数。将\(t=0\)代入上式可解得
\[K=\frac{x(0)}{1-x(0)}.\]
利用这个式子，可以将解重写为
\[\frac{x(0)e^{at}}{1-x(0)+x(0)e^{at}}.\]
此方程有两个平衡解\(x(t)\equiv 0\)和\(x(t)\equiv 1\)。
为了对解的行为有一些定性认识，我们画出该方程的斜率场。从图上易见，所有对应于初值\(x(0)>0\)的解都走向于理想总量\(x(t)\equiv 1\)，这正好与我们的假设吻合。当初值\(x(0)<0\)时，解将趋于-\(\infty\)。当然这些解在物种模型中是没有意义的。
注意，我们还可以从函数\(f_a(x)=ax(1-x)\)的图像上认识这些行为。该函数图像(见图1.4)与\(x\)轴交于\(x=0\)和\(x=1\)两个点，这正对应于两个平衡点。当\(0<x<1\)时，\(f_a(x)>0\)，从而在任何满足\(0<x<1\)时的(\(t,x)\)点处，斜率为正数，从而解在这个区域将增加。而在\(x<0\)或\(x>1\)时，\(f_a(x)<0\)，故解将减小，正如图1.3的解图像和相线所示。
同样，我们还可以从\(f_a\)的图像上看出：\(x=0\)是一源点而\(x=1\)是一汇点。在\(x=0\)附近，当\(x>0\)时，\(f_a(x)>0\)，斜率(即\(x'\))为正，解增加；但当\(x<0\)时，\(f_a(x)<0\)，斜率为负，解减小。这意味着附近的解都要远离0，即0为源点。同样，1为汇点。
此外，我们还可以通过分析得到这些信息。由于\(f'_a(x)=a-2x\)，从而有\(f'_a(0)=a>0\)而\(f'_a(1)=-a<0\)。由于\(f'_a(0)>0\)，当\(x\)通过0时，斜率将单调增加(其实就是说图1.4中，\(f_a(x)\)通过0时，是单调上升的)，于是在\(x=0\)的下方(对应图1.3解的图像，对应图1.4，则是\(x=0\)的左边)，斜率值取负，而在\(x=0\)的上方，斜率值取正值。因此，解要远离\(x=0\)。同样，\(f'_a(1)<0\)将使得解趋于\(x=1\)，从而使得这个平衡点成为一汇点。在今后，我们会遇到很多这样通过计算导数来确定平衡点附近解的定性行为的例子。
[总结：先令\(x'=f(x)=0\)求出平衡点，然后根据平衡点处的导数\(f'(x)\)的符号来判定平衡点的类型，即源点或汇点，若\(f'(x)>0\)，为源点；若\(f'(x)<0\)，为汇点]
1.3 常值收割与分岔
现在将物种的收割考虑进来以修改Logistic模型。假设物种遵循参数\(a=1\)时的Logistic假设，但它们同时以常速率\(h\)被收割。此时，微分方程为
\[x' = x(1 − x) − h\]
其中\(h \ge 0\)为一新参数。
下面我们不解方程而直接利用函数
\[f_h(x) = x(1 − x) − h\]
的图像来“读出”解的定性行为。在\(0 < h < 1/4, h = 1/4,\)及\(h > 1/4\)三种不同的情况下，图1.6作出了\(f_h\)的图像。当\(h\)通过\(h=1/4\)时，另一种分岔现象发生了：当\(h\)单调增加通过1/4时，两个平衡点\(x_l\)和\(x_r\)重合，而当\(h > 1/4\)时，平衡点消失。事实上，当\(h > 1/4\)时，对所有的\(x\)，总有\(f_h(x)<0\)。数学上，这意味着随着时间的增长，微分方程的解都将递减到\(-\infty\)。
我们用分岔图来形象地记录这些变化。在该图中，我们用\(h\)代表横坐标，而对应于每一个\(h\)值，画出相应的相线。图中的曲线代表每个\(h\)值对应的平衡点。这让我们从另一角度看到，当\(h\)通过1/4时，源点和汇点融合成一个平衡点，然后消失。
在生态学上，这种分岔对应于所研究的物种的灾难。当收割率为1/4时或更低时，只要初始总量充分大(\(x(0) \ge x_l)\)，总量就能保持。当\(h=1/4\)时，收割率的微小变化都将导致物种命运的大变化，例如，只要收割率\(h>1/4\)，该物种就要灭绝。
例 我们来看另一种分岔现象。为此，我们考虑微分方程族
\[x' = g_a(x) = x^2 − ax = x(x − a),\]
对应的分岔图如图1.8所示。
1.4 周期收割与周期解
现在不再假设在Logistic模型的收割率为常数。例如，许多种类的鱼在温暖的季节的收割率比在寒冷月份的收割率要高。于是，假设物种的收割率是周期变化的。
\[x' = f (t , x) = ax(1 − x) − h(1 + sin(2\pi t ))\]
就是这样的一个模型，其中\(a\)和\(h\)都是正参数。该方程显式地依赖于时间，从而这是一个非自治的微分方程的例子。该方程不再是可分离变量的，我们无法用通常微积分的方法求出解的解析表达式，这迫使我们去采用更加定性化的方法。
为了描述这种情况下总量的变化趋势，首先注意到，该微分方程的右端对时间变量是以1为周期的，即\(f(t+1,x)=f(t,x)\)。这一事实多少可以简化求解问题。假设我们知道了所有初值问题在\(0 \le t \le 1\)时(不是所有时间!)的解，则事实上我们就知道了任意时刻的解。假设\(x_1(t)\)是定义在\(0 \le t \le 1\)上的一个解，满足\(x_1(0)=x_0\)，而\(x_2(t)\)为满足\(x_2(0)=x_1(1)\)的一个解，则可以通过定义\(x_1(t + 1) = x_2(t )(0 \le t \le 1)\)将解\(x_1(t)\)进行延伸。由于
\[x'_1(t + 1) = x'_2(t ) = f (t , x_2(t ))= f (t + 1, x_1(t + 1)).\]
因而延伸后的函数仍然是一个解。从而，一旦知道了解在区间\(0 \le t \le 1\)上的行为，我们就可以用同样的方法得到解在所有时间区间的行为(可理解为解的拼接)。
其次，假设知道了满足任一初始条件\(x(0)= x_0\)的解在在时刻\(t=1\)时的取值\(x(1)\)，则对于每一个初值\(x_0\)，都对应于满足\(x(0)= x_0\)的解\(x(t)\)的一个值\(x(1)\)。这样就定义了一个函数\(p(x_0) = x(1)\)(相当于给定一个起点，确定一个终点)。将这个函数与自身作一次复合，就得到以\(x_0\)为初值的解在\(t=2\)处的取值，即\(p(p(x_0)) =x(2)\)。将这个函数与自身作\(n\)次复合，就可以算出该解在时刻\(n\)的值，从而也就知道了该解的变化趋势。
上述定义的函数\(p\)称为该微分方程的庞加莱映射。有了这样一个函数，我们就可以从连续动力系统(微分方程)的领域转化到较易理解的离散运动系统的领域(迭代函数)。例如，假设对某一个初值\(x_0\)有\(p(x_0)=x_0\)，即\(x_0\)为函数\(p\)的不动点。则从前面的观察可知，对每一整数\(n\)，都有\(x(n)=x_0\)，进一步，对于满足\(0<t<1\)的任何\(t\)，也有\(x(t ) = x(t + 1)\)，从而对于所有的整数\(n\)都有\(x(t + n) = x(t )\)。这说明满足初始条件\(x(0) = x_0\)的解关于\(t\)是一个以1为周期的周期函数。这样的解称为微分方程的周期解。在图1.10中，我们画出了具有周期收割的Logistic模型的几个周期解(最上面的曲线和最下面的曲线)。
不幸的是，计算微分方程的庞加莱映射通常是一件不可能的事。所幸对于具有周期收割的Logistic方程，我们可以做到。
1.5 计算庞加莱映射
在计算该方程的庞加莱映射之前，我们先引入一些重要的术语。为了强调解对初值$x_0$的依赖关系，将相应的解记为$\phi(t,x_0)$。函数\(\phi :\mathbb{R} \times \mathbb{R} \to \mathbb{R}\)称为该微分方程诱导的流。如果将变量$x_0$固定，则函数
\[t \to \phi(t,x_0)\]
正是该微分方程对应于初值$x_0$的解的另一种表述方式。有时我们也将这个函数记为$\phi_t(x_0)$。
现在我们回到带周期收割的Logistic方程
\[x' = f (t , x) = ax(1 − x) − h(1 + sin(2\pi t ))\]
该方程满足初值条件$x(0)=x_0)$的解由$t \to \phi(t,x_0)$给出。虽然我们不知道具体的表达式，但根据微积分基本定理[附注1]，以及
\[\frac{{\partial \phi }}{{\partial t}}\left( {t,{x_0}} \right) = f\left( {t,\phi \left( {t,{x_0}} \right)} \right),\;\;\;\;\phi \left( {0,{x_{\rm{0}}}} \right){\rm{ = }}{x_0},\]
可知这个解满足
\[\phi(t,x_0) = x_0 + \int_0^tf(s, \phi(s,x_0))ds. \]
将上式对$x_0$求导，由链式法则可得
\[\frac{\partial \phi }{\partial {{x}_{0}}}\left( t,{{x}_{0}} \right)=1+\int_{0}^{t}{\frac{\partial f}{\partial \phi }\left( s,\phi \left( s,{{x}_{0}} \right) \right)\centerdot \frac{\partial \phi }{\partial {{x}_{0}}}\left( s,{{x}_{0}} \right)ds}.\]
现在令
\[z(t)=\frac{\partial \phi}{\partial x_0}\left( t,x_0 \right).\]
注意到
\[z(0) = \frac {\partial \phi}{\partial x_0} (0,x_0) = 1, \]
将$z$对$t$求导可得(此处用到了微积分基本定理：函数$f(x)$积分的导数为$f(x)$)
\[{z}'(t)=\frac{\partial f}{\partial \phi }\left( t,\phi (t,{{x}_{0}}) \right)\centerdot \frac{\partial \phi }{\partial {{x}_{0}}}(t,{{x}_{0}})=\frac{\partial f}{\partial \phi }\left( t,\phi (t,{{x}_{0}}) \right)\centerdot z(t).\]
虽然还不能显式地知道$\phi(t,x_0)$，但上一方程却告诉我们$z(t)$满足微分方程
\[{z}'(t)=\frac{\partial f}{\partial \phi }\left( t,\phi (t,{{x}_{0}}) \right)\centerdot z(t).\]
以及$z(0)=1$。于是，通过分离变量，可求得该方程的解为
\[z(t)=\exp \int_0^t \frac{\partial f}{\partial \phi} \left( s,\phi \left( s,x_0 \right) \right) ds,\]
从而可得(即令上式中的$t=1$)
\[\frac{\partial \phi}{\partial x_0}(1,x_0)= \exp \int_0^1 \frac{\partial f}{\partial \phi} \left( s,\phi \left( s,x_0 \right) \right) ds.\]
由于$p(x_0)=\phi(1,x_0)$，这样我们就确定了庞加莱映射的导数$p'(x_0)$。注意到$p'(x_0)>0$(因为是指数函数)，故$p$为单调增加函数。
再次微分可得
\[{p}''({{x}_{0}})={p}'({{x}_{0}})\left( \int_{0}^{1}{\frac{{{\partial }^{2}}f}{\partial \phi \partial \phi }(s,\phi (s,{{x}_{0}}))\centerdot \exp\left( \int_{0}^{s}{\frac{\partial f}{\partial \phi }(u,\phi (u,{{x}_{0}}))du} \right)ds} \right),\]
这个式子看起来可怕，但是由于
\[f(t,x_0) = ax_0(1 - x_0) - h(1+\sin (2\pi t)),\]
因而
\[\frac{\partial^2f}{\partial \phi \partial \phi} \equiv -2a,\]
故可得$p''(x_0)<0$。于是，庞加莱映射的图像是向下凹的。这蕴涵$p$的图像至多穿过对角线$y=x$两次，即$x$至多只有两个值满足$p(x)=x$。从而庞加莱映射至多只有两个不动点。这些不动点对应于原微分方程的周期解，而且这些解满足，对所有的$t,x(t+1) = x(t)$。我们还可以这样说，当初值$x_0$为这些不动点之一时，流$\phi(t,x_0)$关于$t$是以1为周期的周期解。在$h=0.8$的特例时，在图1.10中，我们看见了两个这样的解。在图1.11中，又有两个解看起来是周期的。注意其中一个周期解看起来吸引附近所有的解，而另一个排斥附近所有的解。后续还会经常讨论这些概念，并使得它们更精确。
回想一下，该微分方程还依赖于收割参数$h$。如图1.11所示，当$h$比较小时，有两个不动点。将$f$对$h$微分可得
\[\frac{\partial f}{\partial h}(t,x_0) = -(1+\sin 2\pi t),\]
于是，$\partial f/ \partial h<0$(除去$t=3/4$的情形)。这蕴涵在任一点$(t,x_0)$处的斜率随着$h$的增加而减小。从而，庞加莱映射的值也随着$h$的增加而减小。从而$h$有唯一的取值$h_*$，使得在该点处庞加莱映射恰好只有一个不动点。而当$h>h_*$时，$p$没有不动点，此时对所有的初值都有$p(x_0)<x_0$。这意味着该物种将灭绝。
[附注1]：微积分基本定理
微积分基本定理描述了微积分的两个主要运算──微分和积分之间的关系。
定理的第一部分，有时称为微积分第一基本定理，表明不定积分是微分的逆运算。这一部分定理的重要之处在于它保证了某连续函数的原函数的存在性。
定理的第二部分，有时称为微积分第二基本定理，表明定积分可以用无穷多个原函数的任意一个来计算。这一部分有很多实际应用，这是因为它大大简化了定积分的计算。
微积分基本定理（FTC）有两个部分，第一部分是关于原函数的导数，第二部分描述了原函数和定积分之间的关系。
第一部分 第一基本定理
设$f$为定义在闭区间$[a,b]$的实函数
$F(x)=\int_a^xf(t)dt,x\in [a,b]$
那么，$F(x)$可导，且$F'(x)=f(x)$。
第二部分 第二基本定理
设$f$为定义在闭区间$[a,b]$的连续实数函数。设$F$为$f$的一个原函数，也就是说，它是使下式成立的无穷多个函数之一，
- $f(x)=F'(x).$
那么
$\int_a^bf(x)dx=F(b)-F(a).$
