查看“︁矩阵微分方程”︁的源代码

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

例如，一阶矩阵[[常微分方程]]为
: <math>\mathbf{\dot{x}}(t) = \mathbf{A}(t)\mathbf{x}(t)</math>
其中<math>\mathbf{x}(t)</math>是基变量<math>t</math>的函数的<math>n \times 1</math>向量，<math>\mathbf{\dot{x}}(t)</math>是函数的一阶导，<math>\mathbf{A}(t)</math>是<math>n \times n</math>系数矩阵。

若<math>\mathbf{A}</math>为常数且有''n''个[[线性无关]]的特征向量，微分方程有如下一般解：
: <math>\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  ~,</math>
其中{{nowrap|''λ''<sub>1</sub>, ''λ''<sub>2</sub>, …, ''λ''<sub>''n''</sub>}}是'''A'''的特征值，{{nowrap|'''u'''<sub>1</sub>, '''u'''<sub>2</sub>, …, '''u'''<sub>''n''</sub>}}是'''A'''相应的特征向量；{{nowrap|''c''<sub>1</sub>, ''c''<sub>2</sub>, …, ''c''<sub>''n''</sub>}}为常数。

更一般地说，若<math>\mathbf{A}(t)</math>等于其积分<math>\int_a^t \mathbf{A}(s)ds</math>，则[[马格努斯展开]]降为前导阶，微分方程的一般解是
: <math>\mathbf{x}(t)=e^{\int_a^t \mathbf{A}(s) ds} \mathbf{c} ~,</math>
其中<math> \mathbf{c} </math>是<math> n \times 1 </math>常向量。

通过使用[[哈密尔顿–凯莱定理]]和类[[范德蒙矩阵]]，这种形式化的[[矩阵指数]]解可简化为一种简单的形式。<ref>{{cite book |first=H. |last=Moya-Cessa |first2=F. |last2=Soto-Eguibar |title=Differential Equations: An Operational Approach |publisher=Rinton Press |location=New Jersey |year=2011 |isbn=978-1-58949-060-4 }}</ref>下面，我们将用普策算法（Putzer's algorithm）来展示这一方法。<ref>{{cite journal |first=E. J. |last=Putzer |title=Avoiding the Jordan Canonical Form in the Discussion of Linear Systems with Constant Coefficients |url=https://archive.org/details/sim_american-mathematical-monthly_1966-01_73_1/page/2 |journal=The American Mathematical Monthly |volume=73 |issue=1 |year=1966 |pages=2–7 |jstor=2313914 |doi=10.1080/00029890.1966.11970714 }}</ref>

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

矩阵方程
:<math>\mathbf{\dot{x}}(t) = \mathbf{Ax}(t) + \mathbf{b}</math> 
当且仅当常数矩阵'''A'''的所有特征值的实部都为负，''n''×1参数常数向量'''b'''才[[稳定性理论|稳定]]。

稳定时收敛到的稳态'''x*'''可置
:<math>\mathbf{\dot{x}}^* (t)=\mathbf{0}~,</math>
找到，因此有
:<math>\mathbf{x}^* = -\mathbf{A}^{-1}\mathbf{b}~,</math>
假设'''A'''可逆。

因此，原方程可用偏离稳态的齐次形式来写
:<math> \mathbf{\dot{x}}(t)=\mathbf{A}[\mathbf{x}(t)-\mathbf{x}^*]~.</math>

一种等效的表达是，'''x*'''是非齐次方程的一个特解，而所有解的形式都是
:<math>\mathbf{x}_h+\mathbf{x}^*   ~,</math> 
其中<math>\mathbf{x}_h</math>是齐次方程（'''b'''='''0'''）的解。

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

==矩阵形式的解==
<math>\mathbf{\dot{x}}(t)=\mathbf{A}[\mathbf{x}(t)-\mathbf{x}^*]</math>的形式解为[[矩阵指数]]形式
:<math>\mathbf{x}(t)=\mathbf{x}^*+e^{\mathbf{A}t}[\mathbf{x}(0)-\mathbf{x}^*] ~,</math> 
使用多种技术中的任何一种进行评估。

===计算{{math|e<sup>'''A'''''t''</sup>}}的普策算法===
给定特征值为<math>\lambda_1,\lambda_2,\dots,\lambda_n</math>的矩阵'''A'''，
:<math>e^{\mathbf{A}t} = \sum_{j=0}^{n-1}r_{j+1}{\left(t\right)}\mathbf{P}_{j}</math>
其中
:<math>\mathbf{P}_0 = \mathbf{I}</math>
:<math>\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</math>
:<math>\dot{r}_1 = \lambda_1 r_1</math>
:<math>r_1{\left(0\right)}=1</math>
:<math>\dot{r}_{j} = \lambda_j r_j + r_{j-1}, \qquad  j=2,3,\dots,n</math>
:<math>r_j{\left(0\right)}=0,  \qquad  j=2,3,\dots,n</math>

<math>r_i (t)</math>的方程是简单的一阶非齐次常微分方程。

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

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

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

: <math>\frac{dx}{dt}=a_1x+b_1y,\quad\frac{dy}{dt}=a_2x+b_2y</math>

其中<math>a_1</math>、<math>a_2</math>、<math>b_1</math>、<math>b_2</math>可为任意标量。
高阶矩阵ODE的形式可能复杂得多。
== 解分解后的矩阵常微分方程 ==

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

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

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

==矩阵ODE已解示例==
{{see also|矩阵指数#线性微分方程_2}}
要按上述3步解矩阵ODE，并在过程中使用简单矩阵，具体来说，现在下面的一阶齐次线性ODE中找到函数{{mvar|x}}、函数{{mvar|y}}，都用单一自变量{{mvar|t}}表示：
: <math>\frac{dx}{dt}=3x-4y,\quad\frac{dy}{dt}=4x-7y~.</math>

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

===第一步===
第一步即找到'''A'''的特征值
: <math>\begin{bmatrix} x'\\y' \end{bmatrix} = \begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix}\begin{bmatrix} x\\y \end{bmatrix}~.</math>
上面的[[导数]]记号''x''′等称为拉格朗日记法（由[[约瑟夫·拉格朗日]]提出，等同于前面方程里的''dx/dt''，这是[[莱布尼兹记法]]，得名于[[戈特弗里德·莱布尼茨]]）。

一旦两个变量的[[系数]]被写为上述[[矩阵]]形式'''A'''，就可估计特征值了。为此，可求矩阵[[行列式]]，即从上述系数矩阵中减去[[单位矩阵]]<math>I_n</math>乘常数{{mvar|λ}}，再得到[[特征多项式]]
: <math>\det\left(\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix} - \lambda\begin{bmatrix} 1 & 0\\0 & 1 \end{bmatrix}\right)~,</math>
再解得其零点。

进一步简化、应用[[矩阵加法]]的基本规则，得出
: <math>\det\begin{bmatrix} 3-\lambda & -4\\4 & -7-\lambda \end{bmatrix}~.</math>

应用求单一2×2矩阵行列式的规则，可得下列一元二次方程
: <math>\det\begin{bmatrix} 3-\lambda & -4\\4 & -7-\lambda \end{bmatrix} = 0</math>
: <math>-21 - 3\lambda + 7\lambda + \lambda^2 + 16 = 0  \,\!</math>
可以进一步简化
: <math>\lambda^2 + 4\lambda - 5 = 0  ~.</math>

应用[[因式分解]]得到给定一元二次方程的两个根<math>\lambda_1</math>、<math>\lambda_2</math>
: <math>\lambda^2 + 5\lambda - \lambda - 5 = 0</math>
: <math>\lambda (\lambda + 5) - 1 (\lambda + 5) = 0</math>
: <math>(\lambda - 1)(\lambda + 5) = 0</math>

: <math>\lambda = 1, -5 ~.</math>
上面算出的<math>\lambda_1 = 1</math>、<math>\lambda_2 = -5</math>即所求'''A'''的特征值。
矩阵ODE的特征值可能是[[复数 (数学)|复数]]，求解过程的下一步及最终形式和解法可能会有巨大变化。

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

对算出的每个特征值，都有单独的特征向量。例如对第一个特征值即<math>\lambda_1 = 1</math>，有
: <math>\begin{bmatrix} 3 & -4\\4 & -7 \end{bmatrix}\begin{bmatrix} \alpha\\\beta \end{bmatrix} = 1\begin{bmatrix} \alpha\\\beta \end{bmatrix}. </math>
应用[[矩阵乘法]]规则简化上式，得到
: <math>3\alpha - 4\beta = \alpha</math>
: <math>\alpha = 2\beta~.</math>

所有计算都是为了得到最后一个式子，本例中就是{{math|1=''α'' = 2''β''}}。现在任取一个无关紧要的小值（这样更容易处理），代入{{math|1=''α'' = 2''β''}}中的{{mvar|α}}或{{mvar|β}}（选哪个并不重要），这样就得到了一个简单的向量，就是这个特定特征值所需的特征向量。在本例中，我们取{{math|1=''α'' = 2}}，得{{math|1=''β'' = 1}}。用标准的向量符号来写，向量是这样的
: <math>\mathbf{\hat{v}}_1 = \begin{bmatrix} 2\\1 \end{bmatrix}. </math>

对第二个特征值<math>\lambda = -5</math>进行相同的计算，得到第二个特征向量，结果为
: <math>\mathbf{\hat{v}}_2 = \begin{bmatrix} 1\\2 \end{bmatrix}. </math>

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

方程包含之前得到的所有信息，形式如下：
: <math>\begin{bmatrix} x\\y \end{bmatrix} = Ae^{\lambda_1t}\mathbf{\hat{v}}_1 + Be^{\lambda_2t}\mathbf{\hat{v}}_2. </math>

代入特征值和特征向量，得到
: <math>\begin{bmatrix} x\\y \end{bmatrix} = Ae^{t}\begin{bmatrix} 2\\1 \end{bmatrix} + Be^{-5t}\begin{bmatrix} 1\\2 \end{bmatrix}. </math>
简化
: <math>\begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 2 & 1\\1 & 2 \end{bmatrix}\begin{bmatrix} Ae^{t}\\Be^{-5t} \end{bmatrix}. </math>
再简化，分别写出函数{{mvar|x}}、{{mvar|y}}的方程
: <math>x = 2Ae^{t} + Be^{-5t} </math>
: <math>y = Ae^{t} + 2Be^{-5t}.</math>

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

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

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

<math>\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} </math>

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

<math>\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} </math>

得到最终解

<math>\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} </math>

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

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

== 另见 ==
* [[齐次微分方程]]
* [[矩阵差分方程]]
* [[冷却定律]]
* [[斐波那契数列]]
* [[差分方程]]
* [[波动方程]]
* [[自治系统 (数学)]]

==参考文献==
{{Reflist}}

[[Category:常微分方程]]