C.2 分步算符法 · 量子力学阶梯

含时薛定谔方程的解 $\Psi(t)=\ee^{-\ii\hat Ht/\hbar}\Psi(0)$ 写起来一行，算起来麻烦： $\hat H=\hat T+\hat V$ 里动能 $\hat T$ 在位置表象是微分算符。 分步算符法（split-operator method）的观察是：

$\ee^{-\ii\hat V\Delta t/\hbar}$ 在位置表象是对角的——逐点乘一个相位，零成本；
$\ee^{-\ii\hat T\Delta t/\hbar}$ 在动量表象是对角的——也是逐点乘相位；
两个表象之间的换基恰好是傅里叶变换（A.2），而 FFT 把它做到了 $O(N\log N)$ 。

于是每一步演化 = 乘相位 → FFT → 乘相位 → 逆 FFT。正文里波包撞势垒、自由波包扩散的动画全是这个算法在跑。

分裂公式与误差阶

麻烦在于 $\hat T$ 与 $\hat V$ 不对易，所以 $\ee^{-\ii(\hat T+\hat V)\Delta t/\hbar}\ne\ee^{-\ii\hat T\Delta t/\hbar}\ee^{-\ii\hat V\Delta t/\hbar}$ 。但 $\Delta t$ 很小时可以近似拆开，拆法的对称性决定精度：

Trotter 与 Strang：为什么对称拆分多赚一阶进阶~6 min

记 $\hat A=-\ii\hat T\Delta t/\hbar$ 、 $\hat B=-\ii\hat V\Delta t/\hbar$ ，两者都是 $O(\Delta t)$ 的小量。 Baker–Campbell–Hausdorff 公式（把两个指数并成一个时，对易子会出来捣乱）：

\ee^{\hat A}\ee^{\hat B}=\exp\left(\hat A+\hat B+\tfrac12[\hat A,\hat B]+\cdots\right)\tag{C.2.1}

朴素拆分（Trotter）： $\ee^{\hat A}\ee^{\hat B}$ 与目标 $\ee^{\hat A+\hat B}$ 差在 $\tfrac12[\hat A,\hat B]=O(\Delta t^2)$ ——单步误差二阶，走完固定总时间要 $t/\Delta t$ 步，全局误差 $O(\Delta t)$ ：一阶方法。

对称拆分（Strang）：把势能劈成两半夹住动能，

\ee^{\hat B/2}\,\ee^{\hat A}\,\ee^{\hat B/2} =\exp\left(\hat A+\hat B+O(\Delta t^3)\right)\tag{C.2.2}

$O(\Delta t^2)$ 的对易子项被对称性抵消（左右互为镜像，奇数阶误差反号相消），单步误差三阶，全局误差 $O(\Delta t^2)$ ：二阶方法。额外代价几乎为零——相邻两步的「半步势能」还可以合并成一步。

一步演化的完整流程（Strang 分裂）：

\Psi\ \xrightarrow{\ \times\,\ee^{-\ii V\Delta t/2\hbar}\ } \ \xrightarrow{\ \text{FFT}\ } \ \xrightarrow{\ \times\,\ee^{-\ii\hbar k^2\Delta t/2m}\ } \ \xrightarrow{\ \text{FFT}^{-1}\ } \ \xrightarrow{\ \times\,\ee^{-\ii V\Delta t/2\hbar}\ }\ \Psi'\tag{C.2.3}

完整代码：波包撞方势垒

$E\approx V_0$ 的波包打在势垒上，一半穿过一半弹回（自然单位 $\hbar=m=1$ ）：

import numpy as np

# ---- 网格 ----
L, N = 40.0, 1024
x = np.linspace(-L/2, L/2, N, endpoint=False)
dx = x[1] - x[0]
k = 2 * np.pi * np.fft.fftfreq(N, d=dx)     # FFT 约定下的动量网格

# ---- 势能：中间一道方势垒 ----
V0, a = 2.0, 1.0
V = np.where(np.abs(x) < a / 2, V0, 0.0)

# ---- 初态：从左边入射的高斯波包（中心动量 k0，动能 k0²/2 = 2 = V0）----
x0, k0, sigma = -10.0, 2.0, 1.5
psi = np.exp(-(x - x0)**2 / (4 * sigma**2) + 1j * k0 * x)
psi /= np.sqrt(np.sum(np.abs(psi)**2) * dx)   # 归一化

# ---- 预先算好演化因子（Strang 分裂：V/2 → T → V/2）----
dt, steps = 0.005, 2000
expV_half = np.exp(-0.5j * V * dt)
expT      = np.exp(-0.5j * k**2 * dt)

for _ in range(steps):
    psi = expV_half * psi                        # 半步势能
    psi = np.fft.ifft(expT * np.fft.fft(psi))    # 整步动能（动量空间）
    psi = expV_half * psi                        # 半步势能

# ---- 检验：模守恒 + 透射/反射概率 ----
norm = np.sum(np.abs(psi)**2) * dx
T = np.sum(np.abs(psi[x >  a/2])**2) * dx
R = np.sum(np.abs(psi[x < -a/2])**2) * dx
print(f"演化 {steps} 步后  norm = {norm:.12f}")
print(f"透射 T = {T:.4f}   反射 R = {R:.4f}   T + R = {T+R:.4f}")

实测输出：

演化 2000 步后  norm = 1.000000000000
透射 T = 0.5119   反射 R = 0.4865   T + R = 0.9984

关键行说明：

np.fft.fftfreq(N, d=dx)：给出 FFT 约定下每个数组下标对应的波数（前一半是正 $k$ 、后一半是负 $k$ ），乘 $2\pi$ 后就是动量网格。 自己手排 $k$ 数组是本算法头号错误来源——顺序排错，波包会瞬间碎掉。
endpoint=False：FFT 隐含周期边界， $x=-L/2$ 与 $x=+L/2$ 是同一个点，不能重复取。
演化因子在循环外预先算好：每步只剩逐点乘法和两次 FFT。
norm 精确到 12 位——保模是结构性的，不依赖 $\Delta t$ 小。
$T+R=0.9984$ 不是误差：还差的 0.0016 是仍滞留在势垒区 $\lvert x\rvert<a/2$ 内的概率。
中途每隔若干步把 np.abs(psi)**2 画出来，就是正文里的隧穿动画。

收敛阶检验

在谐振子势里演化一个相干态，用极小步长的解当参考，看误差随 $\Delta t$ 的变化：

import numpy as np

L, N = 40.0, 1024
x = np.linspace(-L/2, L/2, N, endpoint=False)
dx = x[1] - x[0]
k = 2 * np.pi * np.fft.fftfreq(N, d=dx)
V = 0.5 * x**2                        # 谐振子势

psi0 = np.exp(-(x - 3.0)**2 / 2)      # 相干态：来回摆而不变形
psi0 = psi0 / np.sqrt(np.sum(np.abs(psi0)**2) * dx)

def evolve(psi, dt, steps):
    expV_half = np.exp(-0.5j * V * dt)
    expT = np.exp(-0.5j * k**2 * dt)
    for _ in range(steps):
        psi = expV_half * np.fft.ifft(expT * np.fft.fft(expV_half * psi))
    return psi

T_total = 2.0
ref = evolve(psi0.copy(), T_total / 6400, 6400)   # 极小步长作参考解

print("   dt        误差 ‖ψ - ψ_ref‖    误差比")
prev = None
for steps in [50, 100, 200, 400]:
    psi = evolve(psi0.copy(), T_total / steps, steps)
    err = np.sqrt(np.sum(np.abs(psi - ref)**2) * dx)
    ratio = f"{prev/err:.2f}" if prev else "  --"
    print(f" {T_total/steps:.4f}     {err:.3e}      {ratio}")
    prev = err

实测输出：

   dt        误差 ‖ψ - ψ_ref‖    误差比
 0.0400     1.162e-03        --
 0.0200     2.905e-04      4.00
 0.0100     7.257e-05      4.00
 0.0050     1.809e-05      4.01

步长减半、误差准确变为 1/4——Strang 分裂的 $O(\Delta t^2)$ 全局误差，和推导完全一致。

三个陷阱：都来自「FFT 隐含周期边界」

波包会从右边穿出、左边钻回。 盒子是一个圆环。演化太久，透射波会绕回来污染反射波。对策：盒子取大、及时停表，或在两端加吸收层——一段虚势能 $-\ii\Gamma(x)$ （它故意破坏保模，把到边的概率「吃掉」，模拟粒子一去不返）。
动量也有上限。 网格能表示的最大动量是 $k_{\max}=\pi/dx$ 。波包动量接近它时发生混叠（aliasing）——高动量成分伪装成反向低动量，图像出现莫名其妙的碎波。对策：保证 $\lvert k_0\rvert+3/\sigma\ll\pi/dx$ ，或加密网格。
势能跳变处的假振铃。 方势垒的锐利边缘含无限高的傅里叶分量，网格截断后在边缘附近留下小振荡。物理结论（T、R）通常不受影响，但要画波函数细节时可把势垒边缘平滑化。

本节关键公式

Strang 分裂
$\ee^{-\ii\hat H\Delta t/\hbar}\approx\ee^{-\ii\hat V\Delta t/2\hbar}\,\ee^{-\ii\hat T\Delta t/\hbar}\,\ee^{-\ii\hat V\Delta t/2\hbar}$
对称拆分；全局误差 O(Δt²)，朴素拆分只有 O(Δt)
误差来源
$\ee^{\hat A}\ee^{\hat B}=\ee^{\hat A+\hat B+\frac12[\hat A,\hat B]+\cdots}$
BCH 公式：T 与 V 不对易才需要拆分技巧
每步流程
$\Psi\to\ee^{-\ii V\Delta t/2\hbar}\Psi\to\mathcal{F}\to\ee^{-\ii\hbar k^2\Delta t/2m}\to\mathcal{F}^{-1}\to\ee^{-\ii V\Delta t/2\hbar}$
两次 FFT + 三次逐点乘相位；严格保模
网格动量上限
$k_{\max}=\pi/\dd x$
接近上限即混叠；周期边界 = 盒子是圆环

分步算符法

建议先掌握

学完本节你应该能

分裂公式与误差阶

完整代码：波包撞方势垒

收敛阶检验

本节关键公式

自测共 3 题

分裂公式与误差阶

完整代码：波包撞方势垒

收敛阶检验

∎本节关键公式

本节关键公式