贡献者： addis
- 和一阶常系数线性微分方程组（常微分方程）部分内容重复
1. 常系数情况
一阶常微分方程组
\begin{equation}
\frac{\mathrm{d}{ \boldsymbol{\mathbf{v}} }}{\mathrm{d}{t}} = \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{v}} ~
\end{equation}
的解析解（先假设 $ \boldsymbol{\mathbf{A}} $ 为常矩阵）为
\begin{equation}
\boldsymbol{\mathbf{v}} (t) = \exp\left( \boldsymbol{\mathbf{A}} t\right) \boldsymbol{\mathbf{v}} (0)~.
\end{equation}
其中矩阵（这里是方阵）的指数函数也是一个大小相同的矩阵，有类似于泰勒级数
定义
\begin{equation}
\exp\left( \boldsymbol{\mathbf{A}} t\right) = 1 + \boldsymbol{\mathbf{A}} t + \frac1{2!} ( \boldsymbol{\mathbf{A}} t)^2 + \dots~
\end{equation}
代入即可验证
式 1 。这类似于一阶常系数常微分方程的解（未完成），即 $ \boldsymbol{\mathbf{A}} $ 是一个常数而不是矩阵的情况。
但式 3 不方便直接计算。此时如果 $ \boldsymbol{\mathbf{A}} $ 可以对角化为（链接未完成；任何矩阵都可以对角化吗？）
\begin{equation}
\boldsymbol{\mathbf{A}} = \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{U}} ^{-1}~,
\end{equation}
其中 $ \boldsymbol{\mathbf{U}} $ 是本征列向量排成的矩阵，$ \boldsymbol{\mathbf{\Lambda}} $ 是对应本征值 $\lambda_1, \lambda_2, \dots$ 排成的对角矩阵，
代入
式 3 就有
\begin{equation} \begin{aligned}
\exp\left( \boldsymbol{\mathbf{A}} t\right) &= \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{U}} ^{-1} + \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{U}} ^{-1} t + \frac1{2!} ( \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{U}} ^{-1} t)^2 + \dots\\
&= \boldsymbol{\mathbf{U}} \left(1 + \boldsymbol{\mathbf{\Lambda}} t + \frac1{2!} ( \boldsymbol{\mathbf{\Lambda}} t)^2 + \dots \right) \boldsymbol{\mathbf{U}} ^{-1}\\
&= \boldsymbol{\mathbf{U}} \exp\left( \boldsymbol{\mathbf{\Lambda}} t\right) \boldsymbol{\mathbf{U}} ^{-1}\\
&= \boldsymbol{\mathbf{U}} \begin{pmatrix} \mathrm{e} ^{\lambda_1 t} & 0 & \dots\\ 0 & \mathrm{e} ^{\lambda_2 t} &\dots\\ \vdots & \vdots & \ddots\end{pmatrix} \boldsymbol{\mathbf{U}} ^{-1}~.
\end{aligned} \end{equation}
代入
式 2 就是方程组的解。
2. 含时系数情况的形式解
当式 1 中的系数矩阵 $ \boldsymbol{\mathbf{A}} $ 是 $t$ 的函数 $ \boldsymbol{\mathbf{A}} (t)$ 时，一般没有解析解。我们可以取微小时间步长 $\Delta t$，在每个 $\Delta t$ 内近似认为 $ \boldsymbol{\mathbf{A}} (t_i)$ 为常数，再取极限
\begin{equation}
\boldsymbol{\mathbf{v}} (t) = \lim_{\Delta t \to 0} \prod \exp \left[ \boldsymbol{\mathbf{A}} (t_i)\Delta t \right] \boldsymbol{\mathbf{v}} (0)~,
\end{equation}
如果两个矩阵 $ \boldsymbol{\mathbf{P}} , \boldsymbol{\mathbf{Q}} $ 对易
，就有
\begin{equation}
\exp\left( \boldsymbol{\mathbf{P}} + \boldsymbol{\mathbf{Q}} \right) = \exp\left( \boldsymbol{\mathbf{P}} \right) \exp\left( \boldsymbol{\mathbf{Q}} \right) ~,
\end{equation}
但一般来说 $ \boldsymbol{\mathbf{A}} (t_i)$ 之间不对易，所以我们形式上定义一个时间排序算符 $ \hat{\mathcal T} $ 使例如
\begin{equation}
\hat{\mathcal T} [ \boldsymbol{\mathbf{A}} (t_1) \boldsymbol{\mathbf{A}} (t_3) \boldsymbol{\mathbf{A}} (k_2)] = \boldsymbol{\mathbf{A}} (t_3) \boldsymbol{\mathbf{A}} (t_2) \boldsymbol{\mathbf{A}} (t_1) \qquad ( t_1 < t_2 < t_3 )~,
\end{equation}
这样通解在形式上就可以记为
\begin{equation}
\boldsymbol{\mathbf{v}} (t) = \hat{\mathcal T} \exp\left(\int_0^ t \boldsymbol{\mathbf{A}} (t') \,\mathrm{d}{t} '\right) \boldsymbol{\mathbf{v}} (0)~,
\end{equation}
然而这么做对于数值计算并没有太大意义。
3. 厄米矩阵、反厄米矩阵
量子力学中经常需要求解矩阵形式的薛定谔方程（链接未完成）
\begin{equation}
\mathrm{i} \frac{\mathrm{d}}{\mathrm{d}{t}} \boldsymbol{\mathbf{v}} = \boldsymbol{\mathbf{H}} \boldsymbol{\mathbf{v}} ~.
\end{equation}
其中 $ \boldsymbol{\mathbf{H}} $ 是厄米矩阵（哈密顿矩阵）。若把两边同除以 $ \mathrm{i} $，则系数矩阵 $- \mathrm{i} \boldsymbol{\mathbf{H}} $ 就是一个反厄米矩阵。如果 $ \boldsymbol{\mathbf{H}} $ 不随时间变化，那么 $ \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t}$（$t\in \mathbb R$）是一个酉矩阵（幺正矩阵）
（链接未完成）。所以 $ \boldsymbol{\mathbf{v}} (t)$ 的模 $ \left\lVert \boldsymbol{\mathbf{v}} (t) \right\rVert = \sum \left\lvert v_i \right\rvert ^2$ 将不随时间变化。
证明：
\begin{equation}
\left\lVert \boldsymbol{\mathbf{v}} (t) \right\rVert = \boldsymbol{\mathbf{v}} (t) ^\dagger \boldsymbol{\mathbf{v}} (t) = \left( \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} t} \right) ^\dagger \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} t} \boldsymbol{\mathbf{v}} (0) \left\lVert \boldsymbol{\mathbf{v}} (0) \right\rVert ~.
\end{equation}
两边取厄米共轭得
\begin{equation}
\left( \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} t} \right) ^\dagger = \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} ^\dagger t} = \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t}~,
\end{equation}
由于 $[H,H] = 0$，有 $( \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t}) ^\dagger \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t} = \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} t} \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t} = \boldsymbol{\mathbf{I}} $。证毕。
4. 数值计算
事实上，以上做法相当于分离变量，当 $ \boldsymbol{\mathbf{A}} $ 是厄米矩阵时，令 $ \boldsymbol{\mathbf{v}} (t) = f(t) \boldsymbol{\mathbf{u}} $，代入方程得
\begin{equation}
\frac{f'(t)}{f(t)} \boldsymbol{\mathbf{u}} = \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{u}} ~,
\end{equation}
由于 $ \boldsymbol{\mathbf{A}} $ 和 $ \boldsymbol{\mathbf{u}} $ 都不含时，所以可以令
\begin{equation} \begin{aligned}
& \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{u}} = \lambda \boldsymbol{\mathbf{u}} ~,\\
& f(t)' = \lambda f(t)~.
\end{aligned} \end{equation}
其中第一个方程是 $ \boldsymbol{\mathbf{A}} $ 的本征方程，解为 $N$ 个本征矢 $ \boldsymbol{\mathbf{u}} _i$（即 $ \boldsymbol{\mathbf{U}} $ 的第 $i$ 列）和 $N$ 个本征值 $\lambda_i$（即 $ \boldsymbol{\mathbf{\Lambda}} $ 的第 $i$ 个对角元）。第二条方程的解为 $f(t) = \exp\left(\lambda t\right) $。
在此基础上使用 Lanczos 算法可以进一步提高效率。
5. 算符拆分
有时候我们希望可以在上述计算中把 $ \boldsymbol{\mathbf{A}} $ 写成几个矩阵的和的形式（以两个为例）$ \boldsymbol{\mathbf{A}} = \boldsymbol{\mathbf{B}} + \boldsymbol{\mathbf{C}} $。当 $ \boldsymbol{\mathbf{B}} $ 和 $ \boldsymbol{\mathbf{C}} $ 对易时显然有
\begin{equation}
\exp\left( \boldsymbol{\mathbf{A}} t\right) = \exp\left( \boldsymbol{\mathbf{B}} t\right) \exp\left( \boldsymbol{\mathbf{C}} t\right) ~.
\end{equation}
在程序中这么做可能可以进一步提高速度
1。如果 $ \boldsymbol{\mathbf{B}} $ 和 $ \boldsymbol{\mathbf{C}} $ 不对易，严格来说上式不成立，但可以证明 $t \to 0$ 时近似成立
\begin{equation} \begin{aligned}
& \quad \exp\left( \boldsymbol{\mathbf{B}} t\right) \exp\left( \boldsymbol{\mathbf{C}} t\right) \\
& = \left(1+ \boldsymbol{\mathbf{B}} t+\frac{1}{2!} \boldsymbol{\mathbf{B}} ^2t^2 + \dots \right) \left(1+ \boldsymbol{\mathbf{C}} t+\frac{1}{2!} \boldsymbol{\mathbf{C}} ^2t^2 + \dots \right) \\
& = 1 + ( \boldsymbol{\mathbf{B}} + \boldsymbol{\mathbf{C}} )t + \frac{1}{2!} \left( \boldsymbol{\mathbf{B}} ^2 + \boldsymbol{\mathbf{C}} ^2 + 2 \boldsymbol{\mathbf{B}} \boldsymbol{\mathbf{C}} \right) t^2 + \dots\\
& = \exp\left( \boldsymbol{\mathbf{A}} t\right) + \mathcal{O}\left(t^2 \right) ~,
\end{aligned} \end{equation}
这里的 $ \mathcal{O}\left(t^2 \right) $ 是由于第二个等号后面是 $2 \boldsymbol{\mathbf{B}} \boldsymbol{\mathbf{C}} $ 而不是 $ \boldsymbol{\mathbf{B}} \boldsymbol{\mathbf{C}} + \boldsymbol{\mathbf{C}} \boldsymbol{\mathbf{B}} $。
1. ^ 例如二维波函数的动能算符 $T = T_x + T_y$