4. 打靶法计算任意一维势的束缚态能级

打靶法求解一维束缚态：从试探能量到能级匹配

打靶法把定态薛定谔方程当作初值问题处理：给定试探能量 E，从势阱一侧的物理边界出发，用数值积分把波函数推到另一侧，检查它是否满足另一侧的边界条件。对束缚态而言，波函数在两侧无穷远处都必须衰减到零；在有限计算区间 [-L, L] 上，这等价于要求外侧边界处的对数导数与衰减解匹配。能量偏离本征值时，波函数会在远端“发散”或“过零”，只有本征能量才让远端值趋于零。

本节采用无量纲形式，令 ħ = 1, m = 1，方程为 ψ'' = 2(V(x) - E)ψ。用四阶龙格–库塔法从 x = -L 出发，取 ψ(-L) = ε（小量）、ψ'(-L) = κ ε，其中 κ = sqrt(2(V(-L) - E)) 是外侧衰减指数，这样初始条件本身就满足左边界。积分到 x = +L 后，用右端失配函数 f(E) = ψ'(L) + κ ψ(L) 判断：f(E) = 0 即命中本征能量。由于 f(E) 随 E 单调变号，用二分法在区间内扫描即可稳定收敛。

实现：失配函数与二分扫描

下面给出可直接运行的完整实现。势函数用参数化双势阱 V(x) = a(x² - b²)² / b⁴，在 x = ±b 处有两个极小值，中间 x = 0 处为势垒，高度为 a。通过调节 a 可以观察分裂量 ΔE 的变化。

import numpy as np

def V(x, a=6.0, b=1.5):
    """双势阱：极小值在 ±b，势垒高度 a（无量纲）"""
    return a * (x**2 - b**2)**2 / b**4

def rhs(x, y, E, a, b):
    """y = [psi, psi']，返回导数"""
    psi, dpsi = y
    return np.array([dpsi, 2.0 * (V(x, a, b) - E) * psi])

def rk4_step(x, y, h, E, a, b):
    k1 = rhs(x, y, E, a, b)
    k2 = rhs(x + h/2, y + h/2*k1, E, a, b)
    k3 = rhs(x + h/2, y + h/2*k2, E, a, b)
    k4 = rhs(x + h, y + h*k3, E, a, b)
    return y + h/6.0 * (k1 + 2*k2 + 2*k3 + k4)

def mismatch(E, L=6.0, N=4000, a=6.0, b=1.5):
    """从 -L 积分到 +L，返回右端失配 f(E)"""
    x = -L
    h = 2.0 * L / N
    kappa = np.sqrt(max(2.0 * (V(-L, a, b) - E), 1e-12))
    y = np.array([1e-6, kappa * 1e-6])   # 满足左边界衰减
    for _ in range(N):
        y = rk4_step(x, y, h, E, a, b)
        x += h
    psi_R, dpsi_R = y
    return dpsi_R + kappa * psi_R

def find_levels(Emin, Emax, nscan=400, **kw):
    """扫描失配函数变号区间，二分求根"""
    Es = np.linspace(Emin, Emax, nscan)
    fs = np.array([mismatch(E, **kw) for E in Es])
    roots = []
    for i in range(len(Es) - 1):
        if fs[i] * fs[i+1] < 0:
            lo, hi = Es[i], Es[i+1]
            for _ in range(60):
                mid = 0.5 * (lo + hi)
                if mismatch(lo, **kw) * mismatch(mid, **kw) <= 0:
                    hi = mid
                else:
                    lo = mid
            roots.append(0.5 * (lo + hi))
    return roots

levels = find_levels(0.0, 6.0, a=6.0, b=1.5)
print("能级:", [f"{e:.6f}" for e in levels])
if len(levels) >= 2:
    print("分裂 ΔE =", f"{levels[1] - levels[0]:.6e}")

运行后应得到两个低能级，且 levels[1] - levels[0] 是一个小量。若把 a 从 6 增大到 12，势垒变高，分裂量会显著减小；把 b 增大使两阱远离，同样使 ΔE 减小。这正是双势阱中隧穿耦合导致能级分裂的数值体现。

波函数重建与宇称判读

求得本征能量后，用同一积分器保存整条轨迹即可得到波函数。双势阱的两个低能态分别具有偶宇称（无节点，两阱同号）与奇宇称（中间一个节点，两阱反号）。判读时统计区间内的过零次数：偶宇称态过零 0 次，奇宇称态过零 1 次，这是区分 E₁ 与 E₂ 最直接的数值证据。

def wavefunction(E, L=6.0, N=4000, a=6.0, b=1.5):
    x = -L
    h = 2.0 * L / N
    kappa = np.sqrt(max(2.0 * (V(-L, a, b) - E), 1e-12))
    y = np.array([1e-6, kappa * 1e-6])
    xs, psis = [x], [y[0]]
    for _ in range(N):
        y = rk4_step(x, y, h, E, a, b)
        x += h
        xs.append(x); psis.append(y[0])
    xs, psis = np.array(xs), np.array(psis)
    psis /= np.sqrt(np.trapz(psis**2, xs))   # 归一化
    return xs, psis

for E in levels[:2]:
    xs, psi = wavefunction(E)
    nodes = int(np.sum(np.sign(psi[:-1]) * np.sign(psi[1:]) < 0))
    print(f"E={E:.6f}  节点数={nodes}  宇称={'偶' if nodes % 2 == 0 else '奇'}")

收敛性与参数检查

打靶法的优势在于只需一维积分、无需构造矩阵，对任意形状的势都适用；代价是每次求根都要重新积分，且对深束缚态需要更长的区间与更细的步长。对双势阱这类近简并情形，关键是扫描密度足够高，确保两个相邻零点都被捕获，再用节点数确认宇称归属。