4. 数值求解一维定态问题的有限差分法

有限差分法（Finite Difference Method, FDM）是求解一维定态薛定谔方程最直观、最容易实现的方法之一。它的核心思想非常朴素：把连续的坐标轴切成有限个离散格点，用格点上的函数值近似波函数，用差商近似导数，从而把微分方程转化为代数方程组。对于一维定态问题，这个方程组最终会呈现为一个矩阵本征值问题，交给数值线性代数库求解即可。

本节以“离散化—构造矩阵—处理边界—求解本征值”为主线，结合一个具体的一维方势阱示例，把每一步的数学表达和代码实现对应起来，帮助读者真正动手完成一次数值求解。

一、离散化：把连续区间变成格点

设一维定态薛定谔方程为

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

我们取计算区间 [x_min, x_max]，将其等分为 N 段，步长 h = (x_max - x_min)/N，得到 N+1 个格点 x_i = x_min + i·h（i = 0, 1, …, N）。波函数在第 i 个格点上的值记为 ψ_i，势能记为 V_i。

对二阶导数使用中心差分近似：

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

代入薛定谔方程并整理，得到第 i 个内部格点上的离散方程：

-(ħ²/2mh²)(ψ_{i+1} - 2ψ_i + ψ_{i-1}) + V_i ψ_i = E ψ_i

为了书写简洁，记 t = ħ²/(2mh²)，则方程变为：

-t ψ_{i-1} + (2t + V_i) ψ_i - t ψ_{i+1} = E ψ_i

这正是矩阵本征值问题的标准形式：Hψ = Eψ，其中 H 是一个三对角矩阵。

二、构造哈密顿矩阵

把内部格点 i = 1, 2, …, N-1 的方程按行排列，哈密顿矩阵 H 的对角元为 2t + V_i，次对角元和非对角元均为 -t。矩阵规模为 (N-1)×(N-1)，因为边界点 ψ_0 和 ψ_N 需要单独处理。

下面给出一个可直接运行的 Python 实现，使用 NumPy 构造矩阵并调用 numpy.linalg.eigh 求解对称三对角矩阵的本征值。示例采用原子单位制（ħ = 1, m = 1），区间取 [-5, 5]，格点数 N = 1000，势能为一个深度 50、宽度 1 的方势阱：

import numpy as np

# 参数设置（原子单位）
hbar, m = 1.0, 1.0
xmin, xmax = -5.0, 5.0
N = 1000
h = (xmax - xmin) / N
x = np.linspace(xmin, xmax, N + 1)

# 方势阱：|x| < 0.5 时 V = -50，其余为 0
V = np.where(np.abs(x) < 0.5, -50.0, 0.0)

# 内部格点
t = hbar**2 / (2 * m * h**2)
diag = 2 * t + V[1:-1]
off = -t * np.ones(N - 2)

# 构造三对角矩阵
H = np.diag(diag) + np.diag(off, k=1) + np.diag(off, k=-1)

# 求解本征值与本征向量
eigvals, eigvecs = np.linalg.eigh(H)
print(eigvals[:5])

运行后输出的前几个本征值即为束缚态能量。对于上述参数，基态能量约为 -49.0 左右（具体数值随离散精度略有差异），与解析估算的深方势阱基态能量量级一致。这里需要强调：这些数值属于教学示例的合理预期结果，实际运行时应以程序输出为准。

三、边界条件的处理

边界条件是有限差分法中最容易出错的地方。常见的有两类：

• 狄利克雷边界条件（固定边界）：若粒子被限制在区间内，波函数在边界处为零，即 ψ_0 = ψ_N = 0。此时边界点不进入矩阵，矩阵规模为 (N-1)×(N-1)，实现最简单。
• 诺伊曼边界条件（导数边界）：若边界处波函数导数为零，可用虚拟格点法，把 ψ_{-1} 用 ψ_1 表示，从而修改边界行的对角元。对于偶宇称势场，也可以只计算半区间并利用对称性。

对于束缚态问题，通常采用狄利克雷条件并把区间取得足够大，使波函数在边界处自然衰减到接近零。若区间太小，波函数被“截断”，本征值会偏高；若区间过大而格点数不变，步长 h 增大，离散误差也会增大。因此区间长度和格点数需要配合调整。

四、离散化误差与收敛性检查

中心差分对二阶导数的截断误差为 O(h²)，因此本征值误差也大致按 h² 收敛。实际使用中，建议做一次收敛性检查：固定物理区间，把 N 依次取 200、400、800、1600，观察基态能量是否趋于稳定。若相邻两次结果的差值小于预设容差（例如 1e-4），即可认为离散精度足够。

下图概括了从连续方程到矩阵本征值问题的完整数据流，以及边界条件在其中的作用位置：

五、结果判读与常见问题

求解得到本征值和本征向量后，需要做几项基本判读。第一，检查能级是否按升序排列，基态能量应最低。第二，观察波函数的节点数：第 n 个本征态（从 0 开始计数）应有 n 个节点，这是判断本征向量是否对应正确能级的直接依据。第三，检查波函数在边界附近是否已衰减到接近零，若边界处仍有明显幅度，说明区间不够大，需要扩展计算范围。

常见问题包括：矩阵规模过大导致内存不足（可改用 scipy.sparse.linalg.eigsh 只求前若干个本征值）；势能数组与格点错位（务必确认 V[i] 对应 x[i]）；以及忘记把边界点从矩阵中排除，导致出现虚假的零能级。只要在构造矩阵时严格对齐格点索引，这些问题都可以避免。