一维谐振子定态 Schrödinger 方程的数值解法

Posted 2023-04-01

参考技术A 前几天整理电脑的时候发现了本科上量子力学讨论班时做的一个 Slide，觉得挺有意思的。花了点时间整理成这篇博客。

一个质量为 的粒子，在一维势场 中运动。其哈密顿算符为

其中 为位置算符， 为动量算符。我们需要求解该体系的定态 Schrödinger 方程：

一维谐振子是除了氢原子之外，为数不多的可以解析求解的体系。那么我们为什么要费劲求它的数值解呢？正因为绝大多数的量子体系都无法解析求解，数值方法才显得尤为重要。

回忆一下泰勒公式

令 ，有

两式相加，可得

将 在 区间离散化为

其中 ，则 Schrödinger 方程差分化为

在这里，我们假设 或 时 。这对能量较低的态是成立的。

将差分方程写成矩阵的形式为

这样一来，问题就转化为求差分矩阵的特征值和特征向量。

QR 算法是一种常见的特征值算法。它利用了矩阵的 QR 分解，即将矩阵 分解为一个正交矩阵 和一个上三角矩阵 的乘积：

为什么可以这么分解呢？我们回忆一下 Gram-Schmidt 正交化 ，将矩阵 的列向量看作一组基，则可以通过一系列初等列变换获得一组标准正交基。反过来看，对这组标准正交基所组成的矩阵 作初等列变换也可以得到矩阵 。我们知道，对矩阵作初等列变换相当于右乘初等矩阵，且由于 Gram-Schmidt 正交化不涉及列交换，这里用到的初等矩阵均为上三角矩阵。因此，矩阵 可以由正交矩阵 右乘一个上三角矩阵 得到。在实际应用中，除了 Gram-Schmidt 正交化，还可以用 Householder 变换 、 Givens 旋转 等方法实现 QR 分解。

那么如何利用 QR 分解求解矩阵的特征值呢？

记 ，对 ，

即在每一步，对 进行 QR 分解，再由分解后得到对 和 计算 ，如此迭代。

注意到 是正交矩阵，有 ，故

也就是说， 相似于 。根据递推关系， 全都是相似的，这意味着所有的 都有相同的特征值。在一定条件下， 会收敛为一个三角矩阵，特征值为其主对角元。

特别地，如果 是一个实对称正定矩阵（我们所要求解的差分矩阵刚好是这种情况）， 将会收敛为一个对角矩阵 ，且 依次递减。考虑到

当 收敛时， 的列向量即为属于相应特征值的特征向量。

综上，对于实对称正定矩阵 ，我们令 ，对 ，

当 收敛时，就同时求得 的特征值和特征向量。

我们可以直接调用 numpy.linalg.qr 作 QR 分解：

接下来就可以求解一维谐振子了。为了方便起见，我们令 。

一维谐振子的本征能量为：

对应的本征态为：

其中

为 Hermite 多项式。前三个 Hermite 多项式为：

我们同样令 ，则

画出来看看

八九不离十吧。低能级的误差主要来自截断误差和舍入误差。此外，高能级需要有更大的 来保证 或 时 的假设，因此能级越高，误差越大。