📚 课程目录第4章空间离散化方法

4. 空间离散化方法

在数值求解定态薛定谔方程时，连续的空间被替换为有限个离散格点，微分算符被替换为差分算符，最终问题转化为一个矩阵本征值问题。本节以最常用的一维有限差分法为主线，说明哈密顿矩阵的具体构造方式、边界条件的三种处理策略，以及如何通过网格收敛性分析判断数值结果是否可信。

一、从微分方程到矩阵本征值问题

一维定态薛定谔方程为

-ħ²/(2m) · d²ψ/dx² + V(x)ψ = Eψ

取自然单位 ħ = m = 1，在区间 [x₀, x_N] 上划分 N+1 个等距格点，步长 h = (x_N − x₀)/N，格点坐标 x_i = x₀ + i·h（i = 0, 1, …, N）。对二阶导数采用中心差分：

d²ψ/dx² |_{x_i} ≈ (ψ_{i+1} − 2ψ_i + ψ_{i−1}) / h²

代入方程并整理，得到格点 i 处的离散方程：

−(ψ_{i+1} − 2ψ_i + ψ_{i−1}) / (2h²) + V_i ψ_i = E ψ_i

把 ψ_i 的系数整理成矩阵形式 Hψ = Eψ，其中 H 为三对角对称矩阵：对角元 H_{ii} = 1/h² + V_i，次对角元 H_{i,i±1} = −1/(2h²)。注意系数 1/(2h²) 来自动能项 −(1/2)d²/dx²，若采用 ħ、m 不为 1 的单位制，需相应替换为 ħ²/(2mh²)。

以无限深方势阱（宽度 L = 1，V ≡ 0）为例，取 N = 200，则 h = 0.005，对角元约 40000，次对角元约 −20000。用稠密或稀疏本征值求解器（如 LAPACK 的 dsyevd、SciPy 的 eigh_tridiagonal）即可得到本征值与本征矢。解析解 E_n = n²π²/(2L²)，基态约 4.9348；N = 200 时数值基态通常在 4.935 附近，误差量级约 10⁻⁵，这正是有限差分法 O(h²) 精度的体现。

二、边界条件的三种处理策略

边界条件决定了矩阵首末行如何填写，是构造哈密顿矩阵时最容易出错的地方。常见策略有三类：

狄利克雷边界（固定端点）：要求 ψ(x₀) = ψ(x_N) = 0。此时未知量只取内部格点 i = 1, …, N−1，矩阵维度为 (N−1)×(N−1)，首末行自然不含越界项。无限深方势阱即属此类，实现最简单。
诺伊曼边界（零导数）：要求 ψ′(x₀) = ψ′(x_N) = 0。用虚拟格点法，令 ψ_{−1} = ψ_1，代入 i = 0 的差分方程，得到对角元变为 2/h² + V₀（动能系数翻倍），次对角元仍为 −1/(2h²)。物理上对应“壁外场对称延拓”，适用于对称势或反射边界。
周期边界：要求 ψ(x₀) = ψ(x_N)，矩阵在右上角与左下角各补一个 −1/(2h²) 的耦合项，成为循环三对角矩阵。适用于晶格、环形势等周期结构。

需要强调的是，边界条件不是“数值技巧”，而是物理模型的组成部分。同一势场配不同边界，本征谱可以完全不同。例如宽度为 L 的无限深势阱基态能量为 π²/(2L²)，而同样区间上的零势场加周期边界，基态能量为 0（常数解），二者不可混用。

三、网格收敛性分析

有限差分法的截断误差为 O(h²)，因此本征值误差也应随 h² 线性下降（在双对数坐标下斜率为 2）。验证方法是：固定物理区间，取一系列 N（如 50、100、200、400、800），分别计算目标能级，与解析解或高精度参考值比较，计算相对误差，再拟合误差随 h 的幂次。若拟合斜率明显偏离 2，通常意味着边界处理有误、势场在格点上采样不当，或求解器精度不足。

实际操作中还要注意两点。其一，势场若含尖峰或阶跃（如 δ 势、方势垒），差分法在间断附近的收敛阶会退化，需要局部加密或改用有限元、谱方法。其二，矩阵条件数随 N² 增长，N 过大时双精度求解可能损失有效位数，此时应优先使用针对三对角/稀疏结构的稳定算法，而非直接求逆。

四、收敛性验证的实操流程

下面给出一段可直接运行的 Python 示例，用无限深方势阱检验收敛阶。代码使用 SciPy 的三对角本征值求解器，避免构造稠密矩阵：

import numpy as np
from scipy.linalg import eigh_tridiagonal

L = 1.0
exact = np.pi**2 / (2 * L**2)   # 基态解析值 ≈ 4.934802

for N in [50, 100, 200, 400, 800]:
    h = L / N
    diag = np.full(N - 1, 1.0 / h**2)      # V=0
    off  = np.full(N - 2, -1.0 / (2 * h**2))
    E = eigh_tridiagonal(diag, off, select='i',
                         select_range=(0, 0))[0][0]
    err = abs(E - exact) / exact
    print(f"N={N:4d}  h={h:.5f}  E0={E:.8f}  rel_err={err:.3e}")

预期输出中，N 每翻一倍，相对误差约缩小为原来的 1/4，即误差 ∝ h²。若把边界条件误写成诺伊曼，基态能量会趋近于 0 而非 4.9348，误差不收敛，这是最典型的排障信号。另一个常见问题是势场采样：若 V(x) 在格点上取值而格点未覆盖势的极值点，等效于对势做了粗采样，收敛阶同样会退化。

判断网格是否足够细的实用准则：连续两次加密（N → 2N）后，目标能级的相对变化小于所需精度（如 10⁻⁶）时，可认为该能级已收敛；对高激发态，由于波函数振荡更快，所需 N 通常按量子数平方增长，应单独验证。

不要用“矩阵维度越大越准”作为唯一判断。若边界条件与物理模型不匹配，加密网格只会让错误结果收敛到错误的极限值。构造矩阵后，先检查首末行是否符合所选边界，再谈收敛性。

综上，有限差分法构造哈密顿矩阵的核心是：用三点模板写出三对角结构，用边界条件决定首末行，用 h² 收敛性检验结果可信度。三者缺一不可，任何一环出错都会使本征值失去物理意义。