矩阵微分方程

微分方程是变量的未知函数的数学方程，将函数值与不同阶导数联系起来。矩阵微分方程包含多个函数，以向量形式堆叠在一起，并由一个矩阵将它们与导数联系起来。

例如，一阶矩阵常微分方程为

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

其中$\mathbf{x}(t)$是基变量$t$的函数的$n \times 1$向量，$\mathbf{\dot{x}}(t)$是函数的一阶导，$\mathbf{A}(t)$是$n \times n$系数矩阵。

若$\mathbf{A}$为常数且有n个线性无关的特征向量，微分方程有如下一般解：

$$\mathbf{x}(t) = c_1 e^{\lambda_1 t} \mathbf{u}_1 + c_2 e^{\lambda_2 t} \mathbf{u}_2 + \cdots + c_n e^{\lambda_n t} \mathbf{u}_n ~,$$

其中λ1, λ2, …, λn是A的特征值，u1, u2, …, un是A相应的特征向量；c1, c2, …, cn为常数。

更一般地说，若$\mathbf{A}(t)$等于其积分$\int_a^t \mathbf{A}(s)ds$，则马格努斯展开降为前导阶，微分方程的一般解是

$$\mathbf{x}(t)=e^{\int_a^t \mathbf{A}(s) ds} \mathbf{c} ~,$$

其中$\mathbf{c}$是$n \times 1$常向量。

通过使用哈密尔顿–凯莱定理和类范德蒙矩阵，这种形式化的矩阵指数解可简化为一种简单的形式。下面，我们将用普策算法（Putzer's algorithm）来展示这一方法。

## 矩阵系统的稳定性与稳态

矩阵方程

$$\mathbf{\dot{x}}(t) = \mathbf{Ax}(t) + \mathbf{b}$$

当且仅当常数矩阵A的所有特征值的实部都为负，n×1参数常数向量b才稳定。

稳定时收敛到的稳态x*可置

$$\mathbf{\dot{x}}^* (t)=\mathbf{0}~,$$

找到，因此有

$$\mathbf{x}^* = -\mathbf{A}^{-1}\mathbf{b}~,$$

假设A可逆。

因此，原方程可用偏离稳态的齐次形式来写

$$\mathbf{\dot{x}}(t)=\mathbf{A}[\mathbf{x}(t)-\mathbf{x}^*]~.$$

一种等效的表达是，x*是非齐次方程的一个特解，而所有解的形式都是

$$\mathbf{x}_h+\mathbf{x}^* ~,$$

其中$\mathbf{x}_h$是齐次方程（b=0）的解。

### 双状态变量情形的稳定性

n = 2（2个状态变量）时，稳定条件为：过渡矩阵A的两个特征值均有负实部，等价于A的迹为负、行列式为正。

## 矩阵形式的解

$\mathbf{\dot{x}}(t)=\mathbf{A}[\mathbf{x}(t)-\mathbf{x}^*]$的形式解为矩阵指数形式

$$\mathbf{x}(t)=\mathbf{x}^*+e^{\mathbf{A}t}[\mathbf{x}(0)-\mathbf{x}^*] ~,$$

使用多种技术中的任何一种进行评估。

### 计算eAt的普策算法

给定特征值为$\lambda_1,\lambda_2,\dots,\lambda_n$的矩阵A。

$$e^{\mathbf{A}t} = \sum_{j=0}^{n-1}r_{j+1}{\left(t\right)}\mathbf{P}_{j}$$

其中

$$\mathbf{P}_0 = \mathbf{I}$$

$$\mathbf{P}_j = \prod_{k=1}^{j}\left(\mathbf{A}-\lambda_k \mathbf{I}\right)= \mathbf{P}_{j-1} \left(\mathbf{A}-\lambda_j \mathbf{I}\right), \qquad j=1,2,\dots,n-1$$

$$\dot{r}_1 = \lambda_1 r_1$$

$$r_1{\left(0\right)}=1$$

$$\dot{r}_{j} = \lambda_j r_j + r_{j-1}, \qquad j=2,3,\dots,n$$

$$r_j{\left(0\right)}=0, \qquad j=2,3,\dots,n$$

$r_i (t)$的方程是简单的一阶非齐次常微分方程。

注意该算法并不要求矩阵A可对角化，并绕过了通常使用的若尔当标准形的计算。

## 矩阵常微分方程解构示例

一阶齐次矩阵常微分方程包含两个函数x(t)、y(t)，从矩阵形式解出后有如下形式：

$$\frac{dx}{dt}=a_1x+b_1y,\quad\frac{dy}{dt}=a_2x+b_2y$$

其中$a_1$、$a_2$、$b_1$、$b_2$可为任意标量。

高阶矩阵ODE的形式可能复杂得多。

## 解分解后的矩阵常微分方程

求解上述方程并找到这种特定阶次和形式的所需函数的过程大概分3步。每个步骤的简要说明如下：

- 找到特征值
- 找到特征向量
- 找到所需函数

第三部通常是把前两步的结果代入专门形式的一般方程中，下详。

## 矩阵ODE已解示例

要按上述3步解矩阵ODE，并在过程中使用简单矩阵，具体来说，现在下面的一阶齐次线性ODE中找到函数x、函数y，都用单一自变量t表示：

$$\frac{dx}{dt}=3x-4y,\quad\frac{dy}{dt}=4x-7y~.$$

要解这个常微分方程系统，在过程中的某时刻需要一组两个初始条件（对应起点的两个状态变量）。这时先取1=x(0) = y(0) = 1。

### 第一步

第一步即找到A的特征值

$$\begin{bmatrix} x'\\y' \end{bmatrix} = \begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix}\begin{bmatrix} x\\y \end{bmatrix}~.$$

上面的导数记号x′等称为拉格朗日记法（由约瑟夫·拉格朗日提出，等同于前面方程里的dx/dt，这是莱布尼兹记法，得名于戈特弗里德·莱布尼茨）。

一旦两个变量的系数被写为上述矩阵形式A，就可估计特征值了。为此，可求矩阵行列式，即从上述系数矩阵中减去单位矩阵$I_n$乘常数λ，再得到特征多项式

$$\det\left(\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix} - \lambda\begin{bmatrix} 1 & 0\\0 & 1 \end{bmatrix}\right)~,$$

再解得其零点。

进一步简化、应用矩阵加法的基本规则，得出

$$\det\begin{bmatrix} 3-\lambda & -4\\4 & -7-\lambda \end{bmatrix}~.$$

应用求单一2×2矩阵行列式的规则，可得下列一元二次方程

$$\det\begin{bmatrix} 3-\lambda & -4\\4 & -7-\lambda \end{bmatrix} = 0$$

$$-21 - 3\lambda + 7\lambda + \lambda^2 + 16 = 0 \,\!$$

可以进一步简化

$$\lambda^2 + 4\lambda - 5 = 0 ~.$$

应用因式分解得到给定一元二次方程的两个根$\lambda_1$、$\lambda_2$

$$\lambda^2 + 5\lambda - \lambda - 5 = 0$$

$$\lambda (\lambda + 5) - 1 (\lambda + 5) = 0$$

$$(\lambda - 1)(\lambda + 5) = 0$$

$$\lambda = 1, -5 ~.$$

上面算出的$\lambda_1 = 1$、$\lambda_2 = -5$即所求A的特征值。

矩阵ODE的特征值可能是复数，求解过程的下一步及最终形式和解法可能会有巨大变化。

### 第二步

第二步即找到A的特征向量。

对算出的每个特征值，都有单独的特征向量。例如对第一个特征值即$\lambda_1 = 1$，有

$$\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix}\begin{bmatrix} \alpha\\\beta \end{bmatrix} = 1\begin{bmatrix} \alpha\\\beta \end{bmatrix}.$$

应用矩阵乘法规则简化上式，得到

$$3\alpha - 4\beta = \alpha$$

$$\alpha = 2\beta~.$$

所有计算都是为了得到最后一个式子，本例中就是1=α = 2β。现在任取一个无关紧要的小值（这样更容易处理），代入1=α = 2β中的α或β（选哪个并不重要），这样就得到了一个简单的向量，就是这个特定特征值所需的特征向量。在本例中，我们取1=α = 2，得1=β = 1。用标准的向量符号来写，向量是这样的

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

对第二个特征值$\lambda = -5$进行相同的计算，得到第二个特征向量，结果为

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

### 第三步

最后一步是找到“隐藏”在导数背后的所求函数。有两个函数，因为微分方程涉及两个变量。

方程包含之前得到的所有信息，形式如下：

$$\begin{bmatrix} x\\y \end{bmatrix} = Ae^{\lambda_1t}\mathbf{\hat{v}}_1 + Be^{\lambda_2t}\mathbf{\hat{v}}_2.$$

代入特征值和特征向量，得到

$$\begin{bmatrix} x\\y \end{bmatrix} = Ae^{t}\begin{bmatrix} 2\\1 \end{bmatrix} + Be^{-5t}\begin{bmatrix} 1\\2 \end{bmatrix}.$$

简化

$$\begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 2 & 1\\1 & 2 \end{bmatrix}\begin{bmatrix} Ae^{t}\\Be^{-5t} \end{bmatrix}.$$

再简化，分别写出函数x、y的方程

$$x = 2Ae^{t} + Be^{-5t}$$

$$y = Ae^{t} + 2Be^{-5t}.$$

上述方程就是所求的一般函数，但只是一般形式（A、B的值未指定），但我们想找到它们的精确形式和解。因此现在，考虑问题的给定初始条件（即所谓初值问题）。假设给定了$x(0) = y(0) = 1$，是ODE的起点；条件的应用指定了常数A、B。从条件$x(0) = y(0) = 1$可以看出，1=t = 0时，上述方程的左式等于1，由此可构造下列线性方程组

$$1 = 2A + B$$

$$1 = A + 2B~.$$

求解这些等式，发现常数A、B都等于1/3。因此将这些值代入这两个函数的一般形式，就可以得到它们的精确形式$$x = \tfrac{2}{3}e^{t} + \tfrac{1}{3}e^{-5t}$$ $$y = \tfrac{1}{3}e^{t} + \tfrac{2}{3}e^{-5t}~,$$ 所求的两个函数。

### 使用矩阵指数

上述问题可以直接应用矩阵指数法解决。也就是说，可以说

$$\begin{bmatrix} x(t)\\y(t) \end{bmatrix} = \exp \left(\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix} t\right) \begin{bmatrix} x_0(t)\\y_0(t) \end{bmatrix}$$

给出了（可用MATLAB的`expm`工具包之类，或通过对角化，并利用对角矩阵的矩阵指数与元素的指数化相等这一特性来计算）

$$\exp \left(\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix} t\right) = \begin{bmatrix} 4 e^t/3 - e^{-5t}/3 & 2e^{-5t}/3 - 2e^t/3\\2e^t/3 - 2e^{-5t}/3 & 4e^{-5t}/3 - e^t/3 \end{bmatrix}$$

得到最终解

$$\begin{bmatrix} x(t)\\y(t) \end{bmatrix} = \begin{bmatrix} 4 e^t/3 - e^{-5t}/3 & 2e^{-5t}/3 - 2e^t/3\\2e^t/3 - 2e^{-5t}/3 & 4e^{-5t}/3 - e^t/3 \end{bmatrix} \begin{bmatrix} 1\\1 \end{bmatrix}$$

$$\begin{bmatrix} x(t)\\y(t) \end{bmatrix} = \begin{bmatrix} e^{-5t}/3 + 2e^t/3\\ e^t/3 + 2e^{-5t}/3 \end{bmatrix}$$

这与之前展示的特征向量方法相同。