一阶线性常微分方程组

             

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

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

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}} \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}
然而这么做对于数值计算并没有太大意义.

3. 厄米矩阵

   如果 $ \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$ 将不随时间变化.

   证明:把式 3 两边取厄米共轭得

\begin{equation} ( \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{H}} t}) ^\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{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 算法可以进一步提高效率.

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$

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

         

© 小时科技 保留一切权利