微分方程、动力系统与混沌导论 第6章 高维线性系统[书摘]

第6章　高维线性系统

      在线性代数领域稍作停留后，现在该回到微分方程了，特别地，要回到求解具有常系数的高维线性系统的任务中来。和线性代数那一章一样，我们要讨论大量的不同情形。

 

6.1 不同特征值

      首先考虑线性系统$\boldsymbol X' = \boldsymbol {AX}$，其中$n \times n$矩阵$\boldsymbol A$具有$n$个不同的实特征值$\lambda_1,\cdots,\lambda_n$。根据第5章中的结论，存在坐标变换$\boldsymbol T$使得新系统$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$具有如下相当简单的形式：

\[\begin{array}{l}{{y'}_1} = {\lambda _1}{y_1}\\\;\;\;\;\;\; \vdots \\{{y'}_n} = {\lambda _n}{y_n}\end{array}\]

其中的线性映射$T$将标准基向量$\boldsymbol E_j$映到属于$\lambda_j$的特征向量$\boldsymbol V_j$。显然，下面形式的函数

\[\boldsymbol Y(t)=\left( \begin{array}{l} c_1e^{\lambda_1t} \\  \vdots \\c_ne^{\lambda_nt}\end{array}\right)\]

是$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$满足初值条件$\boldsymbol Y(0) = (c_1,\cdots,c_n)$的一个解。与第3章一样，它是这样的唯一解。

       由此可得$\boldsymbol X(t) = \boldsymbol {TX}(t)$就是$\boldsymbol X' = \boldsymbol {AX}$的通解，而且这个通解可以写成如下形式

\[\boldsymbol X(t) = \sum \limits_{j=1}^nc_je^{\lambda_jt}\boldsymbol V_j.\]

现在假设$\boldsymbol A$的特征值中$\lambda_1,\cdots,\lambda_k$为负的，而$\lambda_{k+1},\cdots,\lambda_n$为正的。因为没有零特征值，故此时系统是双曲的。首先，所有从$\boldsymbol V_1,\cdots,\boldsymbol V_k$张成的子空间出发的解都将永远逗留在这个子空间中，这是因为此时$c_{k+1}=\cdots=c_n=0$。其次，每个这样的解在$t \to \infty$时都趋于原点。与平面系统中引入的术语类似，我们称这个子空间为稳定子空间。类似地，$\boldsymbol V_{k+1},\cdots,\boldsymbol V_n$张成的子空间中的解将远离原点。这个子空间称为不稳定子空间。所有其它的解在时间向后时都趋于稳定子空间，而在时间向前时都趋于不稳定子空间。因而这个系统对应于高维的鞍点。

例  考虑

\[\boldsymbol X' = \left( \begin{array}{l} 1&2&-1 \\ 0&3&-2\\0&2&-2 \end{array} \right) \boldsymbol X.\]

在5.2节中，我们已经看到这个矩阵的特征值为$2,1,-1$，属于它们的特征向量分别为$(3,2,1),(1,0,0),(0,1,2)$。于是矩阵

\[\boldsymbol T = \left( \begin{array}{l} 3&1&0 \\ 2&0&1 \\ 1&0&2 \end{array} \right) \]

就将 $\boldsymbol X' = \boldsymbol {AX}$化为

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y = \left( \begin{array}{l} 2&0&0 \\ 0&1&0 \\ 0&0&-1 \end{array} \right) \boldsymbol Y,\]

我们可以立刻得到它的解，用$\boldsymbol T$乘以这个解就可以得到$\boldsymbol X' = \boldsymbol {AX}$的通解(如果单为求通解的话，无需坐标变换，因为前面已经求得了特征值及对应的特征向量，从而可以直接写出通解，此处作变换的目的是为了化为标准形)为

\[\boldsymbol X(t) = {c_1}{e^{2t}}\left( \begin{array}{l}3\\2\\1\end{array} \right) + {c_2}{e^t}\left( \begin{array}{l}1\\0\\0\end{array} \right) + {c_3}{e^{ - t}}\left( \begin{array}{l}0\\1\\2\end{array} \right).\]

过原点以及$(0,1,2)$的直线就是稳定线，而由$(3,2,1),(1,0,0)$张成的平面是不稳定平面。该系统以及系统$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y$的解如图6.1所示。

例  如果$3 \times 3$矩阵$\boldsymbol A$有3个不同的负特征值，则我们可作坐标变换使得系统化成

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT})\boldsymbol Y = \left( \begin{array}{l} \lambda_1 &0&0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 &\lambda_3 \end{array} \right) \boldsymbol Y,\]

其中$\lambda_3 < \lambda_2 < \lambda_1 <0$。系统的所有解都趋向原点，这样我们就得到了一个高维汇点(见图6.2)。当初值条件$(x_0,y_0,z_0)$中的3个坐标都非零时，对应的解将沿切于$x$轴(特征值绝对值小的对应坐标)的方向趋于原点。

       现在，假设$n \times n$矩阵$\boldsymbol A$具有$n$个不同特征值，其中$k_1$个实的，$2k_2$个非实的，因而$n = k_1+2k_2$。于是，由第5章可知，存在坐标变换使得系统化成

\[\begin{align} x'_j & =\lambda_jx_j \\ u'_l & = \alpha_lu_l + \beta_lv_l \\ v'_l & =-\beta_lu_l + \alpha_lv_l \end{align},\]

这里$j=1,\cdots,k_1,l_1=1,\dots,k_2$。由第3章可知，我们有如下的解：

\[\begin{align} x_j(t) & =c_je^{\lambda_jt} \\ u_l(t) & =p_l e^{\alpha_lt}\cos\beta_lt + q_l e^{\alpha_lt}\sin\beta_lt \\ v_l(t) & =-p_l e^{\alpha_lt}\sin\beta_lt + q_l e^{\alpha_lt}\cos\beta_lt. \end{align}\]

像前面一样，我们可以直接验证这就是通解。于是，我们就证明了：

定理  考虑系统$\boldsymbol X' = \boldsymbol {AX}$，其中$\boldsymbol A$具有不同的特征值，$\lambda_1,\cdots,\lambda_{k_1} \in \mathbb R$，而$\alpha_1+\text i\beta_1,\cdots,\alpha_{k_2}+\text i\beta_{k_2} \in \mathbb C$。假设$\boldsymbol T$将$\boldsymbol A$化为标准形

\[\boldsymbol T^{-1}\boldsymbol {AT} = \left( \begin{array}{l} \lambda_1 \\&\ddots \\ && \lambda_{k_1} \\ &&& B_1 \\ &&&& \ddots \\ &&&&& B_{k_2} \end{array} \right),\]

其中

\[\boldsymbol B_j = \left( \begin{array}{l} \alpha_j & \beta_j \\ -\beta_j & \alpha_j \end{array} \right).\]

则$\boldsymbol X' = \boldsymbol {AX}$的通解为$\boldsymbol {TY}(t)$，而

 

\[\boldsymbol Y(t) = \begin{pmatrix} {c_1}{e^{{\lambda _1}t}}\\ \vdots \\{c_{{k_1}}}{e^{{\lambda _{{k_1}}}t}}\\{a_1}{e^{{\alpha _1}t}}\cos {\beta _1}t + {b_1}{e^{{\alpha _1}t}}\sin {\beta _1}t\\ - {a_1}{e^{{\alpha _1}t}}\sin {\beta _1}t + {b_1}{e^{{\alpha _1}t}}\cos {\beta _1}t\\ \vdots \\{a_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\cos {\beta _{{k_2}}}t + {b_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\sin {\beta _{{k_2}}}t\\ - {a_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\sin {\beta _{{k_2}}}t + {b_{{k_2}}}{e^{{\alpha _{{k_2}}}t}}\cos {\beta _{{k_2}}}t \end{pmatrix}.\]

与通常一样，矩阵$\boldsymbol T$的列向量由对应于每个特征值的特征向量(或特征向量的实部和虚部)构成。依然与前面一样，具有负实部(正实部)的特征值对应的特征向量张成的子空间为稳定(不稳定)子空间。

例  考虑系统

\[\boldsymbol X' = \begin{pmatrix} 1&0&1\\ -1&0&0 \\ 0&0&-1 \end{pmatrix} \boldsymbol X,\]

这里的矩阵已经是标准形，其特征值为$\pm \text i,-1$。满足初值条件$(x_0,y_0,z_0)$的解为

\[\boldsymbol Y(t) = x_0 \begin{pmatrix} \cos t \\ -\sin t \\ 0 \end{pmatrix} + y_0 \begin{pmatrix} \sin t \\ -\cos t \\ 0 \end{pmatrix} + z_0e^{-t} \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix},\]

这就是系统的通解。系统的相图如图6.3所示。稳定线为沿$z$轴的直线，而$xy$平面上的所有解都绕中心为原点的圆周旋转。事实上，其它不在稳定线上的每一个解都位于由“$x^2 + y^2$ =常数 ”所确定的$\mathbb R^3$中的圆柱面上，当初值条件$(x_0,y_0,z_0)$中的$z_0 \ne 0$时，这些解都盘旋地趋向$xy$平面上以$\sqrt {x_0^2 + y_0^2}$为半径的圆形解。

 

例  现在考虑系统 $\boldsymbol X' = \boldsymbol {AX}$，其中

\[\boldsymbol A = \begin{pmatrix} -0.1 & 0 &1 \\ -1 & 1 & –1.1 \\ -1 & 0 & –0.1 \end{pmatrix}.\]

显然它的特征方程为

\[-\lambda^3 + 0.8\lambda^2 -0.81\lambda + 1.01 = 0,\]

求得特征值分别为$1,-0.1\pm \text i$。求解方程$(\boldsymbol A - (0.1 + \text i) \boldsymbol I ) \boldsymbol X = \boldsymbol 0$可得属于$-0.1 + \text i$的特征向量$(-\text i,1,1)$。记$\boldsymbol V_1 = \text {Re} (-\text i,1,1) = (0,1,1), \boldsymbol V_2 = \text {Im} (-\text i,1,1) = (-1,0,0)$。再求解方程$(\boldsymbol A - \boldsymbol I)\boldsymbol X = \boldsymbol 0$可得属于$\lambda = 1$的特征向量$\boldsymbol V_3 = (0,1,0)$。从而以$\boldsymbol V_i$为列向量的矩阵

\[\boldsymbol T = \begin{pmatrix} 0&-1&0 \\ 1&0&1 \\ 1 & 0 &0  \end{pmatrix}\]

就将 $\boldsymbol X' = \boldsymbol {AX}$化为标准形

\[\boldsymbol Y' = \begin{pmatrix} –0.1 & 1 & 0 \\ -1 & –0.1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \boldsymbol Y.\]

这个系统的不稳定线为$z$轴，而$xy$平面为稳定平面。注意，稳定平面上的所有解都盘旋进入原点。我们称这个系统为一个螺线鞍点(见图6.4)。稳定平面外的典型解都将盘旋趋向于$z$轴，而其$z$坐标增加或减小(原文和译文似乎都有点问题？)(见图6.5)。

 

 

 

6.2 调和振子

      考虑一对无阻尼的调和振子，它们的方程是：

\[\begin{array}{l} x''_1 = -\omega_1^2x_1 \\ x''_2 = -\omega_2^2x_2. \end{array} \]

通过检验脑海中的$\sin \omega t,\cos \omega t$这样的函数，我们就几乎可以解出这两个方程。然而，让我们再进一步，首先对上一节中的定理在非实特征值情形做些解释，但更重要的是引入一些有趣的几何。

      先引入两个新变量$y_j=x'_j,j=1,2$，这样方程组就可以改写成下面的系统

\[\begin{align} x'_1 &= y_j \\ y'_j &= -\omega_j^2x_j. \end{align} \]

该系统的矩阵形式为$\boldsymbol X' = \boldsymbol {AX}$，其中$\boldsymbol X = (x_1,y_1,x_2,y_2)$，

\[\boldsymbol A = \begin{pmatrix} 0&1 && \\ -\omega_1^2 & 0 && \\ && 0 & 1 \\ && -\omega_2^2 & 0 \end{pmatrix}.\]

这个系统的特征值为$\pm \text i \omega_1, \pm \text i \omega_2$。属于$ \text i \omega_1, \text i \omega_2$的特征向量分别为$\boldsymbol V_1 = (1,\text i \omega_1,0,0),\boldsymbol V_2 = (0,0,1,\text i \omega_2)$。记$\boldsymbol W_1,\boldsymbol W_2$分别为$\boldsymbol V_1$的实部和虚部，$\boldsymbol W_3,\boldsymbol W_4$分别为$\boldsymbol V_2$的实部和虚部。像通常业样，令$\boldsymbol {TE}_j = \boldsymbol W_j$，则线性映射$\boldsymbol T$就将系统化为标准形，对就的矩阵为

\[\boldsymbol T_{-1}\boldsymbol {AT} = \begin{pmatrix} 0&\omega_1 && \\ -\omega_1 & 0 && \\ && 0 & \omega_2 \\ && -\omega_2 & 0 \end{pmatrix}.\]

这样就可以得到通解(见上一节的定理)

\[\boldsymbol Y(t) = \begin{pmatrix} x_1(t) \\ y_1(t) \\ x_2(t) \\ y_2(t) \end{pmatrix} = \begin{pmatrix} a_1\cos\omega_1t+b_1\sin\omega_1t \\ -a_1\sin\omega_1t+b_1\cos\omega_1t \\ a_2\cos\omega_2t+b_2\sin\omega_2t \\ -a_2\sin\omega_2t+b_2\cos\omega_2t \end{pmatrix}, \]

它的形式正像我们所想像的那样。

      我们也许会认为事情就至此为止了，因为解的表达式都已经得到了，但是，我们将往前进一步。

      显然，每对解$(x_j(t),y_j(t)),j=1,2$都是单个方程的一个以$2\pi / \omega_j$为周期的周期解，但这并不意味着完整的四维解也是一个周期解。事实上，完整的解是一个周期为$\tau$的周期函数当且仅当存在整数$m,n$使得

\[\omega_1 \tau = m·2\pi, \;\;\;\; \omega_2 \tau = n·2\pi .\]

从而，为了得到周期性，我们必须有

\[\tau = \frac{2\pi m}{\omega_1} =\frac {2\pi n}{\omega_2}, \]

这等价于

\[\frac{\omega_2}{\omega_1} = \frac{n}{m}. \]

即，调和振子的两个频率比必须是一个有理数。在图6.6中，我们画出了当频率比为5/2时系统的一个特解$(x_1(t),x_2(t))$(这是一个平面相图，可以这么想像，有一个沙漏沿$x_1$轴来回摆动，同时沿$x_2$轴来回摆动，合成的相图即如图所示)。

      当频率之比为无理数时，情况非常不一样。为了理解这点，我们作另外一个(大家非常熟悉的)坐标变换。用标准形来写，现在的系统是

\[\begin{align} x'_j &= \omega_j y_j \\ y'_j &= -\omega_jx_j \end{align}. \]

我们引入极坐标$(r_j,\theta_j)$来代替变量$x_j,y_j$。对式子

\[r_j^2 = x_j^2 + y_j^2\]

两边求导得得

\[2r_j r'_j = 2x_jx'_j + 2y_jy'_j = 2x_jy_j\omega_j - 2x_jy_j\omega_j = 0.\]

于是，对每一个$j$都有$r'_j = 0$。同样地，对方程

\[\tan \theta_j  = \frac{y_j}{x_j}\]

求导可得

\[(\sec^2 \theta_j) \theta'_j = \frac{y'_jx_j - y_jx'_j}{x_j^2} = \frac{-\omega_jr_j^2}{r_j^2\cos^2 \theta_j},\]

由此可得

\[\theta'_j = -\omega_j.\]

从而在极坐标下，这些方程就变得非常简单：

\[\begin{align} r'_j & =0 \\ \theta'_j &= -\omega_j \end{align}. \]

第一个方程告诉我们，$r_1,r_2$沿任何解都始终为常数。从而，无论我们如何选取$r_1,r_2$的初始值，关于$\theta_j$的方程都是一样的。这样我们就可以将注意力集中到$r_1 = r_2 =1$的情形。在$\mathbb R^4$中得到的点集是一个环面，也就是油炸圈饼的表面(挺形象了，试想一只蚂蚁既可以沿横截面的一个圆作周期运动，也可以沿纵截面的一个圆作周期运动)。虽然在四维空间中这有点难以想像，但不管怎样，这个集合上有2个独立的变量，即$\theta_1,\theta_2$，而且它们都以$2\pi$为周期的。在这一点上，这类似于$\mathbb R^3$中熟知的环面可以被2个独立的循环走向所参数化(蚂蚁可以走到油炸圈饼的任何一点)。

      现在，方程限制在环面上就变成了：

\[\begin{align}\theta'_1  &= -\omega_1 \\ \theta'_2 & = -\omega_2. \end{align} \]

为了方便，我们可将$\theta_1,\theta_2$看成是边长为$2\pi$的正方形中的变量(只要将正方形的对边$\theta_j = 0$和$\theta_j = 2\pi$粘合起来就得到了环面——注：这已经是一种很好的处理方式了，尽管很难将一个正方形变成四维中的环)。在正方形中，向量场的斜率为常数

\[\frac{\theta'_2}{\theta'_1} = \frac{\omega_2}{\omega_1}.\]

于是，在正方形中，解都位于斜率为$\omega_2 / \omega_1$的直线上。当解到达$\theta_1 = 2\pi$这条边时(不妨设在$\theta_2 = c$处)，它就在$\theta_1 = 0$的边马上重新出现(($\theta_2$的坐标依然为$c$)，然后以斜率$\omega_2 / \omega_1$再继续向前。在解到达$\theta_2 = 2\pi$时，将会出现类似的等同。

      这样，我们就对解在环面上的行为有了一个简化的几何图像(什么时候，国内的书能写得如此简明而透彻呢?)。将会出现什么呢？答案依赖于$\omega_2 / \omega_1$的取值。如果这个值为有理数，比如，$n/m$，则解从$(\theta_1(0),\theta_2(0))$出发将垂直经过环面$n$次、水平经过环面$m$次之后再回到出发点。这就是我们上面观察到的周期解。有时，直线解在$\theta_1 \theta_2$平面上的图像和图6.6中所画出的$x_1x_2$平面上的图像并不完全一样。

      在无理情形，将会出现一些很不同的情况(见图6.7)。为了理解现在出现的情况，让我们回到第1章讨论过的庞加莱映射的概念。考虑圆周$\theta_1 = 0$，也就是环面正方形表示中的左边。任给该圆周上的一个初值，比如，$\theta_2 = x_0$，我们跟随从该点出发的解一直到它下次再到达$\theta_1 = 2\pi$。由于环面是等同正方形的对边得到的，解现在实际上就是回到了出发的圆周$\theta_1 = 0$。在这个进行过程中，解有可能会穿过边界$\theta_2 = 2\pi$几次，但它最终会回到$\theta_1 = 0$。因而我们可以定义$\theta_1 = 0$上的庞加莱映射如下：庞加莱映射在圆周上$x_0$处的取值就是首次返回点的相应坐标。假设首次返回出现在点$\theta_2(\tau)$处，其中$\tau$是满足$\theta_1(\tau) = 2\pi$的时刻。由于$\theta_1(t) = \theta_1(0) + \omega_1t$，因而$\tau = 2\pi / \omega_1$。从而$\theta_2(t) = x_0 + \omega_2(2\pi / \omega_1)$。于是在这个圆周上的庞加莱映射可以写成

\[f({x_0}) = {x_0} + 2\pi ({\omega _2}/{\omega _1})\;\bmod \;2\pi \]

其中$x_0 = \theta_2(0)$是初值在这个圆周上的$\theta_2$坐标(见图6.8)(上式可理解为初始位置在$x_0$，步长为$2\pi ( \omega_2 / \omega_1)$的一种运动，只是要进行模$2\pi$，并依次记录点的轨迹，因此也是一种迭代过程)。于是，该圆周上的庞加莱映射就是将圆周上的点旋转角度$2\pi( \omega_2 / \omega_1)$这个函数。由于$\omega_2 / \omega_1$是无理的，这个函数称为圆周上的无理旋转。

定义  点集$x_0,x_1 = f(x_0),x_2 = f(f(x_0)),\cdots, x_n = f(x_{n-1})$称为$x_0$在$f$迭代下的轨道。

      $x_0$的轨道记录了当时间增加时，解是如何穿过$\theta_1 = 2\pi$的(言下之意是另一个变量的状态是什么样的)。

命题  假设$\omega_2 / \omega_1$是无理数，则圆周$\theta_1 = 0$上任何初值$x_0$的轨道在该圆周上都是稠密的(意思是几乎充满整个圆周)。

证明详见书本。

      因为$x_0$的轨道在圆周$\theta_1 = 0$上是稠密的，从而在正方形中将这些点连接起来的直线解也是稠密的，于是原来的解在其所在的环面上不稠密的。这就解释了为什么当$\omega_2 / \omega_1 = \sqrt2$时，解如图6.7所示到$x_1x_2$平面时会稠密地布满。

      回到振子的实际运动，我们就得到当$\omega_2 / \omega_1$为无理数时，两个质点的运动并不是周期的。然而，由于解在环面上的稠密性，当时间增加时，它们又的确会一次次地回到非常靠近初始位置的地方。这种形式的运动称为拟周期运动。

 

6.3 重特征值

      在上一章中，我们已经看到，具有实重特征值系统的求解归结为求解其矩阵含有如下形式的块的系统(如果存在不唯一解的重特征值呢？)：

\[\left( \begin{array}{l} \lambda & 1 \\ & \lambda & 1 \\ & & \ddots \ddots \\ & & & \ddots &1 \\ & & & & \lambda \end{array} \right). \]

例  令

\[\boldsymbol X' = \begin{pmatrix} \lambda & 1 & 0 \\ 0 & \lambda & 1 \\ 0 & 0 & \lambda \end{pmatrix} \boldsymbol X.\]

这个系统唯一的特征值就是$\lambda$，而且唯一的特征向量(1,0,0)。我们像在第3章所做的那样来求解这个系统。首先，由于$x'_3 = \lambda x_3$，故必有

\[x_3(t) = c_3e^{\lambda t},\]

进而有，

\[x'_2 = \lambda x_2 + c_3e^{\lambda t}.\]

与第3章一样，我们猜测它有如下形式的解：

\[x_2(t) = c_2e^{\lambda t}+ \alpha t e^{\lambda t}.\]

将上式代入$x'_2$的微分方程，可以定出$\alpha = c_3$，从而找到

\[x_2(t) = c_2e^{\lambda t}+ c_3 t e^{\lambda t}.\]

最后，方程

\[x'_1 = \lambda x_1 + c_2e^{\lambda t} + c_3 t e^{\lambda t}\]

让人联想起下面形式的解：

\[x_1(t) = c_1 e^{\lambda t} + \alpha t e^{\lambda t} + \beta t^2 e^{\lambda t}. \]

像刚才一样的求解可得

\[x_1(t) = c_1 e^{\lambda t} + c_2 t e^{\lambda t} +  c_3 \frac {t^2}{2} e^{\lambda t}. \]

总之，我们找到了

\[\boldsymbol X(t) = c_1 e^{\lambda t} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + c_2 e^{\lambda t} \begin{pmatrix} t \\ 1 \\ 0 \end{pmatrix} + c_3 e^{\lambda t} \begin{pmatrix} t^2/2 \\ t \\ 1 \end{pmatrix} \]

它就是系统的通解。虽然解的表达式中出现了多项式项，但当$\lambda < 0$时，指数将起主导作用，而且所有的解都趋于零。当$\lambda < 0$时，一些有代表性的解如图6.9所示。注意，该系统只有唯一的直线解(关于直线解的概念，可以参考2.3节和2.6节，对于实特征值对应的特征向量就是直线解)，该解位于$x$轴上，而且$xy$平面是不变的，其上解的形态与平面重特征值情形完全相同。

例  考虑下面的四维系统：

\[\begin{align} x'_1 &= x_1 + x_2 - x_3 \\ x'_2 &= x_2 + x_4 \\ x'_3 &= x_3 + x_4 \\ x'_4 &= x_4 \end{align} \]

我们可以将该系统写成如下的矩阵形式：

\[\boldsymbol X' = \boldsymbol {AX} = \begin{pmatrix} 1&1&-1&0 \\ 0&1&0&1 \\ 0&0&1&1 \\ 0&0&0&1 \end{pmatrix} \boldsymbol X.\]

因为$\boldsymbol A$为上三角的，故它的所有特征值都是1。求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol 0$可找到两个线性无关的特征向量$\boldsymbol V_1 = (1,0,0,0),\boldsymbol W_1 = (0,1,1,0)$。这就将$\boldsymbol A$可能的标准形减少到两种。进一步求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol V_1$可得一个解$\boldsymbol V_2 = (0,1,0,0)$，求解$(\boldsymbol A - \boldsymbol I) \boldsymbol X = \boldsymbol W_1$可得一个解$\boldsymbol W_2 = (0,0,0,1)$(这种方法的理论可参考5.5节的内容)。由此，我们知道系统$\boldsymbol X' = \boldsymbol {AX}$可以变换成

\[\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT}) \boldsymbol Y = \begin{pmatrix} 1&1&0&0 \\ 0&1&0&0 \\ 0&0&1&1 \\ 0&0&0&1 \end{pmatrix} \boldsymbol Y,\]

其中矩阵$\boldsymbol T$由下式给出：

\[\boldsymbol T = \begin{pmatrix} 1&0&0&0 \\ 0&1&1&0 \\ 0&0&1&0 \\ 0&0&0&1 \end{pmatrix} \]

于是$\boldsymbol Y' = (\boldsymbol T^{-1} \boldsymbol {AT}) \boldsymbol Y$的解为(重特征值情形下求通解可参考3.3节的内容 )

\[\begin{align} y_1(t) &= c_1 e^t + c_2 t e^t \\ y_2(t) &= c_2 e^t \\ y_3(t) &= c_3 e^t + c_4 t e^t  \\ y_4(t) &= c_4 e^t. \end{align} \]

再用坐标变换$\boldsymbol T$作用，我们就得到了原系统的通解为

\[\begin{align} x_1(t) &= c_1 e^t + c_2 t e^t \\ x_2(t) &= c_2 e^t + c_3 e^t + c_4 t e^t \\ x_3(t) &= c_3 e^t + c_4 t e^t  \\ x_4(t) &= c_4 e^t. \end{align} \]

 

6.4 矩阵指数

       现在，我们转到矩阵指数求解线性系统，这是解线性系统的一个不同但简练的方法。从某种意义上讲，这是攻克这类系统自然的方法。

      我们先回忆1.1节中是如何求解线性方程$x' = ax$的$1 \times 1$“系统”，其中我们的矩阵就是简单的$(a)$。那时我们并没有进行找特征值和特征向量的过程，(然而，事实上，我们进行了，但过程实在是太简单了。)而是直接取矩阵$(a)$的指数函数就找到了通解$x(t)  = c \text {exp} (at)$。事实上，这个过程对一般$n \times n$矩阵$\boldsymbol A$也是行得通的。我们唯一需要知道的就是如何求一个矩阵的指数。

      办法如下。由微积分可知，指数函数可以表示成一个无穷级数

\[e^x = \sum \limits_{k=0}^\infty \frac {x^k}{k!},\]

而且我们知道对每个$x \in \mathbb R$这个级数都收敛。由于我们可以将矩阵相加，可以取它们的$k$次方，还可以将矩阵的每个元素乘上$1/k!$，这将意味着我们一样可以用上面的级数去取它们的指数。

定义  设$\boldsymbol A$为一个$n \times n$矩阵。我们定义$\boldsymbol A$的指数是由下式给出的矩阵：

\[\text {exp} (\boldsymbol A) = \sum \limits_{k=0}^\infty \frac {\boldsymbol A^k}{k!}.\]

      当然，我们还要担心上面矩阵和式收敛的含义。

例  令

\[\boldsymbol A = \begin{pmatrix} \lambda & 0 \\ 0 & \mu \end{pmatrix} . \]

则我们有

\[\boldsymbol A^k =\begin{pmatrix} \lambda^k & 0 \\ 0 & \mu^k \end{pmatrix},\]

于是

\[\text {exp}(\boldsymbol A) = \begin{pmatrix} \sum \limits_{k=0}^ \infty \lambda^k / k!  & 0 \\ 0 & \sum \limits_{k=0}^\infty \mu^k / k! \end{pmatrix}  = \begin{pmatrix} e^\lambda & 0 \\ 0 & e^\mu \end{pmatrix},\]

这个结果大家可能已经事先猜到。

例  下面来看一个稍等复杂些的例子。令

\[\boldsymbol A = \begin{pmatrix} 0 & \beta \\ -\beta & 0 \end{pmatrix}.\]

由计算可得

\[\text {exp} (\boldsymbol A) =\begin{pmatrix} \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k}}{(2k)!} & \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k+1}}{(2k+1)!} \\ -\sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k+1}}{(2k+1)!} & \sum \limits_{k=0}^\infty (-1)^k \frac{\beta^{2k}}{(2k)!} \end{pmatrix} = \begin{pmatrix} \cos \beta & \sin \beta \\ - \sin\beta & \cos\beta \end{pmatrix}.\]

例  现在令

\[\boldsymbol A = \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix},\]

其中$\lambda \ne 0$。为了照顾随后的应用，我们来计算$\text {exp} (t \boldsymbol A)$，而不是$\text {exp} A)$。我们有

\[(t \boldsymbol A)^k = \begin{pmatrix} (t\lambda)^k & kt^k\lambda^{k-1} \\ 0 & (t\lambda)^k \end{pmatrix}.\]

于是我们得到

\[\text {exp}(t \boldsymbol A)  = \begin{pmatrix} \sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!} & t\sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!}  \\ 0 & \sum \limits_{k=0}^\infty \frac{(t\lambda)^k}{k!} \end{pmatrix}= \begin{pmatrix} e^{t\lambda} & t e^{t\lambda} \\ 0 & e^{t\lambda}.\end{pmatrix} \]

      可以看到，在上面3个例子中，矩阵$\text {exp}(\boldsymbol A)$的元素都是一些无穷级数。因此，我们称矩阵无穷级数$\text {exp}(\boldsymbol A)$绝对收敛，如果每一个单独的项绝对收敛。在上面的几种情形中，收敛性是显然的。遗憾的是，对于一般的矩阵$\boldsymbol A$，收敛性并不清楚。为了证明这里的收敛性，我们需要做一些努力。

      记$a_{ij}(k)$为$\boldsymbol A^k$的$ij$元素。令$a = \text {max}|a_{kj}|$。我们有

\[|a_{ij}(k)| \le n^{k-1}a^k. \]

于是，我们得到了$n \times n$矩阵$\text {exp}(\boldsymbol A)$的$ij$元素的一个上界：

\[\Bigg| \sum \limits_{k=0}^\infty \frac{a_{ij}(k)}{k!} \Bigg| \le \sum \limits_{k=0}^\infty \frac{|a_{ij}(k)|}{k!} \le \sum \limits_{k=0}^\infty \frac{n^{k-1}a^k}{k!} \le \text {exp}(na),\]

因此，根据比较判别法，这个级数绝对收敛。这样，对任何$\boldsymbol A \in L(\mathbb R^n),\text {exp}(\boldsymbol A)$有意义。

      下面的结果表明，矩阵指数也具有通常指数函数的许多熟悉的性质。

命题  设$\boldsymbol A,\boldsymbol B,\boldsymbol T$都是$n \times n$矩阵，则
      (1)如果$\boldsymbol B=\boldsymbol T^{-1} \boldsymbol {AT}$，则$\text {exp}(\boldsymbol B) = \boldsymbol T^{-1} \text {exp}(\boldsymbol A) \boldsymbol T$。
      (2)如果$\boldsymbol {AB} = \boldsymbol {BA}$，则$\text {exp}(\boldsymbol A+\boldsymbol B) = \text {exp} (\boldsymbol A)\text {exp} (\boldsymbol B)$。
      (3)$\text {exp}(- \boldsymbol A) = (\text {exp} (\boldsymbol A))^{-1}$。

证明详见书本。注意，命题中论断(3)意味着，对每一个矩阵$\boldsymbol A,\text {exp}(\boldsymbol A)$都可逆。这一点类似于对每个实数$a,e^a \ne 0$。

      在$\boldsymbol A$和$\text {exp}(\boldsymbol A)$的特征向量之间有一个很简单的关系：

命题  如果$\boldsymbol V \in \mathbb R^n$是$\boldsymbol A$的属于特征值$\lambda$的特征向量，则$\boldsymbol V$也是$\text {exp} (\boldsymbol A)$的特征向量，它属于$e^\lambda$。

       现在我们回到微分方程系统的讨论。令$\boldsymbol A$为$n \times n$矩阵，考虑系统$\boldsymbol X' = \boldsymbol {AX}$。回想我们以前用$L(\mathbb R^n)$记全体$n \times n$矩阵的集合。把每个$t \in \mathbb R$映到$\text {exp}(t\boldsymbol A)$就得到了一个函数$\mathbb R \to L(\mathbb R^n)$。由于我们将$L(\mathbb R^n)$等同于$\mathbb R^{n^2}$，因而谈论这个函数的导数是有意义的。

命题

\[\frac{\text d}{\text d t} \text {exp}(t \boldsymbol A) = \boldsymbol A \text {exp}(t \boldsymbol A) = \text {exp}(t \boldsymbol A) \boldsymbol A. \]

换句话说，矩阵值函数$t \to \text {exp}(t \boldsymbol A)$的导数是另一个矩阵值函数$\boldsymbol A \text {exp} (t \boldsymbol A)$。

      现在我们回到求解微分方程系统。下面的结论可以看成是具有常系数的线性微分方程组的基本定理。

定理  设$\boldsymbol A$为$n \times n$矩阵，则初值问题$\boldsymbol X' = \boldsymbol {AX},\boldsymbol X(0) = \boldsymbol X_0$的解为$\boldsymbol X(t) = \text {exp}(t \boldsymbol A)\boldsymbol X_0$，而且它是唯一解。

例  考虑系统

\[\boldsymbol X' = \begin{pmatrix} \lambda &1 \\ 0 & \lambda \end{pmatrix} \boldsymbol X.\]

根据上面的定理，它的通解为

\[\boldsymbol X(t) = \text {exp}(t \boldsymbol A) \boldsymbol X_0 = \text {exp} \begin{pmatrix} t \lambda & t \\ 0 & t \lambda \end{pmatrix} \boldsymbol X_0.  \]

这其中的矩阵指数在前面已经算过，于是我们有

\[\boldsymbol X(t) = \begin{pmatrix} e^{t\lambda} & t e^{t\lambda} \\ 0 & e^{t\lambda}\end{pmatrix} \boldsymbol X_0.\]

可以看出，这与第3章中的计算结果相吻合。

 

6.5 非自治线性系统

      直到现在，几乎所有碰到的线性微分方程系统都是自治的。但是，在应用中还是经常会出现一些形式的非自治系统。其中之一就是如下形式的系统

\[\boldsymbol X' = \boldsymbol A(t) \boldsymbol X,\] 

其中$\boldsymbol A(t) = [a_{ij}(t)]$为$n \times n$矩阵，它连续地依赖于时间。当以后章节中遇到变分方程时，我们会再深入研究这类系统。

      现在我们只关系另一个不同的非自治线性系统：

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\] 

其中$\boldsymbol A$为常值的$n \times n$矩阵，$\boldsymbol G:\mathbb R \to \mathbb R^n$显式地依赖于$t$，我们称之为强迫项。这是一个一阶线性非自治方程系统。

例  (受迫调和振子)如果调和振子系统中有一个外力作用，决定运动的微分方程就变成

\[x'' + bx' + kx = f(t),\]

其中$f(t)$为外力的大小。一个重要的特殊情形是其中的外力是时间的周期函数，这可对应于，比方说，将质点-弹簧装置所在的桌子前后周期地移动。作为一个系统，受迫调和振子方程变成

\[\boldsymbol X' = \begin{pmatrix} 0 & 1 \\ -k & -b \end{pmatrix} \boldsymbol X + \boldsymbol G(t),\;\;\; \boldsymbol G(t) = \begin{pmatrix} 0 \\ f(t) \end{pmatrix}.\]

      对一个非齐次的系统，如果将其中的时间项扔掉后，所得的系统$\boldsymbol X' = \boldsymbol {AX}$称为齐次方程。我们已经知道如何找到这类系统的通解。利用上节中的记号，齐次系统满足初值条件$\boldsymbol X(0)= \boldsymbol X_0$的解为

\[\boldsymbol X(t) = \text {exp}(t \boldsymbol A) \boldsymbol X_0,\]

而这也是齐次方程的通解。

      为了找到非齐次方程的通解，假设我们已经有了该方程的一个特解$\boldsymbol Z(t)$。于是$\boldsymbol Z'(t) = \boldsymbol {AZ}(t) + \boldsymbol G(t)$。当$\boldsymbol X(t)$为齐次方程的任一个解时，则函数$\boldsymbol Y(t) = \boldsymbol X(t) + \boldsymbol Z(t)$就是非齐次方程的另一个解。这点可由如下的计算得出：

\[\begin{align} \boldsymbol Y' = \boldsymbol X' + \boldsymbol Z' &=\boldsymbol {AX} + \boldsymbol {AZ} +\boldsymbol G(t) \\ &=\boldsymbol A(\boldsymbol X+\boldsymbol Z) +\boldsymbol G(t) \\ &=\boldsymbol {AY}+\boldsymbol G(t).\end{align} \]

因而，由于已经了齐次方程的所有解，一旦找到了非齐次方程的哪怕只是一个特解，我们就马上可以找到非齐次方程的通解。通常，我们可以通过猜解得到这样一个解(微积分中，这种方法通常称为待定系数法)。然而 ，猜解并不总是可行。下面的参数变易法却总是可行的，当然，这并不能保证我们总能计算其中的积分。

定理   (参数变易法)考虑非齐次方程

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\] 

其中$\boldsymbol A$为常值的$n \times n$矩阵，$\boldsymbol G(t)$为$t$的连续函数。则

\[\boldsymbol X(t) = \exp (t \boldsymbol A)\left( {{\boldsymbol X_0} + \int_0^t {\exp ( - s\boldsymbol A) \boldsymbol G(s) \text ds} } \right)\]

是该方程满足$\boldsymbol X(0) = X_0$的一个解。(对方程两边求导即可验证)

      下面我们来给出这个结果在周期受迫调和振子中的几个应用。首先，假设有个外力为$\cos t$的阻尼振子，此时外力项的周期为$2\pi$。这个系统为

\[\boldsymbol X' = \boldsymbol {AX} + \boldsymbol G(t),\]  

其中$\boldsymbol G(t) = (0,\cos t), \boldsymbol A$为矩阵

\[\boldsymbol A = \begin{pmatrix} 0&1 \\ -k & -b \end{pmatrix},\]

其中$b,k>0$。我们断言：这个系统有唯一一个以$2\pi$为周期的周期解。为了证明这个断言，我们首先需要找到一个满足条件$\boldsymbol X(0) = \boldsymbol X_0 =\boldsymbol X(2\pi)$的解。由参数变易法，我们需要找到$\boldsymbol X_)$使得

\[{\boldsymbol X_0} = \exp (2\pi \boldsymbol A){\boldsymbol X_0} + \exp (2\pi \boldsymbol A)\int_0^{2\pi } {\exp ( - s\boldsymbol A)\boldsymbol G(s)\text ds.} \]

现在其中一项

\[\exp (2\pi \boldsymbol A)\int_0^{2\pi }  \exp ( - s\boldsymbol A)\boldsymbol G(s)\text ds\]

为常值向量，将它记为$\boldsymbol W$。这样我们就要求解方程

\[(\exp (2\pi \boldsymbol A) - \boldsymbol I)\boldsymbol X_0 = -\boldsymbol W.\]

由于$\exp (2\pi \boldsymbol A) - \boldsymbol I$可逆，这个方程有唯一解。为了证明这个矩阵的可逆性，假设它是不可逆的，则有非零向量$\boldsymbol V$使得

\[(\exp (2\pi \boldsymbol A) - \boldsymbol I) \boldsymbol V = \boldsymbol 0,\]

或者，换句话说，矩阵$\exp (2\pi \boldsymbol A)$有一个特征值为1。但是，由上节可知，$\exp (2\pi \boldsymbol A)$的特征值都由$\exp (2\pi \lambda_j)$给出，其中$\lambda_j$为$\boldsymbol A$的特征值。由于每个$\lambda_j$的实部都小于0，从而$\exp (2\pi \lambda_j)$小于1，于是矩阵$\exp (2\pi \boldsymbol A) - \boldsymbol I$的确是可逆的。这样导致$2\pi$周期解的唯一初值(这多少有点出人意料，是不是周期解不仅依赖强迫项，同时还依赖初值，并且这样的初值还是唯一的!!)就是

\[\boldsymbol X_0 = (\exp (2\pi \boldsymbol A) - \boldsymbol I)^{-1}(-\boldsymbol W).\]

      令$\boldsymbol X(t)$为满足$\boldsymbol X(0) = \boldsymbol X_0$的周期解。这个解称为稳态解。如果$\boldsymbol Y_0$是其它的初值条件，我们可将它写成$\boldsymbol Y_0 = (\boldsymbol Y_0 -\boldsymbol X_0) +\boldsymbol X_0$，从而过$\boldsymbol Y_0$的解为

\[\begin{align} \boldsymbol Y(t) &= \exp (t \boldsymbol A)({ \boldsymbol Y_0} - { \boldsymbol X_0}) + \exp (t \boldsymbol A){ \boldsymbol X_0} + \exp (t \boldsymbol A)\int_0^t {\exp ( - s \boldsymbol A) \boldsymbol G(s)\text ds} &= \exp (t \boldsymbol A)({ \boldsymbol Y_0} - { \boldsymbol X_0}) + \boldsymbol X(t). \end{align}\]

对于表达式中的第一项，由于它是齐次方程的一个解，因而在$t \to \infty$时，趋于0(只是针对本例而言，因为是带阻尼的) 。于是，这个系统的每个解在$t \to \infty$时都趋于稳态解。在物理上，这一点很明显的：带阻尼的(非受迫)调和振子的运动都趋于平衡态，剩下的运动就是由周期外力引起的。这样，我们就证明了：

定理   考虑受迫带阻尼的调和振子方程

\[x'' + bx' + kx = \cos t,\]

其中$,b>0$。则这个方程的所有解都趋于一个以$2\pi$为周期的稳态解。

      现在考虑一个受迫的无阻尼调和振子的特例：

\[\boldsymbol X' = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \boldsymbol X + \begin{pmatrix} 0 \\ \cos\omega t \end{pmatrix},\]

其中，现在外力的周期为$2\pi/\omega,\omega \ne \pm 1$。记

\[\boldsymbol X = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}.\]

齐次方程的解为

\[\boldsymbol X(t) = \exp (t \boldsymbol A) \boldsymbol X_0 = \begin{pmatrix} \cos t & \sin t \\ -\sin t& \cos t \end{pmatrix} \boldsymbol X_0.\]

由参数变易法可知，非齐次方程从原点出发的解为

\[\begin{align} \boldsymbol Y(t) &= \exp (t\boldsymbol A)\int_0^t \exp ( - s\boldsymbol A) \begin{pmatrix} 0\\ \cos \omega s\end{pmatrix}\text ds \\ &= \exp (t\boldsymbol A)\int_0^t \begin{pmatrix} \cos s &- \sin s \\ \sin s &\cos s \end{pmatrix} \begin{pmatrix} 0\\ \cos \omega s \end{pmatrix} \text ds\\ &= \exp (t\boldsymbol A) \int_0^t \begin{pmatrix} - \sin s\cos \omega s \\ \cos s\cos \omega s \end{pmatrix} \text ds \\ &= \frac{1}{2} \exp (t\boldsymbol A)\int_0^t \begin{pmatrix} \sin (\omega - 1)s - \sin (\omega + 1)s \\ \cos (\omega - 1)s + \cos (\omega + 1)s \end{pmatrix} \text ds. \end{align} \]

回想到

\[\exp (t \boldsymbol A) = \begin{pmatrix} \cos t & \sin t \\ -\sin t& \cos t \end{pmatrix},\]

然后利用$\omega \ne \pm 1$的事实，算出其中的积分值再经过一个长的计算过程后就得到

\[\begin{align} \boldsymbol Y(t) &= \frac{1}{2}\exp (t\boldsymbol A) \begin{pmatrix}\frac{{ - \cos (\omega - 1)t}}{{\omega - 1}} + \frac{{\cos (\omega + 1)t}}{{\omega + 1}}\\\frac{{sin(\omega - 1)t}}{{\omega - 1}} + \frac{{sin(\omega + 1)t}}{{\omega + 1}}\end{pmatrix} + \exp (t\boldsymbol A) \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0\end{pmatrix} \\&= \frac{1}{{{\omega ^2} - 1}} \begin{pmatrix} - \cos \omega t\\\omega \sin \omega t\end{pmatrix} + \exp (t\boldsymbol A) \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0\end{pmatrix} . \end{align}\]

这样就得到了原方程的通解为(因为上式为初值为零时，原方程的通解，加上初值为$X_0$时，相应齐次方程的通解，即得到初值为$X_0$的非齐次方程的通解)

\[\boldsymbol Y(t) = \exp (t\boldsymbol A) \left(\boldsymbol X_0+ \begin{pmatrix}{({\omega ^2} - 1)^{ - 1}}\\0 \end{pmatrix} \right) + \frac{1}{{{\omega ^2} - 1}} \begin{pmatrix} - \cos \omega t\\\omega \sin \omega t\end{pmatrix}. \]

这个表达式中的第一项是以$2\pi$为周期的，而第二项是以$2\pi / \omega$为周期的。与有阻尼的情形不同，现在的解不一定会有周期运动。事实上，这个解为周期的当且仅当$\omega$是一个有理数。当$\omega$为无理数时，这个运动是拟周期的，这正好与6.2节中看到的情形是一样的。

小结：本章主要讨论了高维线性系统解的特性，通过调和振子的例子，引入了一些新的概念，如无理旋转，拟周期。着重讨论了重特征值情形下解的形式和性质，需要结合5.5节中相关内容进行学习，即重特征值情形下如何求特征向量。最后讨论了一个非自治线性系统的例子，用以说明强迫项、稳态解等。