微分方程是變量的未知函數的數學方程,將函數值與不同階導數聯繫起來。矩陣微分方程包含多個函數,以向量形式堆疊在一起,並由一個矩陣將它們與導數聯繫起來。
例如,一階矩陣常微分方程為
![{\displaystyle \mathbf {\dot {x}} (t)=\mathbf {A} (t)\mathbf {x} (t)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/20ddfdfda2aad19899c85fe38e74151c88cb797b)
其中
是基變量
的函數的
向量,
是函數的一階導,
是
係數矩陣。
若
為常數且有n個線性無關的特徵向量,微分方程有如下一般解:
![{\displaystyle \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}~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cb0d853eb3595fdc98efff1b59071506b8419008)
其中λ1, λ2, …, λn是A的特徵值,u1, u2, …, un是A相應的特徵向量;c1, c2, …, cn為常數。
更一般地說,若
等於其積分
,則馬格努斯展開降為前導階,微分方程的一般解是
![{\displaystyle \mathbf {x} (t)=e^{\int _{a}^{t}\mathbf {A} (s)ds}\mathbf {c} ~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1883e6fbe661b2c18293f5081ab4f114225e13f1)
其中
是
常向量。
通過使用哈密爾頓–凱萊定理和類范德蒙矩陣,這種形式化的矩陣指數解可簡化為一種簡單的形式。[1]下面,我們將用普策算法(Putzer's algorithm)來展示這一方法。[2]
矩陣系統的穩定性與穩態[編輯]
矩陣方程
![{\displaystyle \mathbf {\dot {x}} (t)=\mathbf {Ax} (t)+\mathbf {b} }](https://wikimedia.org/api/rest_v1/media/math/render/svg/734b2624382dd9061a647cd3200071404a09ff62)
若且唯若常數矩陣A的所有特徵值的實部都為負,n×1參數常數向量b才穩定。
穩定時收斂到的穩態x*可置
![{\displaystyle \mathbf {\dot {x}} ^{*}(t)=\mathbf {0} ~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a87b6aaee9d4a3a10e5bb3c0abbac4ff534d9574)
找到,因此有
![{\displaystyle \mathbf {x} ^{*}=-\mathbf {A} ^{-1}\mathbf {b} ~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ebb2f9732c3954f2a87fae55128648f9920a6702)
假設A可逆。
因此,原方程可用偏離穩態的齊次形式來寫
![{\displaystyle \mathbf {\dot {x}} (t)=\mathbf {A} [\mathbf {x} (t)-\mathbf {x} ^{*}]~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/15618c92b5bbcaf82c8452dfbece582ee4351e53)
一種等效的表達是,x*是非齊次方程的一個特解,而所有解的形式都是
![{\displaystyle \mathbf {x} _{h}+\mathbf {x} ^{*}~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6ae27172a576935fc1a122f5082b9d2378956d56)
其中
是齊次方程(b=0)的解。
雙狀態變量情形的穩定性[編輯]
n = 2(2個狀態變量)時,穩定條件為:過渡矩陣A的兩個特徵值均有負實部,等價於A的跡為負、行列式為正。
矩陣形式的解[編輯]
的形式解為矩陣指數形式
![{\displaystyle \mathbf {x} (t)=\mathbf {x} ^{*}+e^{\mathbf {A} t}[\mathbf {x} (0)-\mathbf {x} ^{*}]~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f93c0eb9da46c79f6d69a716f41b19871ca9d465)
使用多種技術中的任何一種進行評估。
計算eAt的普策算法[編輯]
給定特徵值為
的矩陣A,
![{\displaystyle e^{\mathbf {A} t}=\sum _{j=0}^{n-1}r_{j+1}{\left(t\right)}\mathbf {P} _{j}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f34a39a9718e26b23f0851a4d5f4a1c5a11422b4)
其中
![{\displaystyle \mathbf {P} _{0}=\mathbf {I} }](https://wikimedia.org/api/rest_v1/media/math/render/svg/897f6dc6e4855e3c68a7dfd7810e23dc979d318e)
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ceabfc84e960b05ece9f9e162c01da8279ba54e6)
![{\displaystyle {\dot {r}}_{1}=\lambda _{1}r_{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9db93c31a9de3add834978c2a40b41ff1e171b2c)
![{\displaystyle r_{1}{\left(0\right)}=1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3592398092098e523ea1cb2d9a89d535681a79d8)
![{\displaystyle {\dot {r}}_{j}=\lambda _{j}r_{j}+r_{j-1},\qquad j=2,3,\dots ,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3228409747f2699c0bfc263bd0901e4dd93b29a9)
![{\displaystyle r_{j}{\left(0\right)}=0,\qquad j=2,3,\dots ,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/759fa131d8481c6a477eda94602627eab4406854)
的方程是簡單的一階非齊次常微分方程。
注意該算法並不要求矩陣A可對角化,並繞過了通常使用的若爾當標準形的計算。
矩陣常微分方程解構示例[編輯]
一階齊次矩陣常微分方程包含兩個函數x(t)、y(t),從矩陣形式解出後有如下形式:
![{\displaystyle {\frac {dx}{dt}}=a_{1}x+b_{1}y,\quad {\frac {dy}{dt}}=a_{2}x+b_{2}y}](https://wikimedia.org/api/rest_v1/media/math/render/svg/92b122a4bba5d83e89757bc7f57e1eb68d6131c6)
其中
、
、
、
可為任意純量。
高階矩陣ODE的形式可能複雜得多。
解分解後的矩陣常微分方程[編輯]
求解上述方程並找到這種特定階次和形式的所需函數的過程大概分3步。每個步驟的簡要說明如下:
第三部通常是把前兩步的結果代入專門形式的一般方程中,下詳。
矩陣ODE已解示例[編輯]
要按上述3步解矩陣ODE,並在過程中使用簡單矩陣,具體來說,現在下面的一階齊次線性ODE中找到函數x、函數y,都用單一自變量t表示:
![{\displaystyle {\frac {dx}{dt}}=3x-4y,\quad {\frac {dy}{dt}}=4x-7y~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/36da497d1f4644227a5499564afae68930798f02)
要解這個常微分方程系統,在過程中的某時刻需要一組兩個初始條件(對應起點的兩個狀態變量)。這時先取x(0) = y(0) = 1。
第一步[編輯]
第一步即找到A的特徵值
![{\displaystyle {\begin{bmatrix}x'\\y'\end{bmatrix}}={\begin{bmatrix}3&-4\\4&-7\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d6f877e95d04a130106d5cdd029681b2719485ea)
上面的導數記號x′等稱為拉格朗日記法(由約瑟夫·拉格朗日提出,等同於前面方程里的dx/dt,這是萊布尼茲記法,得名於戈特弗里德·萊布尼茨)。
一旦兩個變量的係數被寫為上述矩陣形式A,就可估計特徵值了。為此,可求矩陣行列式,即從上述係數矩陣中減去單位矩陣
乘常數λ,再得到特徵多項式
![{\displaystyle \det \left({\begin{bmatrix}3&-4\\4&-7\end{bmatrix}}-\lambda {\begin{bmatrix}1&0\\0&1\end{bmatrix}}\right)~,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/77abb733cdf0e34255e5031d1b4b48d112832a83)
再解得其零點。
進一步簡化、應用矩陣加法的基本規則,得出
![{\displaystyle \det {\begin{bmatrix}3-\lambda &-4\\4&-7-\lambda \end{bmatrix}}~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9f2c2c52fbcc32bb2c39e548f388fe713685ad62)
應用求單一2×2矩陣行列式的規則,可得下列一元二次方程
![{\displaystyle \det {\begin{bmatrix}3-\lambda &-4\\4&-7-\lambda \end{bmatrix}}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fa09e66c5a7e8773a8f76bd9715f746c9969c9c4)
![{\displaystyle -21-3\lambda +7\lambda +\lambda ^{2}+16=0\,\!}](https://wikimedia.org/api/rest_v1/media/math/render/svg/053670fc790cb69d930c985667a5d2844031fb3d)
可以進一步簡化
![{\displaystyle \lambda ^{2}+4\lambda -5=0~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a9b75ef3d260f95d6b00c6a0eb77a1ab7a6f2ece)
應用因式分解得到給定一元二次方程的兩個根
、
![{\displaystyle \lambda ^{2}+5\lambda -\lambda -5=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2bd93fc34999c785370bdc94243439f1e9d69d78)
![{\displaystyle \lambda (\lambda +5)-1(\lambda +5)=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/912c3c18ed6fddefe5d57b1011ad260496f4df70)
![{\displaystyle (\lambda -1)(\lambda +5)=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/689113a1f649e732bfe69cf5219cff64f048f47c)
![{\displaystyle \lambda =1,-5~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ce6849b096d837c60ac23b5c71daacbc693951de)
上面算出的
、
即所求A的特徵值。
矩陣ODE的特徵值可能是複數,求解過程的下一步及最終形式和解法可能會有巨大變化。
第二步[編輯]
第二步即找到A的特徵向量。
對算出的每個特徵值,都有單獨的特徵向量。例如對第一個特徵值即
,有
![{\displaystyle {\begin{bmatrix}3&-4\\4&-7\end{bmatrix}}{\begin{bmatrix}\alpha \\\beta \end{bmatrix}}=1{\begin{bmatrix}\alpha \\\beta \end{bmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/db3682d28082a46dc7f2921a6a407d2c3cad05f0)
應用矩陣乘法規則簡化上式,得到
![{\displaystyle 3\alpha -4\beta =\alpha }](https://wikimedia.org/api/rest_v1/media/math/render/svg/1c5e0b6427aef3466e7828672ebeba6321b4d97a)
![{\displaystyle \alpha =2\beta ~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a2d67a57c87b0fb4ed039bc45445929711bf4ddb)
所有計算都是為了得到最後一個式子,本例中就是α = 2β。現在任取一個無關緊要的小值(這樣更容易處理),代入α = 2β中的α或β(選哪個並不重要),這樣就得到了一個簡單的向量,就是這個特定特徵值所需的特徵向量。在本例中,我們取α = 2,得β = 1。用標準的向量符號來寫,向量是這樣的
![{\displaystyle \mathbf {\hat {v}} _{1}={\begin{bmatrix}2\\1\end{bmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bc2f8a81b3a7920b04af027a55fa2c06d99c2f81)
對第二個特徵值
進行相同的計算,得到第二個特徵向量,結果為
![{\displaystyle \mathbf {\hat {v}} _{2}={\begin{bmatrix}1\\2\end{bmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3d47ccf33a2c0f117a79e910aee462735a88f34c)
第三步[編輯]
最後一步是找到「隱藏」在導數背後的所求函數。有兩個函數,因為微分方程涉及兩個變量。
方程包含之前得到的所有信息,形式如下:
![{\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}=Ae^{\lambda _{1}t}\mathbf {\hat {v}} _{1}+Be^{\lambda _{2}t}\mathbf {\hat {v}} _{2}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f1353aa8d27522df41f0100c39853550fbde198c)
代入特徵值和特徵向量,得到
![{\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}=Ae^{t}{\begin{bmatrix}2\\1\end{bmatrix}}+Be^{-5t}{\begin{bmatrix}1\\2\end{bmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7024fb049c54dcb98157a079c62785753a3dfd3c)
簡化
![{\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}={\begin{bmatrix}2&1\\1&2\end{bmatrix}}{\begin{bmatrix}Ae^{t}\\Be^{-5t}\end{bmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0f7ca5afd738551e8f90e74a1c264035f4c730cc)
再簡化,分別寫出函數x、y的方程
![{\displaystyle x=2Ae^{t}+Be^{-5t}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dcb70f3622754a914ec385ab2b83f4d4c2542dc9)
![{\displaystyle y=Ae^{t}+2Be^{-5t}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c8edfb173a6867411a5a1916a28ae3e234e06a0e)
上述方程就是所求的一般函數,但只是一般形式(A、B的值未指定),但我們想找到它們的精確形式和解。因此現在,考慮問題的給定初始條件(即所謂初值問題)。假設給定了
,是ODE的起點;條件的應用指定了常數A、B。從條件
可以看出,t = 0時,上述方程的左式等於1,由此可構造下列線性方程組
![{\displaystyle 1=2A+B}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d4b5adb357455a69fbe47c24624866e9e0753c6c)
![{\displaystyle 1=A+2B~.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ecb0cdf3526ba33b9571dac489490cee8ac3e60f)
求解這些等式,發現常數A、B都等於1/3。因此將這些值代入這兩個函數的一般形式,就可以得到它們的精確形式
所求的兩個函數。
使用矩陣指數[編輯]
上述問題可以直接應用矩陣指數法解決。也就是說,可以說
給出了(可用MATLAB的expm
工具包之類,或通過對角化,並利用對角矩陣的矩陣指數與元素的指數化相等這一特性來計算)
得到最終解
這與之前展示的特徵向量方法相同。
參考文獻[編輯]