微分方程是变量的未知函数的数学方程,将函数值与不同阶导数联系起来。矩阵微分方程包含多个函数,以向量形式堆叠在一起,并由一个矩阵将它们与导数联系起来。
例如,一阶矩阵常微分方程为
![{\displaystyle \mathbf {\dot {x}} (t)=\mathbf {A} (t)\mathbf {x} (t)}](https://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/734b2624382dd9061a647cd3200071404a09ff62)
当且仅当常数矩阵A的所有特征值的实部都为负,n×1参数常数向量b才稳定。
稳定时收敛到的稳态x*可置
![{\displaystyle \mathbf {\dot {x}} ^{*}(t)=\mathbf {0} ~,}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/a87b6aaee9d4a3a10e5bb3c0abbac4ff534d9574)
找到,因此有
![{\displaystyle \mathbf {x} ^{*}=-\mathbf {A} ^{-1}\mathbf {b} ~,}](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/15618c92b5bbcaf82c8452dfbece582ee4351e53)
一种等效的表达是,x*是非齐次方程的一个特解,而所有解的形式都是
![{\displaystyle \mathbf {x} _{h}+\mathbf {x} ^{*}~,}](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/f93c0eb9da46c79f6d69a716f41b19871ca9d465)
使用多种技术中的任何一种进行评估。
给定特征值为
的矩阵A,
![{\displaystyle e^{\mathbf {A} t}=\sum _{j=0}^{n-1}r_{j+1}{\left(t\right)}\mathbf {P} _{j}}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/f34a39a9718e26b23f0851a4d5f4a1c5a11422b4)
其中
![{\displaystyle \mathbf {P} _{0}=\mathbf {I} }](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/ceabfc84e960b05ece9f9e162c01da8279ba54e6)
![{\displaystyle {\dot {r}}_{1}=\lambda _{1}r_{1}}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/9db93c31a9de3add834978c2a40b41ff1e171b2c)
![{\displaystyle r_{1}{\left(0\right)}=1}](https://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/92b122a4bba5d83e89757bc7f57e1eb68d6131c6)
其中
、
、
、
可为任意标量。
高阶矩阵ODE的形式可能复杂得多。
求解上述方程并找到这种特定阶次和形式的所需函数的过程大概分3步。每个步骤的简要说明如下:
第三部通常是把前两步的结果代入专门形式的一般方程中,下详。
要按上述3步解矩阵ODE,并在过程中使用简单矩阵,具体来说,现在下面的一阶齐次线性ODE中找到函数x、函数y,都用单一自变量t表示:
![{\displaystyle {\frac {dx}{dt}}=3x-4y,\quad {\frac {dy}{dt}}=4x-7y~.}](https://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/77abb733cdf0e34255e5031d1b4b48d112832a83)
再解得其零点。
进一步简化、应用矩阵加法的基本规则,得出
![{\displaystyle \det {\begin{bmatrix}3-\lambda &-4\\4&-7-\lambda \end{bmatrix}}~.}](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/fa09e66c5a7e8773a8f76bd9715f746c9969c9c4)
![{\displaystyle -21-3\lambda +7\lambda +\lambda ^{2}+16=0\,\!}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/053670fc790cb69d930c985667a5d2844031fb3d)
可以进一步简化
![{\displaystyle \lambda ^{2}+4\lambda -5=0~.}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/a9b75ef3d260f95d6b00c6a0eb77a1ab7a6f2ece)
应用因式分解得到给定一元二次方程的两个根
、
![{\displaystyle \lambda ^{2}+5\lambda -\lambda -5=0}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/2bd93fc34999c785370bdc94243439f1e9d69d78)
![{\displaystyle \lambda (\lambda +5)-1(\lambda +5)=0}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/912c3c18ed6fddefe5d57b1011ad260496f4df70)
![{\displaystyle (\lambda -1)(\lambda +5)=0}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/689113a1f649e732bfe69cf5219cff64f048f47c)
![{\displaystyle \lambda =1,-5~.}](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/db3682d28082a46dc7f2921a6a407d2c3cad05f0)
应用矩阵乘法规则简化上式,得到
![{\displaystyle 3\alpha -4\beta =\alpha }](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/1c5e0b6427aef3466e7828672ebeba6321b4d97a)
![{\displaystyle \alpha =2\beta ~.}](https://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/bc2f8a81b3a7920b04af027a55fa2c06d99c2f81)
对第二个特征值
进行相同的计算,得到第二个特征向量,结果为
![{\displaystyle \mathbf {\hat {v}} _{2}={\begin{bmatrix}1\\2\end{bmatrix}}.}](https://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/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://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/0f7ca5afd738551e8f90e74a1c264035f4c730cc)
再简化,分别写出函数x、y的方程
![{\displaystyle x=2Ae^{t}+Be^{-5t}}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/dcb70f3622754a914ec385ab2b83f4d4c2542dc9)
![{\displaystyle y=Ae^{t}+2Be^{-5t}.}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/c8edfb173a6867411a5a1916a28ae3e234e06a0e)
上述方程就是所求的一般函数,但只是一般形式(A、B的值未指定),但我们想找到它们的精确形式和解。因此现在,考虑问题的给定初始条件(即所谓初值问题)。假设给定了
,是ODE的起点;条件的应用指定了常数A、B。从条件
可以看出,t = 0时,上述方程的左式等于1,由此可构造下列线性方程组
![{\displaystyle 1=2A+B}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/d4b5adb357455a69fbe47c24624866e9e0753c6c)
![{\displaystyle 1=A+2B~.}](https://backend.710302.xyz/https/wikimedia.org/api/rest_v1/media/math/render/svg/ecb0cdf3526ba33b9571dac489490cee8ac3e60f)
求解这些等式,发现常数A、B都等于1/3。因此将这些值代入这两个函数的一般形式,就可以得到它们的精确形式
所求的两个函数。
上述问题可以直接应用矩阵指数法解决。也就是说,可以说
给出了(可用MATLAB的expm
工具包之类,或通过对角化,并利用对角矩阵的矩阵指数与元素的指数化相等这一特性来计算)
得到最终解
这与之前展示的特征向量方法相同。