一阶线性系统与矩阵方法

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

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

多变量模型汇聚成矩阵方程的示意图

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

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

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

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

本章承接线性代数 I，但重点不是重新复习矩阵运算。真正要看清的是：矩阵 $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$ 说明 $_{}$ 会直接影响。

二维线性系统的速度场示意图

在每个状态点，矩阵 $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$ 分量时，速度会增加多少”。矩阵乘法把这些贡献加起来。

常数矩阵系统是自治系统：右边不显含 $t$ 。所以同一个状态点永远给出同一个速度向量。解曲线不能在同一点交叉后走向不同未来。

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

解 $\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}.

先求特征值。因为 $A$ 是上三角矩阵，特征值就是对角线上的数： $\lambda_1=3$ ， $\lambda_2=1$ 。

如果初始条件是

\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}.

两个特征模式线性组合成系统解的示意图

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

常见错误是只写出 $e^{\lambda t}$ ，忘记乘对应的特征向量。系统的解是向量函数，不能只给一个标量指数函数。

可对角化情形的矩阵写法

如果 $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}.

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

矩阵指数不是把矩阵中每个元素逐项求指数。一般也不能把 $e^{(A+B)t}$ 拆成 $e^{At}e^{Bt}$ ，除非相关矩阵满足可交换条件。矩阵乘法的顺序会改变结果。

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

对于系统

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

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

\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$ 。衰减较慢的特征值对应长期占主导的模式。

线性系统的特征值不仅用于求公式，还能直接解释长期趋势。负实部意味着对应模式衰减；正实部意味着对应模式增长；接近 $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'$ 。则

练习二

求系统

\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.

练习三

设

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$ 是任意 $2\times 2$ 矩阵，则 $e^{At}$ 可以通过对 $A$ 的每个元素分别求指数得到。

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

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

x

2

x_2

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

常微分方程 I：一阶方程、线性系统与建模｜第 16 章：一阶线性系统与矩阵方法 | 自在学