3. 有限差分法求解一维定态薛定谔方程

问题设定与离散化思路

一维定态薛定谔方程在原子单位（ℏ=1, m=1）下写作

−½·d²ψ/dx² + V(x)ψ = Eψ

本节取无限深方势阱作为验证对象：宽度 L=1，势阱内 V=0，阱外 V=∞，边界条件 ψ(0)=ψ(L)=0。解析能级为 Eₙ = n²π²/(2L²)，取 L=1 时 E₁≈4.9348，E₂≈19.7392，E₃≈44.4132（原子单位）。

把区间 [0, L] 划成 N 个内部格点，步长 h = L/(N+1)，格点 xᵢ = i·h（i=1…N）。用中心差分近似二阶导数：

d²ψ/dx² ≈ (ψᵢ₊₁ − 2ψᵢ + ψᵢ₋₁)/h²

代入方程并整理，得到每个内部格点上的离散方程：

−(1/(2h²))·ψᵢ₋₁ + (1/h² + Vᵢ)·ψᵢ − (1/(2h²))·ψᵢ₊₁ = E·ψᵢ

这正是矩阵本征值问题 H·ψ = E·ψ，其中 H 是三对角对称矩阵：主对角线为 1/h² + Vᵢ，上下次对角线恒为 −1/(2h²)。边界格点 ψ₀=ψ_{N+1}=0 不进入未知量，因此矩阵规模为 N×N。由于 H 对称三对角，用 scipy.linalg.eigh_tridiagonal 求解既快又稳。

Python 实现：构造三对角矩阵并求解

下面代码只依赖 NumPy 与 SciPy，直接构造主对角与次对角数组，调用 eigh_tridiagonal 得到升序本征值与本征向量。波函数按最大绝对值归一化，便于与解析解比较形状。

import numpy as np
from scipy.linalg import eigh_tridiagonal

def solve_well(N, L=1.0):
    """无限深方势阱：返回本征值(原子单位)与本征向量矩阵(列为本征态)"""
    h = L / (N + 1)
    diag = np.full(N, 1.0 / h**2)          # 主对角 a = 1/h^2 + V, 阱内 V=0
    off = np.full(N - 1, -1.0 / (2.0 * h**2))  # 次对角 b
    E, vec = eigh_tridiagonal(diag, off)
    # 归一化：按最大绝对值缩放，便于比较形状
    vec = vec / np.max(np.abs(vec), axis=0)
    return E, vec, h

# 解析能级
def exact_levels(n_max, L=1.0):
    n = np.arange(1, n_max + 1)
    return n**2 * np.pi**2 / (2 * L**2)

for N in [19, 39, 79, 159, 319]:
    E, vec, h = solve_well(N)
    E_ex = exact_levels(3)
    err = np.abs(E[:3] - E_ex) / E_ex
    print(f"N={N:4d}  h={h:.5f}  E1={E[0]:.6f}  "
          f"相对误差 E1/E2/E3 = {err[0]:.2e} / {err[1]:.2e} / {err[2]:.2e}")

运行后应观察到：随 N 增大（h 减小），E₁、E₂、E₃ 均从上方单调逼近解析值，相对误差按 O(h²) 下降。这是中心差分二阶精度的直接体现。若把 N 翻倍，误差约缩小为原来的 1/4。

收敛性验证与结果判读

为了定量确认二阶收敛，可对误差取对数，用最小二乘拟合斜率：

Ns = np.array([19, 39, 79, 159, 319, 639])
hs, errs = [], []
for N in Ns:
    E, _, h = solve_well(N)
    hs.append(h)
    errs.append(abs(E[0] - exact_levels(1)[0]) / exact_levels(1)[0])

slope, _ = np.polyfit(np.log(hs), np.log(errs), 1)
print("拟合斜率 ≈", slope)   # 预期接近 2

判读要点：

• 斜率接近 2 说明离散格式为二阶收敛；若明显偏离，检查是否误用了单侧差分或边界处理有误。
• 基态波函数应呈单个半正弦，第一激发态一个节点，第二激发态两个节点；节点数等于 n−1，可作为本征态排序正确性的快速检查。
• 能量误差始终为正（离散能级偏高），这是变分性质的体现，可作为符号方向的辅助判断。