图

一阶线性常微分方程组

         

预备知识 一阶线性微分方程

解析解

   形式为

\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{M}} \right) = 1 + \boldsymbol{\mathbf{M}} + \frac1{2!} \boldsymbol{\mathbf{M}} ^2 + \dots \end{equation}
式 2 写成矩阵级数的形式即可验证式 1

形式解

   当式 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}} \right) \exp\left( \boldsymbol{\mathbf{Q}} \right) = \exp\left( \boldsymbol{\mathbf{P}} + \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}
然而这么做对于数值计算并没有太大意义.

   Expokit 是一个数学包可以用于计算矩阵的指数函数(Fortran 和 Matlab).

数值计算

   如果矩阵 $ \boldsymbol{\mathbf{M}} $ 是厄米矩阵, 则可以先做对角化 $ \boldsymbol{\mathbf{M}} = \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{U}} ^\dagger $, 其中 $ \boldsymbol{\mathbf{\Lambda}} $ 是对角矩阵, $ \boldsymbol{\mathbf{U}} $ 是酉矩阵. 这样就有

\begin{equation} \exp\left( \boldsymbol{\mathbf{M}} \right) = \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{U}} ^\dagger + \boldsymbol{\mathbf{U}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{U}} ^\dagger + \boldsymbol{\mathbf{U}} \frac1{2!} \boldsymbol{\mathbf{\Lambda}} ^2 \boldsymbol{\mathbf{U}} ^\dagger + \dots = \boldsymbol{\mathbf{U}} \exp\left( \boldsymbol{\mathbf{\Lambda}} \right) \boldsymbol{\mathbf{U}} ^\dagger \end{equation}
由于对角矩阵相乘等于每个对角元分别相乘, 把 $ \boldsymbol{\mathbf{\Lambda}} $ 的每个矩阵元求指数函数就可以得到 $ \exp\left( \boldsymbol{\mathbf{\Lambda}} \right) $. 这样做可以减少计算量.

   事实上, 以上做法相当于分离变量, 当 $ \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) $. 所以方程的通解为
\begin{equation} \boldsymbol{\mathbf{v}} (t) = \sum_{i = 1}^N c_i f_i(t) \boldsymbol{\mathbf{u}} _i = \sum_{i = 1}^N \boldsymbol{\mathbf{u}} _i ^\dagger \boldsymbol{\mathbf{v}} (0) \exp\left(\lambda_i t\right) \boldsymbol{\mathbf{u}} _i = \boldsymbol{\mathbf{U}} \exp\left(\Lambda t\right) \boldsymbol{\mathbf{U}} ^\dagger \boldsymbol{\mathbf{v}} (0) \end{equation}

   在此基础上使用 Lanczos 算法可以进一步提高效率.

Split Operator

   有时候我们希望可以在上述计算中把 $ \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$

致读者: 小时物理百科一直以来坚持所有内容免费无广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费, 甚至被收购。 因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们一个星期内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获得知识, 我们在此表示感谢。

         

© 小时物理百科 保留一切权利