贡献者： addis
预备知识 氢原子球坐标薛定谔方程数值解
，
正交归一基底
1. 背景介绍
1在计算含时薛定谔方程时，若已知某时刻 $t$ 的波函数 $\Psi( \boldsymbol{\mathbf{r}} ,t)$，要求 $\Psi( \boldsymbol{\mathbf{r}} ,t+\Delta t)$，通常使用传播算符
\begin{equation}
U(\Delta t) = \exp\left(-\frac{ \mathrm{i} }{\hbar} \int_{t}^{t+\Delta t} \,\mathrm{d}{t} \boldsymbol\cdot H(t)\right)
\approx \exp\left(- \mathrm{i} H(t+\Delta t/2) \Delta t\right) ~,
\end{equation}
其中 $H$ 为哈密顿算符。若 $H$ 不含时，约等号变为等号
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} , t+\Delta t) = \exp\left(- \mathrm{i} H\Delta t/\hbar\right) \Psi( \boldsymbol{\mathbf{r}} , t)~.
\end{equation}
通常情况下，用有限个正交归一基底 $\Psi_0,\dots\Psi_{K-1}$ 近似展开 $\Psi( \boldsymbol{\mathbf{r}} , t+\Delta t)$，这时，$H$ 可以表示成矩阵的形式。
\begin{equation}
H_{ij} = \left\langle \Psi_i \right\rvert H \left\lvert \Psi_j \right\rangle ~.
\end{equation}
这样，幺正算符 $ \exp\left(- \mathrm{i} H\Delta t/\hbar\right) $ 也可以表示成矩阵。根据定义
\begin{equation}
\exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t/\hbar\right) = 1+ (- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t/\hbar) + \frac{1}{2!} (- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t/\hbar)^2\dots~
\end{equation}
此时若把 $ \boldsymbol{\mathbf{H}} $ 对角化（求解本征方程），得到本征矢列矩阵 $ \boldsymbol{\mathbf{P}} $（单位正交阵）以及本征值矩阵 $ \boldsymbol{\mathbf{\Lambda}} $（对角矩阵），则可进行基底变换变到 $H$ 矩阵的本征空间求解上式
\begin{equation} \begin{aligned}
\exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t/\hbar\right) &= \boldsymbol{\mathbf{P}} \ \ \mathrm{diag}( \mathrm{e} ^{- \mathrm{i} E_1 t/\hbar}, \mathrm{e} ^{- \mathrm{i} E_2 t/\hbar}\dots \mathrm{e} ^{- \mathrm{i} E_K t/\hbar}) \boldsymbol{\mathbf{P}} ^{-1}\\
&= \sum_{j=0}^{K-1} P_{ij} \mathrm{e} ^{- \mathrm{i} E_j t/\hbar} P_{jk}^{-1}
= \sum_{j=0}^{K-1} P_{ij}P_{kj} \mathrm{e} ^{- \mathrm{i} E_j t/\hbar} ~.
\end{aligned} \end{equation}
要说明的是，这种算法的误差（除数值计算误差）来源于两个方面，第一是用 $ \exp\left( \mathrm{i} H\Delta t/\hbar\right) $ 代替 $ \exp\left(-\frac{ \mathrm{i} }{\hbar} \int_{t}^{t+\Delta t} \,\mathrm{d}{t} \cdot H(t)\right) $，但如果不含时就没有该误差，第二是有限数量的基底不具有完备性，这个误差可以随着基底数量增加而减小。
2. Lanczos 算法
在以上的算法中，选取施密特正交归一化 的 Krilov 基底作为基底。令 $\Psi_0 \equiv \Psi( \boldsymbol{\mathbf{r}} , t)$， $K$ 阶的 Krilov 基底为
\begin{equation}
\left\{ \left\lvert \Psi_0 \right\rangle , H \left\lvert \Psi_0 \right\rangle , H^2 \left\lvert \Psi_0 \right\rangle \dots H^{K-1} \left\lvert \Psi_0 \right\rangle \right\} ~.
\end{equation}
我们把正交归一化后的基底记为
\begin{equation}
\left\{ \left\lvert \Psi_0 \right\rangle , \left\lvert \Psi_1 \right\rangle , \left\lvert \Psi_2 \right\rangle \dots \left\lvert \Psi_{K-1} \right\rangle \right\} ~.
\end{equation}
这样，用其展开哈密顿算符
\begin{equation}
H_{ij}= \left\langle \Psi_i \middle| H \middle| \Psi_j \right\rangle ~.
\end{equation}
在算法上，本来我们可以按部就班地按照以上步骤做，然而由于 Krilov 基底的性质，可以通过一些定理（见文章最后）大大减少计算量。正交归一化步骤如下
【1】如果 $\Psi( \boldsymbol{\mathbf{r}} ,t)$ 没有归一化，将其进行归一化
\begin{equation}
\left\lvert \Psi_0 \right\rangle =\Psi( \boldsymbol{\mathbf{r}} ,t)/\sqrt{ \left\langle{\Psi( \boldsymbol{\mathbf{r}} ,t)}\middle| \Psi( \boldsymbol{\mathbf{r}} ,t) \right\rangle }~.
\end{equation}
【2】把 $\Psi_1$ 进行施密特正交化，令
\begin{equation}
\lvert{\tilde\Psi_1}\rangle = (1 - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert )H \left\lvert \Psi_0 \right\rangle ~,
\qquad
\left\lvert \Psi_1 \right\rangle = \lvert{\tilde\Psi_1}\rangle /\sqrt{ \langle{\tilde\Psi_1}|{\tilde\Psi_1}\rangle }~.
\end{equation}
【3】现在起我们令 $\beta_j = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle }$。
定理 1
对以上定义的基底
\begin{equation}
\left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \ H^{j + 1} \left\lvert \Psi_0 \right\rangle ~
\end{equation}
与
\begin{equation}
\left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \ H \left\lvert \Psi_j \right\rangle ~.
\end{equation}
共线。
故现在起采用后者进行正交化。
定理 2
\begin{equation}
\left\langle \Psi_{j-1} \right\rvert H \left\lvert \Psi_j \right\rangle = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle } = \beta_j~.
\end{equation}
【4】把 $H^2 \left\lvert \Psi_0 \right\rangle $ 进行施密特正交化，令
\begin{equation}
\alpha_i = \left\langle \Psi_i \middle| H \middle| \Psi_i \right\rangle ~.
\end{equation}
根据以上两个定理，正交化结果也可以写成
\begin{equation} \begin{aligned}
\lvert{\tilde\Psi_2}\rangle &= \left(1 - \left\lvert \Psi_1 \right\rangle \left\langle \Psi_1 \right\rvert - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_1 \right\rangle \\
&= H \left\lvert \Psi_1 \right\rangle - \alpha_1 \left\lvert \Psi_1 \right\rangle - \beta_1 \left\lvert \Psi_0 \right\rangle ~.\\
\end{aligned} \end{equation}
【5】再次归一化 $ \left\lvert \Psi_2 \right\rangle = \lvert{\tilde\Psi_2}\rangle / \sqrt{ \langle{\tilde\Psi_2}|{\tilde\Psi_2}\rangle } = \lvert{\tilde\Psi_2}\rangle /\beta_2~.$
定理 3 定理 3
\begin{equation} \begin{aligned}
\lvert{\tilde\Psi_{j+1}}\rangle &= \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\
&= \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \left\lvert \Psi_{j-1} \right\rangle \left\langle \Psi_{j-1} \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\
&= H \left\lvert \Psi_j \right\rangle - \alpha_j \left\lvert \Psi_j \right\rangle - \beta_j \left\lvert \Psi_{j-1} \right\rangle ~,
\end{aligned} \end{equation}
这就是最简洁的正交化公式。比起最原始的正交化
\begin{equation}
\lvert{\tilde\Psi_{j+1}}\rangle = \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H^{j+1} \left\lvert \Psi_0 \right\rangle ~.
\end{equation}
定理 3 只用到了一个矩阵-矢量乘法，两个矢量数乘和两个矢量减法。
接下来只要对 $j = 2,4,5...K - 1$ 不断重复定理 3 中的正交化和归一化，就可以得到所有正交归一基底 $ \left\lvert \Psi_j \right\rangle $。
现在我们来求该基底下的 $ \boldsymbol{\mathbf{H}} $ 矩阵。根据 $\alpha_i$ 的定义以及定理 2 可以发现它们分别是矩阵的主对角元和副对角元。
定理 4
在正交归一基底 $ \left\lvert \Psi_j \right\rangle $ 中，矩阵 $ \boldsymbol{\mathbf{H}} $ 是对称三对角矩阵。所以我们已经顺便求出了所有的矩阵元！
\begin{equation}
\boldsymbol{\mathbf{H}} =
\begin{pmatrix}
\alpha_0 & \beta_1 &&& \\
\beta_1 & \alpha_1 & \beta_2 && \\
& \beta_2 & \alpha_2 & \ddots& \\
& & \ddots& \ddots & \beta_{K-1} \\
&&&\beta_{K-1} & \alpha_{K-1}
\end{pmatrix}~.\end{equation}
该矩阵具有维数小，易求本征问题的优势。令
\begin{equation}
\boldsymbol{\mathbf{H}} = \boldsymbol{\mathbf{P}} \boldsymbol{\mathbf{\Lambda}} \boldsymbol{\mathbf{P}} ^{-1}~,
\end{equation}
那么
\begin{equation}
\mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t} = \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) \boldsymbol{\mathbf{P}} ^\dagger ~.
\end{equation}
注意在 $ \left\lvert \Psi_i \right\rangle $ 基底下，$\Psi( \boldsymbol{\mathbf{r}} , t)$ 就是 $(1, 0, 0, \dots)$。所以在该基底下
\begin{equation}
\boldsymbol{\mathbf{a}} \equiv \mathrm{e} ^{- \mathrm{i} \boldsymbol{\mathbf{H}} \Delta t} \Psi( \boldsymbol{\mathbf{r}} , t)= \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) ( \boldsymbol{\mathbf{P}} ^\dagger )_{\text{第 1 列}} = \boldsymbol{\mathbf{P}} \exp\left(- \mathrm{i} \boldsymbol{\mathbf{\Lambda}} \Delta t\right) ( \boldsymbol{\mathbf{P}} _{\text{第 1 行}}) ^\dagger ~.
\end{equation}
这里的 $ \boldsymbol{\mathbf{a}} $ 其实就是 $ \exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \,\mathrm{d}{t} \right) $ 的第一列。
3. 误差分析
要判断 $K$ 是否足够，只需要计算式 21 中的矢量 $ \boldsymbol{\mathbf{a}} $ 的最后几个矩阵元的平方和是否可以忽略不计，如果可以就说明 $K$ 是足够的。
第一个被截去的系数 $a_{K+1}$ 的模长可以用下式估计（Feist 毕业论文 4.51）
\begin{equation}
\left\lvert a_{K+1} \right\rvert \approx \frac{\Delta t^K}{K!} \prod_{i=1}^K \beta_i~.
\end{equation}
推导：把
式 18 做 $K-1$ 次幂，左下角的矩阵元就是 $\prod_i \beta_i$。此时如果把
式 18 拓展到无限大的矩阵，那么做 $K-1$ 次幂后第一列也是相同的。也就是说，$ \exp\left(- \mathrm{i} \boldsymbol{\mathbf{H}} \,\mathrm{d}{t} \right) $ 的泰勒展开的前 $K$ 项的第一列
2都是精确的。若考察泰勒展开的 $K+1$ 项，则无限大矩阵的结果中第一列最后会多出一个元素 $\beta_1\beta_2\dots\beta_K$（可以借助 Mathematica）。接下来的推导就很简单了。
对于固定的 $K$，该式可以用于薛定谔方程的变步长演化。Feist 用的 Sort Iterative Lanczos 就是使用 $K$ 等于 "a few"，以及较小的 $\Delta t$。
4. 定理 1 证明
根据施密特正交化的性质，Krilov 基底的前 $j$ 项的与 $ \left\lvert \Psi_0 \right\rangle \dots \left\lvert \Psi_j \right\rangle $ 展开同一空间。所以 $H^j \left\lvert \Psi_0 \right\rangle = c_j \left\lvert \Psi_j \right\rangle +\ldots + c_0 \left\lvert \Psi_0 \right\rangle $，所以（以下所有系数 只是用来表示线性组合，具体值不重要）
\begin{equation} \begin{aligned}
H^{j+1} \left\lvert \Psi_0 \right\rangle &= c_j H \left\lvert \Psi_j \right\rangle + H \left(c_{j-1} \left\lvert \Psi_{j-1} \right\rangle \dots + c_0 \left\lvert \Psi_0 \right\rangle \right) \\
& = c_j H \left\lvert \Psi_j \right\rangle + H \left(c'_{j-1} H^{j-1} \left\lvert \Psi_0 \right\rangle + c'_{j-2} H^{j - 2} \left\lvert \Psi_0 \right\rangle \dots + c'_0 \left\lvert \Psi_0 \right\rangle \right) \\
&= c_j H \left\lvert \Psi_j \right\rangle + \left(c'_{j-1} H^j \left\lvert \Psi_0 \right\rangle + c'_{j-2} H^{j-1} \left\lvert \Psi_0 \right\rangle \dots + c'_0 \left\lvert \Psi_0 \right\rangle \right) \\
&= c_j H \left\lvert \Psi_j \right\rangle + \left(c''_{j-1} \left\lvert \Psi_j \right\rangle + c''_{j-2} \left\lvert \Psi_{j-1} \right\rangle + \ldots + c''_0 \left\lvert \Psi_0 \right\rangle \right) ~.\\
\end{aligned} \end{equation}
把 $H^{j+1} \left\lvert \Psi_0 \right\rangle $ 施密特正交归一化
\begin{equation} \begin{aligned}
&\phantom{=} \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H^{j+1} \left\lvert \Psi_0 \right\rangle \\
&= c_j \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle \\
&\phantom{=} + \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) \left(c''_{j-1} \left\lvert \Psi_j \right\rangle + c''_{j - 2} \left\lvert \Psi_{j - 1} \right\rangle + \ldots + c''_0 \left\lvert \Psi_0 \right\rangle \right) \\
&= c_j \left(1 - \left\lvert \Psi_j \right\rangle \left\langle \Psi_j \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H \left\lvert \Psi_j \right\rangle ~.
\end{aligned} \end{equation}
5. 定理 2 证明
要证 $ \left\langle \Psi_{j-1} \right\rvert H \left\lvert \Psi_j \right\rangle = \sqrt{ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle }$，即证 $ \langle{\Psi_{j - 1}}|{H}|{\tilde \Psi_j}\rangle = \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle $，即证 $ \langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{H\Psi_{j-1}}\rangle $ 而
\begin{equation}
\langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{ \left(1- \left\lvert \Psi_{j-1} \right\rangle \left\langle \Psi_{j-1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) H}|{\Psi_{j-1}}\rangle ~.
\end{equation}
由于上式中 $ \langle{{\tilde \Psi }_j}\rvert \left( \left\lvert \Psi_{j - 1} \right\rangle \left\langle \Psi_{j - 1} \right\rvert - \ldots - \left\lvert \Psi_0 \right\rangle \left\langle \Psi_0 \right\rvert \right) = 0~,$
\begin{equation}
\langle{\tilde\Psi_j}|{\tilde\Psi_j}\rangle = \langle{\tilde\Psi_j}|{H}|{\Psi_{j-1}}\rangle ~.
\end{equation}
6. 定理 3 证明
要证明定理 3，即证，对 $n \leqslant j - 2$，有 $ \left\langle \Psi_n \right\rvert H \left\lvert \Psi_j \right\rangle = 0$，$H$ 为厄米算符时即证 $ \left\langle \Psi_j \right\rvert H \left\lvert \Psi_n \right\rangle = 0$，即证对 $m \geqslant j + 2$， $ \left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = 0$。在定理 1 的证明类似，我们知道
\begin{equation} \begin{aligned}
H \left\lvert \Psi_j \right\rangle &= H \left(c_j H^j \left\lvert \Psi_0 \right\rangle + \ldots + c_0 \left\lvert \Psi_0 \right\rangle \right) = c_j H^{j+1} \left\lvert \Psi_0 \right\rangle + \ldots + c_0 H \left\lvert \Psi_0 \right\rangle \\
&= c'_j \left\lvert \Psi_{j+1} \right\rangle + \ldots + c'_0 \left\lvert \Psi_0 \right\rangle ~,\\
\end{aligned} \end{equation}
所以对 $m \geqslant j + 2$ 有
\begin{equation}
\left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = \left\langle \Psi_m \right\rvert \left(c'_j \left\lvert \Psi_{j+1} \right\rangle + \ldots + c'_0 \left\lvert \Psi_0 \right\rangle \right) = 0~.
\end{equation}
7. 定理 4 证明
首先证明 $ \boldsymbol{\mathbf{H}} $ 是三对角矩阵。即证 $m \geqslant j + 2$ 或 $m \leqslant j - 2$ 时 $ \left\langle \Psi_m \right\rvert H \left\lvert \Psi_j \right\rangle = 0$。
$m \geqslant j + 2$ 的情况在证明 3 中已经证明，只需证明另一种情况。对厄米矩阵，$ \langle{H\Psi_m}|{\Psi_j}\rangle = \left\langle \Psi_m \middle| H \middle| \Psi_j \right\rangle = 0$，取复共轭，即 $ \left\langle \Psi_j \middle| H \middle| \Psi_m \right\rangle = 0$。可见左边的角标的确小于等于右边的角标减 2，证毕。
然后证明 $ \boldsymbol{\mathbf{H}} $ 是实矩阵，考虑到厄米矩阵 $H_{mn} = H_{nm}^* $ 的性质，只需要证明 $H_{j, j+1}$ 是实数。对任意 $j = 2 \dots K-1$
\begin{equation} \begin{aligned}
H_{j, j+1} &= \left\langle \Psi_j \middle| H \middle| \Psi_{j+1} \right\rangle = \left\langle \Psi_j \right\rvert H \left(H \left\lvert \Psi_j \right\rangle - \alpha_j \left\lvert \Psi_j \right\rangle - \beta_j \left\lvert \Psi_{j-1} \right\rangle \right) \\
&= \left\langle \Psi_j \middle| H^2 \middle| \Psi_j \right\rangle - \alpha_j \left\langle \Psi_j \right\rvert H \left\lvert \Psi_j \right\rangle - \beta_j \left\langle \Psi_j \right\rvert H \left\lvert \Psi_{j-1} \right\rangle ~.
\end{aligned} \end{equation}
易证 $ \left\langle \Psi_j \middle| H^2 \middle| \Psi_j \right\rangle $，$\alpha_j$ 及 $\beta_j$ 是实数，只需要证明 $ \left\langle \Psi_j \middle| H \middle| \Psi_{j-1} \right\rangle $ 是实数即可，即证 $ \left\langle \Psi_{j-1} \middle| H \middle| \Psi_j \right\rangle $ 是实数，重复上述论证，即证 $ \left\langle \Psi_{j-2} \middle| H \middle| \Psi_{j - 1} \right\rangle $ 是实数，……，即证 $ \left\langle \Psi_1 \middle| H \middle| \Psi_0 \right\rangle $ 是实数，即证 $ \left\langle \Psi_0 \middle| H \middle| \Psi_1 \right\rangle $ 是实数，即证 $ \langle{\Psi_0}|{H}|{\tilde\Psi_1}\rangle $ 是实数。而
\begin{equation}
\langle{\Psi_0}|{H}|{\tilde\Psi_1}\rangle = \left\langle \Psi_0 \right\rvert H \left(H \left\lvert \Psi_0 \right\rangle - \alpha_0 \left\lvert \Psi_0 \right\rangle \right) = \left\langle \Psi_0 \middle| H^2 \middle| \Psi_0 \right\rangle - \alpha_0 \left\langle \Psi_0 \middle| H \middle| \Psi_0 \right\rangle ~.
\end{equation}
显然是实数。以上的思路是数学归纳法的逆过程。
1. ^ 参考 Renate Pazourek 的论文
2. ^ 为什么只看第一列？因为当前波函数就是第一个 Krylov 基底。