6.4 线性系统的稳定性与相图

在前面几节里面, 我们学习了如何去求解一个线性微分方程组. 但值得注意的是能够求解的微分方程仅为少数. 很多非线性微分方程组有着非常复杂的结构, 也因此没有通解. 比如Lorenz方程:

$$\begin{pmatrix} x'(t) \\ y'(t) \\ z'(t) \end{pmatrix} = \begin{pmatrix} \sigma(y-x) \\ x(\rho-z) - y \\ xy-\beta z\end{pmatrix},\quad \sigma,\rho,\beta\in\mathbb{R}.$$

在这个系统中, 微小的参数变化便可能会对方程的解产生很大的影响. 比如在下图中, 我们取$\sigma=10,\beta=\frac{8}{3}$, 以及相同的初始条件$x(0)=y(0)=z(0)=1$. 我们随后选取不同的参数$\rho$, 得到下图中$x(t),z(t)$的变化关系:

在这种情况下, 我们去尝试求解便不再现实了. 对于 Lorenz 方程这样的非线性系统,我们通常无法像常系数线性系统那样利用特征值, 特征向量和矩阵指数写出通解. 然而, 我们仍然可以从几何角度研究解的行为: 例如,解是否会趋近某个点? 是否会远离某个点?不同初始条件下的轨迹是否会呈现相似的形状, 等等. 如果我们回到上图, 我们发现当$\rho=10$时, 随着时间增大, 解也在逐渐趋近于某一个点; 当$\rho$取值变大后我们便发现解就变得更加无序. 研究这些问题并不要求我们知道解的精确公式, 而是关注解在空间中的整体运动趋势, 这便引出了稳定性与相图的研究.

为了更好地理解稳定性的概念, 我们先从一维的微分方程$x'(t)=f(x)$出发.

当$f(x)>0$时, 我们有$x'(t)>0$, 即$x(t)$单调递增, 解向右运动. 我们不妨考虑$x'(t)=1$, 在这里面我们知道$f(x)=1>0$, 因此我们在不求解方程的情况下便知道解向右运动. 如果对原方程进行求解我们有$x(t)=t+C$, 那么不难发现$x(t)$随着$t$的增大而增大, 因此在$x(t)$坐标轴上随着$t$的增大$x(t)$向右运动. 同理若$f(x)<0$, $x(t)$单调递减, 解向左运动. 如果我们再考虑$x'(t)=-1$, 那么此时我们知道解向左运动. 求解得$x(t)=-t+C$. 因此随着时间$t$的增大$x(t)$在减小, 那么在$x(t)$坐标轴上随着$t$的增大$x(t)$向左运动.

我们如果考虑一个稍微复杂一点的情形: $x'(t) = x^2-1$. 那么我们知道$x=\pm 1$为该方程的平衡点. 当$x \in (-\infty,-1)$时, $x'(t)>0$, 因此在这个区间上解向右运动; 当$x \in (-1,1)$时, $x'(t)<0$, 因此解向左运动; 当$x \in (1,\infty)$时, $x'(t)>0$, 因此解向右运动. 如果我们在坐标轴上用箭头画出不同位置下$x$的运动方向, 我们得到如下的图:

我们发现, 在$x=-1$附近, 解会逐渐趋近于该点. 我们称这样的点为稳定平衡点; 而在$x=1$附近, 解会逐渐远离该点. 我们称这样的点为不稳定平衡点. 我们如果同$f(x)=x^2-1$的图像进行对比:

不难发现, 稳定平衡点对应了其导函数$f'(x)$的极大值点; 不稳定平衡点对应了其导函数$f'(x)$的极小值点. 对于稳定性我们进而给出更加严谨的定义:

Lyapunov稳定描述的是“初始值靠近则轨迹始终靠近”; 而准渐进稳定描述的是“初始值靠近则轨迹最终收敛”.若二者同时成立，则称平衡点渐近稳定. 在上面的例子$x'(t) = x^2-1$中, 点$x=-1$即为渐近稳定.

随后, 我们开始讨论线性系统的稳定性与相图. 为方便读者理解和简化概念, 我们在本节将重点学习$2\times 2$矩阵$\vec A$. 我们首先假设$\vec A \in M_2(\mathbb{R})$可逆, 这样一来 系统$\vec y'(t) = \vec A\vec y(t)$的平衡点便是$\vec A\vec y(t) = \vec 0$的解, 即$\vec y(t)=\vec 0$为唯一的平衡点. 在此条件下, 矩阵$\vec A$涵盖了以下几种情况:

(i) 矩阵$\vec A$可对角化, 且特征值$\lambda_1,\lambda_2$均为实数;

(ii) 矩阵$\vec A$有一对互为共轭复数的特征值;

(iii) 矩阵$\vec A$不可对角化.

针对不同的情况, 平衡点$\vec y(t)=\vec 0$有着不同的相图类型, 我们将平衡点附近轨迹的形状称为$\vec y(t) = \vec 0$的节点. 我们下面将逐个分析各个情况下节点的形状:

• 情况一: 矩阵$\vec A$可对角化, 且特征值为实数: 记$\lambda_1,\lambda_2$为特征值，$\vec u_1,\vec u_2$为与之对应的特征向量. 那么根据我们之前的推论, 微分方程$\vec y'(t) = \vec A\vec y(t)$的通解为

$$\vec y(t) = c_1 e^{\lambda_1 t}\vec u_1+c_2e^{\lambda_2 t}\vec u_2, \quad c_1,c_2\in\mathbb{R}.$$

在此情况下, 我们可以通过特征值的符号来判断$\vec y(t)=\vec 0$这个平衡点的稳定性.

• 若$\lambda_1,\lambda_2<0, \lambda_1 \neq \lambda_2$: 那么$e^{\lambda t} \to 0$. 因此

$$\lim_{t \to \infty} \vec y(t) = \lim_{t \to \infty} c_1 e^{\lambda_1 t}\vec u_1 + c_2 e^{\lambda_2 t}\vec u_2 = \vec 0,$$

故$\vec y(t)=\vec 0$为准渐进稳定. 同时结合$e^{\lambda t}$的单调性, 我们得出$\vec y(t)=\vec 0$为Lyapunov稳定. 因此当$\vec A$的所有特征值均为负数时, $\vec y(t)=\vec 0$为渐近稳定节点. 此时平衡点附近的轨迹有向原点移动的趋势.

• 若$\lambda_1,\lambda_2 >0, \lambda_1 \neq \lambda_2$: 那么不难发现

$$\lim_{t \to \infty} \vec y(t)$$

发散. 因此当$\vec A$的所有特征值均为正数时, $\vec y(t)=\vec 0$不稳定节点, 此时平衡点附近的轨迹方向有远离原点的趋势.

• 若$\vec A$有两个相同的特征值$\lambda_1=\lambda_2$: 此时$\vec y(t) = \vec 0$被称作是一个真节点 (Proper Node). 此时其稳定性取决于特征值的符号. 若特征值为正数, 则$\vec y(t)=\vec 0$不稳定; 若特征值为负数, 则$\vec y(t)=\vec 0$渐近稳定. 其对应的节点分别为渐近稳定节点和不稳定节点.

• 若两个特征值异号($\lambda_1\lambda_2<0$): 此时$\vec y(t)=\vec 0$被称作是鞍点 (Saddle Point). 此时沿着负特征值对应的特征向量方向的解会趋向原点; 沿着正特征值对应的特征向量方向的解会远离原点. 此时的平衡点不稳定.

回顾我们先前提到的方程$x'(t)=f(x(t))$, 其状态空间是一条数轴. 由于 $x'(t)>0$ 表示解向右运动; 而 $x'(t)<0$ 表示解向左运动. 我们便可以在数轴上画出向图,从而判断平衡点附近解的运动趋势. 那么同理, 对于二维系统

$$\begin{pmatrix} y_1'(t)\\ y_2'(t) \end{pmatrix} = \begin{pmatrix} f(y_1(t),y_2(t))\\ g(y_1(t),y_2(t)) \end{pmatrix},$$

系统的状态由点$(y_1(t),y_2(t))^\top\in\mathbb R^2$表示. 因此,我们把 $\mathbb R^2$ 称为该系统的相平面. 对于相平面中的每一点 $(y_1(t),y_2(t))^\top$, 我们都可以得到一个速度向量$(y_1'(t), y_2'(t))^\top$. 这个向量表示解经过点 $(y_1(t),y_2(t))$ 时的瞬时运动方向. 由这些向量组成的图像称为方向场, 而由不同初始条件得到的解轨迹所组成的图像称为相图. 因此我们可以类比在一维平面里$x'(t)=x^2-1$的相图, 去绘制$\vec y'(t) = \vec A\vec y(t)$的相图.

下面的几个相图便分别对应了我们上面讨论过的$4$种情况: 在下面的四幅图中, 蓝色实线表示特征向量的方向, 红色曲线表示点的运动轨迹, 粉色箭头即为速度方向向量.

• 情况二: 矩阵$\vec A$含有一对共轭特征值: 此时我们假设$\lambda_1=a+bi, \lambda_2=a-bi, b\neq 0$为$\vec A$的特征向量, $\vec u + i \vec v, \vec u - i\vec v$分别为$\lambda_1,\lambda_2$对应的特征向量. 那么方程$\vec y'(t)=\vec A\vec y$的解为

$$\vec y'(t) = c_1e^{at}\cdot\Big( \cos(bt)\vec u-\sin(bt)\vec v\Big) + c_2 e^{at} \cdot\Big( \sin(bt)\vec u+\cos(bt)\vec v\Big).$$

如果我们令$\vec p=c_1\vec u+c_2\vec v, \vec q = c_2\vec u - c_1\vec v$, 那么我们就有

$$\vec y'(t) = e^{at}\Big(\vec p \cos(bt) + \vec q \sin(bt) \Big).$$

不难发现, $\vec p \cos(bt) + \vec q \sin(bt) = \begin{pmatrix} p_1\cos(bt)+q_1\sin(bt) \\ p_2\cos(bt) + q_2\sin(bt) \end{pmatrix}$表示了一个经过旋转变换之后的椭圆:

$$\begin{pmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{pmatrix}\begin{pmatrix} a\cos\theta \\ b\sin\theta \end{pmatrix} = \begin{pmatrix} a\cos\phi \cos\theta - b\sin\phi\sin\theta\\ a\sin\phi\cos\theta + b\cos\phi\sin\theta \end{pmatrix},$$

然后通过选取合适的角度及参数, 我们便可以将其表示为$\vec p\cos(bt)+\vec q\sin(bt)$的形式. 至于$e^{at}$, 我们可以看成是一个控制椭圆放大或缩小的参数. 当$a>0$时, $e^{at}$的模随着时间增大而增大, 也就是说其离原点越来越远, 那么$\vec y(t)=\vec 0$便是一个不稳定点; 若$a<0$, 那么$e^{at}$的模随着时间的增大而逐渐趋于$0$, 也就是说其离远点越来越近, $\vec y(t) = \vec 0$便是一个渐近稳定点. 最后当$a=0$的时候我们发现其解集便是一个个相似的椭圆. 因此我们将目前的发现总结成以下三点:

下面我们分别设$\vec A_1=\begin{pmatrix} 6 & 15 \\ -3 & -6 \end{pmatrix}, \vec A_2 = \begin{pmatrix} -1 & 6 \\ -\frac{3}{2} & -1 \end{pmatrix}, \vec A_3 = \begin{pmatrix} 7 & -10 \\ 5 & -3 \end{pmatrix}$. 它们分别对应了矩阵$\vec A$的复数特征值的实部分别为$0,$负数, 正数的情况. 下面展示了它们的相图.

• 情况三: 矩阵$\vec A$不可对角化: 在前面的几种情形中, 我们都可以找到足够多的特征向量, 从而利用特征向量方向来描述相图. 然而, 当矩阵$\vec A$有重复特征值但缺少足够多的特征向量时,矩阵无法对角化.此时我们需要使用广义特征向量,而相图中也只会出现唯一的特征向量方向.这类平衡点附近的相图类型被称为退化节点 (Improper Node).在$2\times 2$矩阵$\vec A$中, $\vec A$不可对角化的条件是存在代数重数为$2$, 但几何重数为$1$的特征值$\lambda$. 我们记其广义特征向量为$\vec u \in K_1,\vec v \in K_2 \setminus K_1$. 利用广义特征向量的知识, 我们知道此时微分方程组$\vec y'(t) = \vec A\vec y(t)$的解满足

$$\vec y(t) = c_1 e^{\lambda t}\vec u+c_2 e^{\lambda t}(\vec v + t\vec u),\quad c_1,c_2\in\mathbb{R}.$$

当$\lambda<0$且$t \to \infty$ 时, 虽然$\vec y(t)$趋近于$\vec 0$, 但是包含$te^{\lambda t}\vec u$一项会影响轨线接近原点的方向. 大多数轨线在靠近原点时会变得越来越接近唯一的特征向量方向, 即$\vec u$的方向. 所以在相图上我们会看到很多轨线弯曲着靠近原点,并且最终都几乎沿着同一个方向进入原点:

下面的两幅图展示了此时系统$\vec y'(t) = \vec A\vec y(t)$在原点附近的相图:

此时读者不妨回忆第三章第$4$节的课后习题第一题: 现在如果再回去看这一道题, 是不是就小菜一碟了呢?

我们把本节里讨论过的平衡点$\vec y(t)=\vec 0$的稳定性以及其对应的相图类型总结成了以下的表格, 方便读者快速记忆:

{6.4 练习}

1. [分岔论初步] 回顾我们最开始提到的微分方程$x'(t) = x^2-1$, 我们绘制出了如下的相图:

现在, 我们引入一个变量$\mu$, 然后研究微分方程$x'(t) = f(x,\mu)= x^2-\mu, \mu\geq0$, 那么不难发现$x = \pm \sqrt{\mu}$即为两个平衡点. 此时我们想要知道:$\mu$的取值会如何影响平衡点的位置及性质呢? 因此在$(\mu,x)$平面中平衡点便是曲线$x=\sqrt{\mu}$与曲线$x=-\sqrt{\mu}$. 根据我们前面的讨论, $x=\sqrt{\mu}$为不稳定平衡点, 因此我们将曲线$x=\sqrt{\mu}$绘制成虚线; $x=-\sqrt{\mu}$为渐近稳定的平衡点, 我们将曲线绘制成实线. 如此一来, 我们便得到了一个分岔图 (Bifurcation Diagram):

我们称其为鞍结分岔 (Saddle Node Bifurcation). 其中粉色的竖直箭头即可表示当$\mu$为确定值时系统解的运动方向, $(x,\mu)=(0,0)$的位置便是一个分岔点.

(i) [跨临界分岔 (Transcritical Bifurcation)] 类比我们上面的推理, 现在我们设$x'(t) = f(x,\mu)=x^2 - \mu x$, $\mu \in \mathbb{R}$. 先求出这个微分方程的平衡点, 讨论其稳定性, 然后尝试画出这个系统的分岔图.

(ii) [Pitchfork分岔 (Pitchfork Bifurcation)] 类比我们上面的推理, 现在我们设$x'(t) = f(x,\mu)=x^3 - \mu x$, $\mu \in \mathbb{R}$. 先求出这个微分方程的平衡点, 讨论其稳定性, 然后尝试画出这个系统的分岔图.

在一维系统$x'(t)=f(x,\mu)$中,平衡点$x^*$的稳定性由$f_x(x^*,\mu)$的符号决定. 当这个导数在某个参数值处变为$0$时,系统可能发生分岔. 在高维系统$\vec x'(t)=F(\vec x,\mu)$中,一维导数$f_x(\vec x^*,\mu)$ 便由其Jacobi矩阵$D_xF(\vec x_*,\mu)$所取代. 此时平衡点的稳定性由该矩阵的特征值决定: 当特征值穿过虚轴时, 系统的相图结构可能发生改变. 比如当实特征值穿过$0$ 时可能出现鞍结分岔(saddle-node bifurcation); 当一对共轭复特征值穿过虚轴时可能出现Hopf分岔(Hopf bifurcation).

2. 考虑二维线性系统 $\vec y'(t)=\vec A\vec y(t), \quad A= \begin{pmatrix} a&b\\ c&d \end{pmatrix}.$ 令

$$p=\operatorname{tr}(A)=a+d, \qquad q=\det(A)=ad-bc, \qquad \Delta=p^2-4q,$$

(i) 将下列条件与平衡点$\vec y(t)=\vec 0$的稳定性配对:

$$\begin{array}{c|c} \text{条件} & \text{稳定性} \\ \hline \text{(1) } q>0,\ p<0 & \text{(A) 不稳定} \\[1mm] \text{(2) } p=0,\ q>0 & \text{(B) 渐近稳定} \\[1mm] \text{(3) } q<0\ \text{或}\ p>0 & \text{(C) Lyapunov 稳定但非渐近稳定} \end{array}$$

(ii) 将下列条件与平衡点$\vec y(t) = \vec 0$处相应的相图类型配对:

$$\begin{array}{c|c} \text{条件} & \text{相图类型} \\ \hline \text{(1) } q<0 & \text{(A) 中心点} \\[1mm] \text{(2) } q>0,\ \Delta>0 & \text{(B) 螺旋点} \\[1mm] \text{(3) } q>0,\ \Delta<0,\ p\neq 0 & \text{(C) 鞍点} \\[1mm] \text{(4) } q>0,\ p=0 & \text{(D) 结点} \end{array}$$

3. [带电粒子在匀强磁场中的运动] 设一带电粒子在$xy$平面内运动, 已知存在垂直于直面向外的匀强磁场, 且粒子的重力和空气阻力忽略不计.

(i) 若粒子在$x,y$方向上的速度分量$v_x(t), v_y(t)$满足

$$\begin{cases} v_x'(t) &= \omega v_y(t)\\ v_y'(t) &= - \omega v_x(t) \end{cases}, \quad \omega \text{为常数.}$$

求出$v_x(t), v_y(t)$的通解, 并判断平衡点的稳定性, 绘制出平衡点附近的相图;

(ii) 若粒子在$x,y$方向上的速度分量$v_x(t), v_y(t)$满足

$$\begin{cases} v_x'(t) &= -\alpha v_x(t) + \omega v_y(t)\\ v_y'(t) &= - \omega v_x(t) -\alpha v_y(t) \end{cases}, \quad \alpha,\omega \text{为常数.}$$

求出$v_x(t), v_y(t)$的通解, 并判断平衡点的稳定性, 绘制出平衡点附近的相图.