矩阵指数解线性微分方程

最新推荐文章于 2024-05-04 00:19:27 发布

BPSSY

最新推荐文章于 2024-05-04 00:19:27 发布

阅读量6.3k

点赞数 1

分类专栏：工控

工控专栏收录该内容

3 篇文章

订阅专栏

矩阵指数

维基百科，自由的百科全书

矩阵指数是方块矩阵的一种矩阵函数，与指数函数类似。矩阵指数给出了矩阵李代数与对应的李群之间的关系。

设X为n×n的实数或复数矩阵。X的指数，用e^X或exp(X)来表示，是由以下幂级数所给出的n×n矩阵：

$e^X = \sum_{k=0}^\infty{1 \over k!}X^k.$

以上的级数总是收敛的，因此X的指数是定义良好的。注意，如果X是1×1的矩阵，则X的矩阵指数就是由X的元素的指数所组成的1×1矩阵。

性质

设X和Y为n×n的复数矩阵，并设a和b为任意的复数。我们把n×n的单位矩阵记为I，把零矩阵记为0。矩阵指数满足以下性质：

e⁰ = I。
e^aXe^bX = e⁽^{a + b)X}。
e^Xe^−X = I。
如果XY = YX，则e^Xe^Y = e^Ye^X = e^(X + Y)。
如果Y是可逆的，则e^YXY⁻¹ = Ye^XY⁻¹。
det(e^X) = e^tr(X)，其中tr(X)是X的迹数。
exp(X^T) = (exp X)^T，其中X^T表示X的转置。从中可以推出，如果X是对称矩阵，则e^X也是对称矩阵；如果X是斜对称矩阵，则e^X是正交矩阵。
exp(X*) = (exp X)*，其中X*表示X的共轭转置。可以推出，如果X是埃尔米特矩阵，则e^X也是埃尔米特矩阵；如果X是斜埃尔米特矩阵，则e^X是酉矩阵。

线性微分方程

矩阵指数的一个重要性，是它可以用来解微分方程。从(1)可知，以下微分方程

$\frac{d}{dt} y(t) = Ay(t), \quad y(0) = y_0,$

其中A是矩阵，具有解

$y(t) = e^{At} y_0. \,$

矩阵指数也可以用来解非齐次方程：

$\frac{d}{dt} y(t) = Ay(t) + z(t), \quad y(0) = y_0.$

参见以下的例子。

当A不是常数时，以下形式的微分方程没有闭式解：

$\frac{d}{dt} y(t) = A(t) \, y(t), \quad y(0) = y_0,$

但马格努斯级数可以给出无穷级数形式的解。

指数相加

我们知道，对于任何实数（标量）x和y，指数函数都满足公式e^{x + y} = e^xe^y。类似的等式对于可交换矩阵也成立：如果矩阵X和Y是可交换的（意思是说XY = YX），则：

$e^{X+Y} = e^Xe^Y. \,$

但是，如果它们不是可交换的，则以上的等式不一定成立。在这种情况下，我们可以用Baker-Campbell-Hausdorff公式来计算e^X + Y。

这个命题反过来不成立：e^X + Y = e^Xe^Y并不一定就意味着X和Y是可交换的。但是，如果X和Y只含有代数数，而且它们的大小至少为2×2，则反过来也成立。（Horn & Johnson 1991，pp.435–437）.

指数映射

注意矩阵的指数总是非奇异方阵。e^X的逆矩阵由e^−X给出。这与复数的指数总是非零的事实类似。这样，矩阵指数就给出了一个映射：

$\exp \colon M_n(\mathbb C) \to \mbox{GL}(n,\mathbb C)$

这是从所有n×n矩阵的空间到一般线性群（所有非奇异方阵所组成的群）的映射。实际上，这个映射是满射，就是说每一个非奇异方阵都可以写成某个矩阵的指数。矩阵对数就是这个映射的逆映射。

对于任何两个矩阵X和Y，我们有：

$\| e^{X+Y} - e^X \| \le \|Y\| e^{\|X\|} e^{\|Y\|},$

其中|| · ||表示任意的矩阵范数。从中可以推出，指数映射在M_n(C)的紧子集内是连续和利普希茨连续的。

以下的映射

$t \mapsto e^{tX}, \qquad t \in \mathbb R$

定义了一般线性群中的一条光滑曲线，当t = 0时穿过单位元。实际上，这给出了一般线性群的一个单参数子群，这是由于：

$e^{tX}e^{sX} = e^{(t+s)X}.\,$

这条曲线在点t的导数（或切向量）由以下等式给出：

$\frac{d}{dt}e^{tX} = Xe^{tX}. \qquad (1)$

t = 0时的导数就是矩阵X，所以我们可以说，X是这个单参数子群的推广。

更加一般地：

$\frac{d}{dt}e^{X(t)} = \int_0^1 e^{(1-\alpha) X(t)} \frac{dX(t)}{dt} e^{\alpha X(t)}\,d\alpha.$

矩阵指数的计算

寻找计算矩阵指数的可靠和准确的方法是困难的，目前在数学和数值分析领域中仍然是一个正在研究的话题。有些方法列举如下。

可对角化矩阵

如果矩阵是对角的：

$A=\begin{bmatrix} a_1 & 0 & \ldots & 0 \\0 & a_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & a_n \end{bmatrix},$

则把主对角线上的所有元素取指数，就是原矩阵的指数：

$e^A=\begin{bmatrix} e^{a_1} & 0 & \ldots & 0 \\0 & e^{a_2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & e^{a_n} \end{bmatrix}.$

这也允许了我们计算可对角化矩阵的指数。如果 $A = UDU^{-1}$ ，且D是对角矩阵，则 $e^A = U e^D U^{-1}$ 。用西尔维斯特公式，也可以得到相同的结果。

幂零矩阵

如果对于某个整数q，有N^q = 0，则矩阵N称为幂零矩阵。在这种情况下，矩阵指数e^N可以直接从级数展开式来计算，这是因为级数在有限个项后就终止了：

$e^N = I + N + \frac{1}{2}N^2 + \frac{1}{6}N^3 + \cdots + \frac{1}{(q-1)!}N^{q-1}.$

推广

当矩阵X的最小多项式可以分解为一次多项式的积时，它就可以表示为以下的和：

$X = A + N \,$

其中：

A是可对角化矩阵；
N是幂零矩阵；
A与N是可交换的（也就是说， AN = NA）。

这称为Dunford分解。

这就是说，我们可以通过化为前两种情况，来计算X的指数：

$e^X = e^{A+N} = e^A e^N. \,$

注意为了让最后一步成立， A和N必须是可交换的。

另外一个密切相关的方法，是利用X的若尔当标准型。假设X = PJP⁻¹，其中J是X的若尔当标准型。那么：

$e^{X}=Pe^{J}P^{-1}.\,$

另外，由于

$J=J_{a_1}(\lambda_1)\oplus J_{a_2}(\lambda_2)\oplus\cdots\oplus J_{a_n}(\lambda_n),$

$\begin{align}e^{J} & {} = \exp \big( J_{a_1}(\lambda_1)\oplus J_{a_2}(\lambda_2)\oplus\cdots\oplus J_{a_n}(\lambda_n) \big) \\& {} = \exp \big( J_{a_1}(\lambda_1) \big) \oplus \exp \big( J_{a_2}(\lambda_2) \big) \oplus\cdots\oplus \exp \big( J_{a_k}(\lambda_k) \big).\end{align}$

因此，我们只需要知道怎样计算若尔当块的矩阵指数。但是，每一个若尔当块都具有形式

$J_{a}(\lambda) = \lambda I + N \,$

其中N是幂零矩阵。则这个区块的矩阵指数由下式给出：

$e^{\lambda I + N} = e^{\lambda}e^N. \,$

计算

假设我们想要计算以下矩阵的指数。

$B=\begin{bmatrix}21 & 17 & 6 \\-5 & -1 & -6 \\4 & 4 & 16 \end{bmatrix}.$

它的若尔当型为：

$J = P^{-1}BP = \begin{bmatrix}4 & 0 & 0 \\0 & 16 & 1 \\0 & 0 & 16 \end{bmatrix},$

其中矩阵P由下式给出：

$P=\begin{bmatrix}-\frac14 & 2 & \frac54 \\\frac14 & -2 & -\frac14 \\0 & 4 & 0 \end{bmatrix}.$

我们首先来计算exp(J)。我们有：

$J=J_1(4)\oplus J_2(16)$

1×1矩阵的指数仅仅是该矩阵的元素的指数，因此exp(J₁(4)) = [e⁴]。 $J_2(16)$ 的指数可以用以上提到的公式exp(λ $I$ +N) = e^λ exp(N)来算出：

$\exp \left( \begin{bmatrix} 16 & 1 \\ 0 & 16 \end{bmatrix} \right) = e^{16} \exp \left( \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \right) = e^{16} \left(\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} + \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} + {1 \over 2!}\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix} + \cdots \right)= \begin{bmatrix} e^{16} & e^{16} \\ 0 & e^{16} \end{bmatrix}.$

因此，原矩阵B的指数为：

$\exp(B) = P \exp(J) P^{-1} = P \begin{bmatrix} e^4 & 0 & 0 \\ 0 & e^{16} & e^{16} \\ 0 & 0 & e^{16} \end{bmatrix} P^{-1} = {1\over 4} \begin{bmatrix} 13e^{16} - e^4 & 13e^{16} - 5e^4 & 2e^{16} - 2e^4 \\ -9e^{16} + e^4 & -9e^{16} + 5e^4 & -2e^{16} + 2e^4 \\ 16e^{16} & 16e^{16} & 4e^{16} \end{bmatrix}.$

应用

线性微分方程

矩阵指数在解线性微分方程时十分有用。前面曾提到，以下形式的微分方程

$\mathbf{y}' = C\mathbf{y}$

具有解e^Cty(0)。如果我们考虑以下向量

$\mathbf{y}(t) = \begin{pmatrix} y_1(t) \\ \vdots \\y_n(t) \end{pmatrix}$

我们就可以把线性微分方程表示为：

$\mathbf{y}'(t) = A\mathbf{y}(t)+\mathbf{b}(t).$

如果我们作一个猜想，把两边乘以一个积分因子 e^−At，便得到：

$e^{-At}\mathbf{y}'-e^{-At}A\mathbf{y} = e^{-At}\mathbf{b}$

$\frac{d}{dt} (e^{-At}\mathbf{y}) = e^{-At}\mathbf{b}.$

如果我们可以计算e^At，那么就得到了微分方程的解。

例子（齐次）

假设我们有以下的微分方程组：

$\begin{matrix}x' &=& 2x&-y&+z \\y' &=& &3y&-1z \\z' &=& 2x&+y&+3z \end{matrix}$

非齐次的情况──参数变换

对于非齐次的情况，我们可以用积分因子的方法（类似于参数变换的方法）。我们找到形为y_p(t) = exp(tA)z(t)一个特解：

$\mathbf{y}_p' = (e^{tA})'\mathbf{z}(t)+e^{tA}\mathbf{z}'(t)$

$= Ae^{tA}\mathbf{z}(t)+e^{tA}\mathbf{z}'(t)$

$= A\mathbf{y}_p(t)+e^{tA}\mathbf{z}'(t)$

为了让y_p为方程的解，必须有：

$e^{tA}\mathbf{z}'(t) = \mathbf{b}(t)$

$\mathbf{z}'(t) = (e^{tA})^{-1}\mathbf{b}(t)$

$\mathbf{z}(t) = \int_0^t e^{-uA}\mathbf{b}(u)\,du+\mathbf{c}$

因此，

$\begin{align}\mathbf{y}_p & {} = e^{tA}\int_0^t e^{-uA}\mathbf{b}(u)\,du+e^{tA}\mathbf{c} \\& {} = \int_0^t e^{(t-u)A}\mathbf{b}(u)\,du+e^{tA}\mathbf{c}\end{align}$

其中c由问题的初始条件决定。

例子（非齐次）

假设我们有以下的微分方程组：

$\begin{matrix}x' &=& 2x&-y&+z&+e^{2t} \\y' &=& &3y&-1z& \\z' &=& 2x&+y&+3z&+e^{2t} \end{matrix}$

那么我们有

$M=\begin{bmatrix}2 & -1 & 1 \\0 & 3 & -1 \\2 & 1 & 3 \end{bmatrix}$

以及

$\mathbf{b}=e^{2t}\begin{bmatrix}1 \\0\\1\end{bmatrix}.$

用前面的方法，我们可以得出齐次微分方程的解。由于齐次方程的通解与非齐次方程的特解的和就是非齐次方程的通解，因此我们只需要找到一个特解（用参数变换法）。

我们有：

$\mathbf{y}_p = e^{t}\int_0^t e^{(-u)A}\begin{bmatrix}e^{2u} \\0\\e^{2u}\end{bmatrix}\,du+e^{tA}\mathbf{c}$

$\mathbf{y}_p = e^{t}\int_0^t\begin{bmatrix} 2e^u - 2ue^{2u} & -2ue^{2u} & 0 \\ \\-2e^u + 2(u+1)e^{2u} & 2(u+1)e^{2u} & 0 \\ \\ 2ue^{2u} & 2ue^{2u} & 2e^u\end{bmatrix}\begin{bmatrix}e^{2u} \\0\\e^{2u}\end{bmatrix}\,du+e^{tA}\mathbf{c}$

$\mathbf{y}_p = e^{t}\int_0^t\begin{bmatrix}e^{2u}( 2e^u - 2ue^{2u}) \\ \\ e^{2u}(-2e^u + 2(1 + u)e^{2u}) \\ \\ 2e^{3u} + 2ue^{4u}\end{bmatrix}+e^{tA}\mathbf{c}$

$\mathbf{y}_p = e^{t}\begin{bmatrix}-{1 \over 24}e^{3t}(3e^t(4t-1)-16) \\ \\{1 \over 24}e^{3t}(3e^t(4t+4)-16) \\ \\{1 \over 24}e^{3t}(3e^t(4t-1)-16)\end{bmatrix}+\begin{bmatrix} 2e^t - 2te^{2t} & -2te^{2t} & 0 \\ \\-2e^t + 2(t+1)e^{2t} & 2(t+1)e^{2t} & 0 \\ \\ 2te^{2t} & 2te^{2t} & 2e^t\end{bmatrix}\begin{bmatrix}c_1 \\c_2 \\c_3\end{bmatrix}$