线性微分方程组

本篇主要对应《工科数学分析》中关于线性微分方程组内容的证明和解释。《工科数学分析》只讲了计算的办法，但是证明却都是略，造成了大量没头没脑的公式。

下面我们考虑求解线性微分方程组：

\[ \boldsymbol x'(t) = \boldsymbol{Ax}(t) + \boldsymbol f(t) \]

前置：Jordan 标准型

Jordan 块

Jordan 块是如下形式的方阵：

\[ \boldsymbol J = \begin{bmatrix} \lambda & 1 \\ & \lambda & 1 \\ & & \dots & \dots\\ & & & \lambda & 1 \\ & & & & \lambda \end{bmatrix} \]

即主对角线为相等数，副对角线全为 1。特别地，一阶 Jordan 标准型只有元素 \(\lambda\)。

Jordan 标准型

Jordan 标准型是 Jordan 块构成的对角阵。

\[ \boldsymbol J = \text{diag}(\boldsymbol J_{1},\boldsymbol J_{2},\dots,\boldsymbol J_{m}) \]

可以证明，任何实方阵都相似于一个 Jordan 标准型。即存在可逆矩阵 \(\boldsymbol P\)，使得 \(\boldsymbol A = \boldsymbol {PJP}^{-1}\)。

特别地，如果 \(\boldsymbol A\) 可对角化，则 \(\boldsymbol J\) 是由一阶 Jordan 块构成的。

\(\boldsymbol J\) 的主对角线就是 \(A\) 的特征值。而 \(\boldsymbol P\) 的列向量是所谓广义特征向量。如果特征值 \(\lambda_{i}\) 的代数重数是 \(c_{i}\)，则其对应的广义特征向量 \(u\) 是如下方程的非零解：

\[ (\boldsymbol A - \lambda_{i} \boldsymbol E)^{c_{i}}\boldsymbol u = \boldsymbol 0 \]

可以证明总有 \(c_{i}\) 个线性无关的解，所以总可以找到 \(n\) 个线性无关的广义特征向量，进而可以构成 \(\boldsymbol P\)。特别地，一个特征值的广义特征向量至少有一个是特征向量。

更多性质参考 特征值分解与 Jordan 标准型。只要记住 \(\boldsymbol A = \boldsymbol {PJP} ^{-1}\)。

前置：矩阵指数

定义

考虑 \(\exp(x) = e^x\) 的级数形式：

\[ \exp(x) = \sum_{n=0}^{\infty} \frac{x^{n}}{n!} \]

只要 \(x\) 可加，可求幂，我们就可以在形式上写出 \(\exp(x)\)。

不难想到方阵满足这些条件。所以我们定义，对于方阵 \(\boldsymbol A\)，有：

\[ \exp(\boldsymbol A) = \sum_{n=0}^{\infty} \frac{\boldsymbol A^n}{n!} \]

规定任意矩阵的零次幂为 \(\boldsymbol E\)。

根据这个形式，我们知道如下性质：

• \(\exp(\boldsymbol O)=\boldsymbol E\)。
• \(\exp(\boldsymbol {PAP}^{-1}) = \boldsymbol P\exp(\boldsymbol A)\boldsymbol P^{-1}\)。这可以由分配律证得。
• 如果 \(\boldsymbol {XY}=\boldsymbol{YX}\)，则有 \(\exp(\boldsymbol X + \boldsymbol Y) = \exp(\boldsymbol X)\exp(\boldsymbol Y)\)。这是因为当 \(\boldsymbol{XY}=\boldsymbol{YX}\) 成立时，二项式定理就可以用来计算 \((\boldsymbol X + \boldsymbol Y)^n\)。将级数展开得到：

\[ \begin{aligned} \exp(\boldsymbol X)\exp(\boldsymbol Y) &= \sum_{n=0}^\infty \sum_{m=0}^\infty \frac{\boldsymbol X^n \boldsymbol Y^m}{n!m!}\\ &= \sum_{k=0}^{\infty}\sum_{n=0}^k \frac{\boldsymbol X^n\boldsymbol Y^{k-n}}{n!(k-n)!}\\ &= \sum_{k=0}^\infty \sum_{n=0}^k \frac{k!}{n!(k-n)!}\cdot \frac{\boldsymbol X^n \boldsymbol Y^{k-n}}{k!}\\ &= \sum_{k=0}^\infty \frac{1}{k!}\sum_{n=0}^k \binom{k}{n}\boldsymbol X^n \boldsymbol Y^{k-n}\\ &= \sum_{k=0}^\infty \frac{(\boldsymbol X + \boldsymbol Y)^k}{k!}\\ &= \exp(\boldsymbol X+\boldsymbol Y) \end{aligned} \]

• 从上一条不难得到，因为 \((-\boldsymbol X)\boldsymbol X =\boldsymbol X(-\boldsymbol X)\)，则 \(\boldsymbol E = \exp(\boldsymbol O) = \exp(\boldsymbol X + (- \boldsymbol X)) = \exp(\boldsymbol X)\exp(-\boldsymbol X)\)，所以 \(\exp^{-1}(\boldsymbol X) = \exp(-\boldsymbol X)\)。

计算

如何把这个无穷级数变成收敛的形式呢？如果 \(\boldsymbol A\) 可对角化，即 \(\boldsymbol A = \boldsymbol {PDP}^{-1} = \boldsymbol P\text{diag}(\lambda_{1},\dots,\lambda_{n}) \boldsymbol P^{-1}\)，则我们有：

\[ \begin{aligned} \exp(\boldsymbol A) &= \boldsymbol P \exp(\boldsymbol D)\boldsymbol P^{-1}\\ &=\boldsymbol P\sum_{i=0}^{\infty} \text{diag}\left( \frac{\lambda_{1}^i}{i!},\dots, \frac{\lambda_{n}^i}{i!} \right) \boldsymbol P^{-1}\\ &=\boldsymbol P \text{diag}(e^{\lambda_{1}},\dots,e^{\lambda_{n}})\boldsymbol P^{-1} \end{aligned} \]

而如果 \(\boldsymbol A\) 只是可以 Jordan 标准化，即 \(\boldsymbol A = \boldsymbol {PJP}^{-1}\)，同样有：

\[ \begin{aligned} \exp(\boldsymbol A) &= \boldsymbol P \exp(\boldsymbol J) \boldsymbol P^{-1}\\ &= \boldsymbol P \text{diag}(\exp(\boldsymbol J_{1}),\dots,\exp(\boldsymbol J_{m}))\boldsymbol P^{-1} \end{aligned} \]

对于 Jordan 块 \(\boldsymbol J_{i}\)，它的 \(\exp\) 有通式：

\[ \exp(\boldsymbol J_{i}) = \begin{bmatrix} \lambda_{i} & \frac{\lambda_{i}}{1!} & \frac{\lambda_{i}}{2!} & \dots & \frac{\lambda_{i}^{c_{i}-1}}{(c_{i}-1)!} \\ & \lambda _{i} & \frac{\lambda_{i}}{1!} & \dots & \frac{\lambda_{i}^{c_{i} - 2}}{(c_{i} - 2)!} \\ & & \dots & \dots & \dots \\ & & & \lambda_{i} & \frac{\lambda_{i}}{1!} \\ & & & & \lambda_{i} \end{bmatrix} \]

其中 \(c_{i}\) 是它的阶数。所以我们可以算出任何方阵的 \(\exp\)。

自变量与求导

特别地，如果是带自变量的矩阵 \(\boldsymbol At\)，则：

\[ \begin{aligned} \exp(\boldsymbol At) &= \sum_{n=0}^{\infty} \frac{\boldsymbol A^n}{n!}t^n\\ &= \boldsymbol P \exp(\boldsymbol Jt) \boldsymbol P^{-1} \end{aligned} \]

如果对 \(t\) 求导，我们得到如下结果：

\[ \exp'(\boldsymbol At) = \sum_{n=1}^{\infty} \frac{\boldsymbol A^n}{(n-1)!}t^{n-1} = \boldsymbol A\sum_{n=0}^{\infty} \frac{\boldsymbol A^n}{n!}t^n = \boldsymbol A\exp(\boldsymbol At) \]

齐次方程组

《工科数学分析》中首先分析了齐次方程组的情况。即：

\[ \boldsymbol x'(t) = \boldsymbol {Ax}(t) \]

基解矩阵

定义

首先要明白，\(n\) 个变量的线性微分方程组总有 \(n\) 个线性无关的特解。因为将矩阵展开后会有 \(n\) 个形如 \(x_{i}' = f(\boldsymbol x)\) 的形式，而每个等式在积分之后都会带一个常数。这些常数是互相独立的，所以解空间总维度是 \(n\) 维的。我们将这样 \(n\) 个线性无关的特解叫做基解，因为它们是解空间的基。

这样，任何一个解都可以表示成基解的线性组合。如果将 \(n\) 个基解作为列向量组成一个矩阵 \(\boldsymbol X\)，则所有解都可以表示成 \(\boldsymbol{Xc}\) 的形式，其中 \(\boldsymbol c\) 是一个 \(\mathbb{R}^n\) 上的常数向量。我们把基解构成的矩阵称作基解矩阵。

性质

注意到基解矩阵满足 \(\boldsymbol X' = \boldsymbol {AX}\)。（每个列向量都满足）

不难知道，只要满足这个式子的满秩矩阵都是基解矩阵。

由于 \((\boldsymbol {XP})' = \boldsymbol X' \boldsymbol P\)，所以基解矩阵右乘一个常数可逆阵还是基解矩阵。

求解

指数解

求微分方程组就是求基解矩阵。

从我们的前置知识可以知道，\(\exp(\boldsymbol At)' = \boldsymbol A\exp (\boldsymbol At)\)，则 \(\exp(\boldsymbol At) = \boldsymbol P \exp(\boldsymbol Jt) \boldsymbol P^{-1}\) 是基解矩阵。

具体计算

《工科数学分析》中对 \(\boldsymbol A\) 的情况进行了分类讨论，分成了可对角化和不可对角化的情况。

可对角化

在可对角化时，记 \(\boldsymbol A = \boldsymbol{PDP}^{-1}\)，书上说明了特征值 \(\lambda_{i}\) 和特征向量 \(\boldsymbol u_{i}\) 可以构成一个特解 \(\boldsymbol u_{i}e^{\lambda_{i}t}\)。

这是因为：

\[ (\boldsymbol u_{i}e^{\lambda_{i}t})' = \lambda_{i}\boldsymbol u_{i}e^{\lambda_{i}t} = \boldsymbol {Au}_{i}e^{\lambda t} \]

同时注意到 \(\boldsymbol u_{i}\) 线性无关。那么此时，基解矩阵是 \(\boldsymbol u_{i}e^{\lambda_{i}t}\) 构成的矩阵，即 \(\boldsymbol P \exp(\boldsymbol Dt)\)。

和我们求出的 \(\exp(\boldsymbol At)\) 对比，其实就是前面说的性质：基解矩阵右乘可逆阵还是基解矩阵。《工科数学分析》上给出的答案，其实是 \(\exp (\boldsymbol At) \boldsymbol P\)。用 \(\exp\) 固然是简洁的通解，但是在具体计算，特别是维度较小时，采用书上的办法不需要求 \(\boldsymbol P^{-1}\)，还更快捷。

不可对角化

不可对角化的情况，使用 \(\exp\) 也是一样的，通解矩阵是 \(\exp(\boldsymbol At) = \boldsymbol P \exp (\boldsymbol Jt) \boldsymbol P^{-1}\)。而《工科数学分析》则不加证明地给出 3.18：对于重数是 \(c_{i}\) 的特征值 \(\lambda_{i}\)，总有 \(c_{i}\) 个如下形式的线性无关的特解：

\[ \boldsymbol x_{i} = e^{\lambda_{i} t}\left(\sum_{n=0}^{c_{i}-1} \frac{\boldsymbol r_{n}}{n!} \right) \]

其中 \(\boldsymbol r_{0}\) 是 \((\boldsymbol A - \lambda_{i} \boldsymbol E)^{c_{i}}\boldsymbol x=0\) 的解，而 \(\boldsymbol r_{n+1} = (\boldsymbol A - \lambda_{i}\boldsymbol E)\boldsymbol r_{n}\)。

如果仔细学习了 Jordan 标准型的知识，就知道 \(\boldsymbol r_{i}\) 都是所谓广义特征向量，而 \(\boldsymbol r_{n+1} = (\boldsymbol A - \lambda_{i}\boldsymbol E)\boldsymbol r_{n}\) 是广义特征向量的性质，来自所谓的循环子空间。而 \(\frac{1}{n!}\) 的阶乘的系数，其实是 Jordan 块 \(\exp(\boldsymbol J_{i})\) 的系数。

总之，和齐次的情况相同，《工科数学分析》不加证明地给出了基解矩阵 \(\exp(\boldsymbol At)\boldsymbol P = \boldsymbol P \exp(\boldsymbol Jt)\) 的计算公式，在具体计算，特别是维度较小时，采用书上的办法不需要求 \(\boldsymbol P^{-1},\boldsymbol J\)，还更快捷。

非齐次方程组

对于非齐次方程组 \(\boldsymbol x' = \boldsymbol {Ax} + \boldsymbol f\)，我们先求出齐次方程组的基解矩阵 \(\boldsymbol X = \exp(\boldsymbol At)\)。则齐次的通解为 \(\boldsymbol {Xc}\)。

参考一元的常数变易法，我们假设非齐次的通解是 \(\boldsymbol {XYc}\)，其中 \(\boldsymbol Y = \boldsymbol Y(t)\) 是 \(t\) 的函数矩阵。则有：

\[ \begin{aligned} &&(\boldsymbol{XYc})' &= \boldsymbol {AXYc} + \boldsymbol f \\ &\implies&\boldsymbol X'\boldsymbol {Yc} + \boldsymbol X (\boldsymbol {Yc})' &= \boldsymbol {AXYc} + \boldsymbol f\\ &\implies& \boldsymbol X (\boldsymbol {Yc})' &= \boldsymbol f\\ &\implies&\boldsymbol {Yc} &= \int \boldsymbol X^{-1}(t)\boldsymbol f(t) \mathrm{d} t\\ \end{aligned} \]

其中第二步运用了齐次解的性质。

则我们最终的答案是：

\[ \begin{aligned} \boldsymbol{XYc} &= \boldsymbol X(t) \int \boldsymbol X^{-1}(t) \boldsymbol f(t) \mathrm{d}t&(1)&\\ &= \exp(\boldsymbol At)\int \exp(-\boldsymbol At)\boldsymbol f(t) \mathrm{d}t &(2)&\\ &= \boldsymbol P \exp(\boldsymbol Jt)\boldsymbol P^{-1} \int \boldsymbol P \exp(-\boldsymbol Jt)\boldsymbol P^{-1} \boldsymbol f(t)\mathrm{dt} &&\\ &= \boldsymbol P \exp(\boldsymbol Jt) \int \exp(-\boldsymbol Jt)\boldsymbol P^{-1} \boldsymbol f(t) \mathrm{dt} &(3)&\\ \end{aligned} \]

其中的 \((1)\) 是对于无论什么方法求出的基解矩阵都成立，对应《工科数学分析》的 3.12。注意这里是不定积分，而 3.12 是定积分，所以 3.12 多一个 \(\boldsymbol {Xc}\) 项。

而 \((2)\) 是我们采用的形式化的比较无脑的指数解。\((3)\) 则是带入《工科数学分析》中方便计算的基解矩阵 \(\exp (\boldsymbol At) \boldsymbol P=\boldsymbol P \exp(\boldsymbol Jt)\)。

这个时候，\((3)\) 的积分项内不得不计算 \(\boldsymbol P^{-1}\)，原先的方便之处不存在了。《工科数学分析》在定理 3.9 中提出了这一点。

但是 \(\boldsymbol X(0)=\boldsymbol E\) 这个条件令人莫名其妙。实际上，\(\boldsymbol X = \exp (\boldsymbol At)\) 正好满足这一点，因为 \(\exp^{-1}(\boldsymbol At) = \exp(-\boldsymbol At),\exp(\boldsymbol O) = \boldsymbol E\)。

《工科数学分析》想要表达的意思是，原来的 \(\exp(\boldsymbol At)\boldsymbol P\) 在这里不好用了，要右乘一个 \(\boldsymbol P^{-1}\) 变成 \(\exp (\boldsymbol At)\)。但是又没办法直接说明，才用了 \(\boldsymbol X(0) = \boldsymbol E\) 的要求。事实上，接下来的例题也是这么做的：

它所谓的 \(\boldsymbol X\) 其实是 \(\exp (\boldsymbol At) \boldsymbol P\)，而 \(\boldsymbol X^{-1}(0)\) 就是 \(\boldsymbol P^{-1}\)。

特解

最后，作为结尾，给出在知道 \(\boldsymbol x(t_{0})=\boldsymbol x_{0}\) 时的特解。

齐次方程组的特解为：

\[ \boldsymbol x(t) = \exp (\boldsymbol A(t - t_{0})) \boldsymbol x_{0} \]

而非齐次方程组的特解为：

\[ \begin{aligned} \boldsymbol X(t) &= \exp(\boldsymbol At)\\ \boldsymbol x(t) &= \boldsymbol X(t - t_{0}) \boldsymbol x_{0} + \int_{t_{0}}^t \boldsymbol X(t - u)\boldsymbol f(u)\mathrm{du} \end{aligned} \]