一阶线性常微分方程组(简明微积分)

                     

贡献者: 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)$ 之间不对易,所以我们形式上定义一个时序算符子节 2 )$ \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$

                     

© 小时科技 保留一切权利