矩陣微分方程

微分方程是變量的未知函數的數學方程，將函數值與不同階導數聯繫起來。矩陣微分方程包含多個函數，以向量形式堆疊在一起，並由一個矩陣將它們與導數聯繫起來。

例如，一階矩陣常微分方程為

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

這與之前展示的特徵向量方法相同。