矩陣微分方程

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

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

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

其中λ₁, λ₂, …, λ_n是A的特徵值，u₁, u₂, …, u_n是A相應的特徵向量；c₁, c₂, …, c_n為常數。

更一般地說，若$\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} ^{*}]~,$

使用多種技術中的任何一種進行評估。

計算e^At的普策算法

給定特徵值為$\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_{1}x+b_{1}y,\quad {\frac {dy}{dt}}=a_{2}x+b_{2}y$

其中$a_{1}$、$a_{2}$、$b_{1}$、$b_{2}$可為任意標量。 高階矩陣ODE的形式可能複雜得多。

解分解後的矩陣常微分方程

求解上述方程並找到這種特定階次和形式的所需函數的過程大概分3步。每個步驟的簡要說明如下：

• 找到特徵值
• 找到特徵向量
• 找到所需函數

第三部通常是把前兩步的結果代入專門形式的一般方程中，下詳。

矩陣ODE已解示例

參見：矩陣指數 § 線性微分方程_2

要按上述3步解矩陣ODE，並在過程中使用簡單矩陣，具體來說，現在下面的一階齊次線性ODE中找到函數x、函數y，都用單一自變量t表示：

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

要解這個常微分方程系統，在過程中的某時刻需要一組兩個初始條件（對應起點的兩個狀態變量）。這時先取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 ~.$

所有計算都是為了得到最後一個式子，本例中就是α = 2β。現在任取一個無關緊要的小值（這樣更容易處理），代入α = 2β中的α或β（選哪個並不重要），這樣就得到了一個簡單的向量，就是這個特定特徵值所需的特徵向量。在本例中，我們取α = 2，得β = 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 _{1}t}\mathbf {\hat {v}} _{1}+Be^{\lambda _{2}t}\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$可以看出，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}4e^{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}4e^{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}}$

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

另見

• 齊次微分方程
• 矩陣差分方程
• 冷卻定律
• 斐波那契數列
• 差分方程
• 波動方程
• 自治系統 (數學)