数值求解薛定谔方程

上一节末尾把话挑明了：真实的势是光滑的，分段常数的台阶逼近不是长久之计。而回望第 2 章，我们能精确解出的势屈指可数——无限阱、有限阱、谐振子、方垒，每个都靠一件专属的巧劲（驻波条件、超越方程、升降算符）。换一个稍微不规则的势，比如两个不等深的阱、或者加了电场的斜底阱，所有巧劲同时失灵。

但第 3 章早就把出路写好了，只是当时听着像哲学：波函数是矢量的分量，算符是矩阵（3.1 节、3.4 节）。那里的「矩阵」是无穷维的，写不进计算机。这一节做一件不客气的事：把无穷维砍成有限维——$x$ 轴换成一排格点，$\psi(x)$ 换成一列数，$\hat H$ 就真的变成一个你能打印出来看的矩阵。能级？矩阵的本征值。波函数？本征矢量。剩下的交给线性代数库。

第一步：把函数变成数组

选一个足够大的区间 $[-L/2,\,L/2]$（大到束缚态波函数在端点处早已衰减殆尽），撒 $N$ 个等距格点：

$$x_n=-\frac{L}{2}+nh,\qquad h=\frac{L}{N+1},\qquad n=1,2,\dots,N\tag{4.4.1}$$

波函数从连续曲线降格为一列采样值 $\psi_n\equiv\psi(x_n)$。势能同样变成一列数 $V_n\equiv V(x_n)$——这就是为什么任何势都不再特殊：势再古怪，也不过是数组里换一列数字。

定态薛定谔方程里动能项含 $\psi''$。麻烦来了：只剩一排离散的点，导数怎么办？

第二步：用邻居算二阶导

推导：三点差分公式及其误差基础~7 min展开

思路：一个点的二阶导数衡量曲线在该点的弯曲程度，而弯曲程度可以用「两侧邻居的平均值比本点高多少」来探测。把这句话变成公式，工具是泰勒展开。

**第一步：向两边展开。**假设 $\psi$ 足够光滑，在 $x_n$ 附近展开到四阶：

$$\psi_{n\pm1}=\psi(x_n\pm h)=\psi_n\pm h\psi_n'+\frac{h^2}{2}\psi_n''\pm\frac{h^3}{6}\psi_n'''+\frac{h^4}{24}\psi_n''''+\cdots\tag{4.4.2}$$

第二步：相加。两式相加时，所有奇数阶项（$\pm h\psi'$、$\pm h^3\psi'''/6$）符号相反，成对抵消——这就是取「左右对称的邻居」的全部动机：

$$\psi_{n+1}+\psi_{n-1}=2\psi_n+h^2\psi_n''+\frac{h^4}{12}\psi_n''''+\cdots\tag{4.4.3}$$

**第三步：解出 $\psi''$。**移项、除以 $h^2$：

$$\psi_n''=\frac{\psi_{n+1}-2\psi_n+\psi_{n-1}}{h^2}-\frac{h^2}{12}\psi_n''''+\cdots\tag{4.4.4}$$

丢掉修正项，就得到三点差分公式：

$$\boxed{\ \psi''(x_n)\approx\frac{\psi_{n+1}-2\psi_n+\psi_{n-1}}{h^2}\ }\tag{4.4.5}$$

误差是 $O(h^2)$：主导误差项 $-\frac{h^2}{12}\psi''''$ 随格距的平方缩小。把 $h$ 减半，误差缩到四分之一——这个「减半验四」的指纹，等会儿在数值实验里要亲眼核对。

读一眼公式的模样：分子正是「左右邻居之和减去本点两倍」，即邻居均值与本点之差的两倍——和开头的直觉严丝合缝。$\psi$ 在局部是直线时分子为零：直线不弯曲，二阶导为零，正确。

第三步：哈密顿量现出矩阵原形

把差分公式塞进定态方程 $-\frac{\hbar^2}{2m}\psi''+V\psi=E\psi$，第 $n$ 个格点上：

$$-\frac{\hbar^2}{2mh^2}\psi_{n+1}+\left(\frac{\hbar^2}{mh^2}+V_n\right)\psi_n-\frac{\hbar^2}{2mh^2}\psi_{n-1}=E\,\psi_n\tag{4.4.6}$$

左边是对列向量 $(\psi_1,\dots,\psi_N)$ 的线性操作——它就是矩阵乘法 $H\boldsymbol\psi=E\boldsymbol\psi$，其中

$$H=\begin{pmatrix} \frac{\hbar^2}{mh^2}+V_1&-\frac{\hbar^2}{2mh^2}&&\\[2pt] -\frac{\hbar^2}{2mh^2}&\frac{\hbar^2}{mh^2}+V_2&-\frac{\hbar^2}{2mh^2}&\\[2pt] &\ddots&\ddots&\ddots \end{pmatrix}\tag{4.4.7}$$

一个三对角矩阵：对角线装着「动能基底 + 当地势能」，紧邻对角线装着格点间的耦合，其余全为零——因为三点公式只跟左右邻居说话。

两处细节交代清楚：

• 边界条件：$n=1$ 的方程里出现 $\psi_0$、$n=N$ 里出现 $\psi_{N+1}$，即区间端点的值。我们直接令 $\psi_0=\psi_{N+1}=0$——相当于把系统关进一个宽 $L$ 的无限深盒子。只要 $L$ 够大、束缚态在墙边早就指数衰减到零，这两堵假墙就无关痛痒（多大才「够大」，见下面的警告）。
• 矩阵是实对称的——第 3 章厄米算符的离散化身。3.5 节的定理原样兑现：本征值全是实数，本征矢量相互正交。理论保证与数值验证在下面的代码里当场对账。

第四步：跑起来

用谐振子练手（取无量纲单位 $\hbar=m=\omega=1$，能级应为 $n+\frac12$）。三对角矩阵有专用的高效例程 eigh_tridiagonal，不必存储 $N^2$ 个零：

import numpy as np
from scipy.linalg import eigh_tridiagonal

def solve(V, L, N):
    """解 -psi''/2 + V(x) psi = E psi，区间 [-L/2, L/2]，N 个内部格点。"""
    x = np.linspace(-L/2, L/2, N + 2)[1:-1]     # 去掉两端的边界点
    h = x[1] - x[0]
    diag = 1.0 / h**2 + V(x)                    # 对角元：动能基底 + 势能
    off  = -0.5 / h**2 * np.ones(N - 1)         # 次对角元：邻居耦合
    E, psi = eigh_tridiagonal(diag, off, select='i', select_range=(0, 5))
    return x, E, psi / np.sqrt(h)               # 除以 √h：让 Σ|ψ|²h = 1

x, E, psi = solve(lambda x: 0.5 * x**2, L=16, N=1000)
print(E)        # [0.499992  1.499960  2.499896  3.499800  4.499673  5.499513]

六个能级与精确值 $0.5,1.5,\dots,5.5$ 的偏差都在小数点后第四、五位。顺手验证第 3 章的承诺：psi[:,0] @ psi[:,1] * h 输出 $0.0$——不同本征态严格正交，机器精度以内。

第五步：怀疑你的结果

数值计算最重要的纪律：结果必须证明自己收敛了。上面的数字有两个独立的误差来源，各配一个检查。

**误差一：格距 $h$ 有限。**加密网格，看结果怎么变：

$N$	50	100	200	400	800
$E_0$ 的误差	$3.1\times10^{-3}$	$7.9\times10^{-4}$	$2.0\times10^{-4}$	$5.0\times10^{-5}$	$1.2\times10^{-5}$

$N$ 每翻倍（$h$ 减半），误差几乎精确地缩到 $1/4$——推导预言的 $O(h^2)$ 指纹清清楚楚。这条纪律的实操版本：把 $N$ 翻倍重算一遍，头几位数字不动了，才可以引用它们。

**误差二：盒子 $L$ 有限。**假墙离得太近会把波函数「夹疼」。检查方式相同：把 $L$ 加大重算，结果不变才算数。

∎本节关键公式

• 三点差分

  $$\psi''(x_n)\approx\frac{\psi_{n+1}-2\psi_n+\psi_{n-1}}{h^2},\quad \text{误差}\ -\frac{h^2}{12}\psi''''$$

  奇数阶项对称相消；h 减半误差缩 1/4

• 离散哈密顿量

  $$H_{nn}=\frac{\hbar^2}{mh^2}+V_n,\qquad H_{n,n\pm1}=-\frac{\hbar^2}{2mh^2}$$

  实对称三对角矩阵；本征值实、本征矢正交

• 离散阱精确谱

  $$E_j^{(h)}=\frac{\hbar^2}{mh^2}\left[1-\cos\frac{j\pi h}{L}\right]\approx E_j\left[1-\frac{(j\pi h/L)^2}{12}\right]$$

  差分法把能级算低；高能级误差按 j² 增长

• 离散内积

  $$\braket{\phi}{\psi}\to\sum_n\phi_n^*\psi_n\,h$$

  归一化与正交性检查都要带上测度 h

?自测共 4 题

1. 1.

   三点差分公式的推导中，把 ψ(x+h) 与 ψ(x−h) 的展开式相加的目的是：

   • 让所有奇数阶导数项成对抵消，直接暴露二阶导数
   • 提高数值稳定性
   • 保证矩阵是三对角的
   • 消去势能项

   提交
2. 2.

   离散哈密顿矩阵是实对称的。这在数值上兑现了第 3 章的哪些承诺？（多选）

   多选题

   • 本征值全部是实数
   • 不同本征值的本征矢量相互正交
   • 本征矢量可以构成完备基（N 维空间的 N 个正交矢量）
   • 本征值一定非负

   提交
3. 3.

   用 N = 400 算出的基态能量误差是 5.0×10⁻⁵。其他条件不变，把 N 提到 800，误差大约变成：

   • 1.2×10⁻⁵（约 1/4）
   • 2.5×10⁻⁵（约 1/2）
   • 5.0×10⁻⁵（不变）
   • 6×10⁻⁶（约 1/8）

   提交
4. 4.

   宽 L 的无限深势阱用格距 h = L/10 的差分法求解。用离散精确谱 E₁⁽ʰ⁾/E₁ ≈ 1 − (πh/L)²/12 估算基态能量的相对误差，以百分比作答。

   %允许 10% 相对误差

   提交

接下来

方法齐了，代码也有了，但每次换势都要重新摆弄网格、单位、检查项，像每顿饭都从生火开始。下一节把这些打包成一个通用求解器：喂给它任意一条势能曲线，它吐出能级和波函数，并自带质量检查清单。然后用它做一件解析方法办不到的事——解双势阱，看着一对能级挤到只差百分之一，从中读出隧穿、氨分子微波激射器、以及化学键的雏形。