一阶线性系统与矩阵方法

前面几章里，一个未知函数通常代表一个量：温度、位移、电流、药物浓度。很多模型却不能只看一个量。一个水箱里的盐量会受另一个水箱影响，弹簧振子的位移和速度互相决定，电路中多个电流或电压也会彼此耦合。

这时我们把若干未知函数放进一个向量。矩阵的作用不再只是线性代数里的“变换一个向量”，而是变成一个动态规则：现在的状态在哪里，下一瞬间的速度就由矩阵算出来。

多个变量的变化规律可以统一写成状态向量和矩阵方程。

本章只讨论最核心的一类系统：

$$\mathbf{x}'(t)=A\mathbf{x}(t)$$

其中 $\mathbf{x}(t)$ 是状态向量，$A$ 是常数矩阵。我们先把这种写法读懂，再用特征值、特征向量、矩阵指数和基本矩阵来组织解。

---

从一个方程到一个系统

一个 $n$ 维一阶系统由 $n$ 个未知函数组成：

$$\begin{aligned} x_1' &= a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n,\\ x_2' &= a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n,\\ &\ \vdots\\ x_n' &= a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n. \end{aligned}$$

把未知函数写成状态向量

$$\mathbf{x}(t)= \begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots\\ x_n(t) \end{bmatrix}$$

再把系数写成矩阵

$$A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix},$$

整个系统就压缩成

$$\mathbf{x}'=A\mathbf{x}.$$

这个写法很短，但信息没有减少。第 $i$ 行描述的是 $x_i$ 的变化率怎样由所有状态变量共同决定。非对角元素表示耦合：例如 $a_{12}\neq 0$ 说明 $x_2$ 会直接影响 $x_1'$。

在每个状态点，矩阵 $A$ 把当前位置变成速度向量。

---

把高阶方程改写成状态方程

一阶系统不只用来描述“本来就有多个量”的模型。任何高阶微分方程也可以改写成一阶系统，这就是很多数值算法和系统理论偏爱一阶系统的原因。

以二阶方程为例：

$$y''+ay'+by=0.$$

令

$$x_1=y,\qquad x_2=y'.$$

那么

$$x_1'=y'=x_2,$$

而

$$x_2'=y''=-by-ay'.$$

用 $x_1$ 和 $x_2$ 表示就是

$$\begin{aligned} x_1' &= x_2,\\ x_2' &= -bx_1-ax_2. \end{aligned}$$

写成矩阵形式：

$$\begin{bmatrix} x_1'\\ x_2' \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -b & -a \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix}.$$

这里的状态不是只记录位置 $y$，还记录速度 $y'$。知道位置但不知道速度，通常不能决定下一步会怎样走。

---

矩阵怎样给出几何运动

在二维系统中，状态向量可以看成平面上的点：

$$\mathbf{x}(t)= \begin{bmatrix} x(t)\\ y(t) \end{bmatrix}.$$

方程 $\mathbf{x}'=A\mathbf{x}$ 的意思是：点在当前位置 $\mathbf{x}$ 时，它的速度向量就是 $A\mathbf{x}$。把平面上许多点都代入 $A\mathbf{x}$，就得到一个速度场。

例如

$$A= \begin{bmatrix} 2 & 1\\ 0 & -1 \end{bmatrix}.$$

在点 $(1,2)$ 处，

$$A \begin{bmatrix} 1\\ 2 \end{bmatrix} = \begin{bmatrix} 4\\ -2 \end{bmatrix}.$$

所以经过 $(1,2)$ 的解曲线，在那一瞬间向右且向下运动。换一个点，速度向量也换。解曲线就是始终沿着这些速度箭头走出来的轨迹。

矩阵的每一列也有可解释的意义。第一列告诉我们“只改变 $x_1$ 分量时，速度会增加多少”；第二列告诉我们“只改变 $x_2$ 分量时，速度会增加多少”。矩阵乘法把这些贡献加起来。

---

特征向量给出不会转弯的方向

解 $\mathbf{x}'=A\mathbf{x}$ 时，最有用的方向是特征向量方向。若

$$A\mathbf{v}=\lambda\mathbf{v},$$

则 $\mathbf{v}$ 是特征向量，$\lambda$ 是对应特征值。假设解沿着 $\mathbf{v}$ 的方向：

$$\mathbf{x}(t)=u(t)\mathbf{v}.$$

代入系统：

$$u'(t)\mathbf{v}=A(u(t)\mathbf{v})=u(t)A\mathbf{v}=u(t)\lambda\mathbf{v}.$$

只要 $\mathbf{v}\neq \mathbf{0}$，这就变成标量方程

$$u'=\lambda u.$$

于是

$$u(t)=Ce^{\lambda t},$$

所以得到一个向量解：

$$\mathbf{x}(t)=Ce^{\lambda t}\mathbf{v}.$$

沿特征方向出发，方向保持不变，只是按指数因子伸缩。

如果 $\lambda>0$，沿这个方向的分量会指数增长；如果 $\lambda<0$，它会指数衰减；如果 $\lambda=0$，这个特征方向上的分量保持常数。这里的“增长”和“衰减”说的是特征方向上的模式，不一定等于每个坐标都单调变化。

---

例题：由特征值写出通解

求解系统

$$\mathbf{x}'= \begin{bmatrix} 3 & 1\\ 0 & 1 \end{bmatrix} \mathbf{x}.$$

这里

$$A= \begin{bmatrix} 3 & 1\\ 0 & 1 \end{bmatrix}.$$

如果初始条件是

$$\mathbf{x}(0)= \begin{bmatrix} 4\\ -2 \end{bmatrix},$$

则

$$c_1 \begin{bmatrix} 1\\ 0 \end{bmatrix} +c_2 \begin{bmatrix} 1\\ -2 \end{bmatrix} = \begin{bmatrix} 4\\ -2 \end{bmatrix}.$$

由第二行得 $c_2=1$，再由第一行得 $c_1=3$。因此特解为

$$\mathbf{x}(t)= 3e^{3t} \begin{bmatrix} 1\\ 0 \end{bmatrix} +e^t \begin{bmatrix} 1\\ -2 \end{bmatrix}.$$

可对角化时，先把初值分解到特征方向，再让每个模式独立演化。

---

可对角化情形的矩阵写法

如果 $A$ 有 $n$ 个线性无关特征向量

$$\mathbf{v}_1,\mathbf{v}_2,\ldots,\mathbf{v}_n,$$

对应特征值

$$\lambda_1,\lambda_2,\ldots,\lambda_n,$$

则通解可以写成

$$\mathbf{x}(t)=c_1e^{\lambda_1t}\mathbf{v}_1+c_2e^{\lambda_2t}\mathbf{v}_2+\cdots+c_ne^{\lambda_nt}\mathbf{v}_n.$$

把特征向量作为列放进矩阵

$$P= \begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_n \end{bmatrix},$$

把特征值放进对角矩阵

$$D= \begin{bmatrix} \lambda_1 & 0 & \cdots & 0\\ 0 & \lambda_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \lambda_n \end{bmatrix},$$

就有

$$A=PDP^{-1}.$$

这说明在特征向量坐标中，系统被拆成了互不耦合的标量方程：

$$\mathbf{z}'=D\mathbf{z}.$$

也就是

$$z_i'=\lambda_i z_i.$$

原坐标里的耦合，换到特征向量坐标后变成几个独立指数模式。这个视角比单纯套公式更重要。

---

矩阵指数与状态推进

一阶标量方程

$$x'=ax$$

的解是

$$x(t)=e^{at}x(0).$$

矩阵系统的对应形式是

$$\mathbf{x}(t)=e^{At}\mathbf{x}(0).$$

这里的 $e^{At}$ 称为矩阵指数，它用幂级数定义：

$$e^{At}=I+At+\frac{A^2t^2}{2!}+\frac{A^3t^3}{3!}+\cdots.$$

这个定义看起来抽象，但它有一个直接意义：$e^{At}$ 是把初始状态推进到时间 $t$ 的状态转移矩阵。

矩阵指数把初始状态沿系统轨线推进到时刻 $t$。

如果 $A=PDP^{-1}$，那么

$$e^{At}=Pe^{Dt}P^{-1},$$

其中

$$e^{Dt}= \begin{bmatrix} e^{\lambda_1t} & 0 & \cdots & 0\\ 0 & e^{\lambda_2t} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & e^{\lambda_nt} \end{bmatrix}.$$

这和前面“分解初值、各模式独立演化、再相加”的解释是同一件事。

---

基本矩阵把所有解装进一个矩阵

对于系统

$$\mathbf{x}'=A\mathbf{x},$$

若 $\mathbf{x}_1(t),\mathbf{x}_2(t),\ldots,\mathbf{x}_n(t)$ 是 $n$ 个线性无关的向量解，把它们作为列排成矩阵：

$$\Phi(t)= \begin{bmatrix} \mathbf{x}_1(t) & \mathbf{x}_2(t) & \cdots & \mathbf{x}_n(t) \end{bmatrix}.$$

这个矩阵称为一个基本矩阵。因为每一列都满足系统，所以

$$\Phi'(t)=A\Phi(t).$$

只要 $\Phi(t)$ 的列线性无关，任意解都可以写成

$$\mathbf{x}(t)=\Phi(t)\mathbf{c},$$

其中 $\mathbf{c}$ 是常数向量。

基本矩阵的每一列都是一个独立解，乘常数向量后得到任意解。

矩阵指数本身就是一个特殊的基本矩阵，因为

$$e^{At}$$

的每一列可以看成从标准基向量出发的解。它还满足

$$e^{A\cdot 0}=I.$$

因此，当初始时刻取 $0$ 时，

$$\mathbf{x}(t)=e^{At}\mathbf{x}(0)$$

是初值问题最紧凑的写法。

---

建模例子：两个连通水箱

设两个水箱中的盐量分别为 $x_1(t)$ 和 $x_2(t)$，单位是千克。每个水箱体积保持不变，水流使盐在两个水箱之间交换，并且第二个水箱还有排出流。假设混合均匀，模型简化为

$$\begin{aligned} x_1' &= -0.3x_1+0.1x_2,\\ x_2' &= 0.3x_1-0.4x_2. \end{aligned}$$

矩阵形式是

$$\mathbf{x}'= \begin{bmatrix} -0.3 & 0.1\\ 0.3 & -0.4 \end{bmatrix} \mathbf{x}.$$

这个矩阵里，$-0.3$ 表示第一个水箱自己的盐量会因流出而降低，$0.1$ 表示第二个水箱的盐会流入第一个水箱。第二行同理。行内各项不是抽象符号，而是“流入减流出”的账本。

求特征值：

$$\det(A-\lambda I)= \begin{vmatrix} -0.3-\lambda & 0.1\\ 0.3 & -0.4-\lambda \end{vmatrix} =0.$$

展开得

$$(\lambda+0.3)(\lambda+0.4)-0.03=0,$$

也就是

$$\lambda^2+0.7\lambda+0.09=0.$$

两个特征值是

$$\lambda_1=\frac{-0.7+\sqrt{0.13}}{2},\qquad \lambda_2=\frac{-0.7-\sqrt{0.13}}{2}.$$

它们都为负，所以两个特征模式都会衰减。模型预测：如果没有新的盐输入，两个水箱中的盐量最终都会趋向 $0$。衰减较慢的特征值对应长期占主导的模式。

---

常见误区与自查

把状态变量选错

把二阶方程改写成系统时，常见错误是令 $x_2$ 等于另一个看起来方便的量，却没有保证 $x_1'=x_2$。标准选择是按导数链建立状态：

$$x_1=y,\quad x_2=y',\quad x_3=y'',\quad \ldots$$

这样每个新变量的导数都有明确来源。

把特征值当成完整解

特征值只给出指数倍率，特征向量给出方向。对系统来说，$e^{\lambda t}$ 只是模式的时间因子，完整模式是

$$e^{\lambda t}\mathbf{v}.$$

忽略初始条件的分解

通解写出后，初始条件不是随便代一个常数，而是要把 $\mathbf{x}(0)$ 表示成特征向量的线性组合。系数 $c_i$ 来自这个分解。

逐项计算矩阵指数

如果

$$A= \begin{bmatrix} 0 & 1\\ -2 & -3 \end{bmatrix},$$

那么 $e^{At}$ 不是

$$\begin{bmatrix} e^0 & e^t\\ e^{-2t} & e^{-3t} \end{bmatrix}.$$

矩阵指数由矩阵幂级数或对角化等方法计算，矩阵乘法中的耦合会进入结果。

---

练习

练习一

把方程

$$y''+4y'+3y=0$$

改写成一阶矩阵系统。

显示答案

令 $x_1=y$，$x_2=y'$。则

$$x_1'=x_2,\qquad x_2'=y''=-4y'-3y=-3x_1-4x_2.$$

所以

$$\begin{bmatrix} x_1'\\ x_2' \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -3 & -4 \end{bmatrix} \begin{bmatrix} x_1\\ x_2 \end{bmatrix}.$$

练习二

求系统

$$\mathbf{x}'= \begin{bmatrix} 2 & 0\\ 0 & -1 \end{bmatrix} \mathbf{x},\qquad \mathbf{x}(0)= \begin{bmatrix} 3\\ 4 \end{bmatrix}$$

的解。

显示答案

这个矩阵已经是对角矩阵，两个坐标互不耦合：

$$x_1'=2x_1,\qquad x_2'=-x_2.$$

由初值可得

$$x_1(t)=3e^{2t},\qquad x_2(t)=4e^{-t}.$$

所以

$$\mathbf{x}(t)= \begin{bmatrix} 3e^{2t}\\ 4e^{-t} \end{bmatrix}.$$

练习三

设

$$A= \begin{bmatrix} 1 & 1\\ 0 & 2 \end{bmatrix}.$$

求 $A$ 的特征值和一组对应特征向量，并写出 $\mathbf{x}'=A\mathbf{x}$ 的通解。

显示答案

因为 $A$ 是上三角矩阵，特征值为 $\lambda_1=1$ 和 $\lambda_2=2$。

对 $\lambda_1=1$，

$$A-I= \begin{bmatrix} 0 & 1\\ 0 & 1 \end{bmatrix},$$

所以可取

$$\mathbf{v}_1= \begin{bmatrix} 1\\ 0 \end{bmatrix}.$$

对 $\lambda_2=2$，

$$A-2I= \begin{bmatrix} -1 & 1\\ 0 & 0 \end{bmatrix},$$

所以可取

$$\mathbf{v}_2= \begin{bmatrix} 1\\ 1 \end{bmatrix}.$$

通解为

$$\mathbf{x}(t)= c_1e^t \begin{bmatrix} 1\\ 0 \end{bmatrix} +c_2e^{2t} \begin{bmatrix} 1\\ 1 \end{bmatrix}.$$

练习四

判断下面说法是否正确：若 $A$ 是任意 $2\times 2$ 矩阵，则 $e^{At}$ 可以通过对 $A$ 的每个元素分别求指数得到。

显示答案

这个说法不正确。矩阵指数由

$$e^{At}=I+At+\frac{A^2t^2}{2!}+\frac{A^3t^3}{3!}+\cdots$$

定义，里面出现的是矩阵幂 $A^2,A^3,\ldots$。这些矩阵幂会把不同元素通过矩阵乘法混合起来，所以一般不能逐项求指数。只有在非常特殊的情形下，例如对角矩阵，结果才表现为对角线元素分别取指数。