1. 一维谐振子波包演化：从解析解到数值验证

问题设定与无量纲化

一维谐振子哈密顿量为 H = p²/2m + ½mω²x²。为便于数值计算，取 m = 1、ħ = 1、ω = 1，则 H = -½∂²/∂x² + ½x²。初始态取偏离中心的相干态（即被平移的基态），其解析演化严格保持高斯形状，中心按 x_c(t) = x₀cos(t) 振荡，宽度不变。这正好用来检验数值演化是否正确。

初始波函数取 ψ(x,0) = (1/π)^(1/4) exp(-(x-x₀)²/2)，取 x₀ = 2.0。解析解为 ψ(x,t) = (1/π)^(1/4) exp(-(x-x_c(t))²/2) · exp(i·p_c(t)·x) · phase(t)，其中 p_c(t) = -x₀sin(t)。概率密度 |ψ|² 是中心在 x_c(t) 处、方差恒为 1/2 的高斯，周期 T = 2π。

分裂算符法：为什么适合本问题

含时薛定谔方程 i∂ψ/∂t = Hψ 的形式解为 ψ(t+Δt) = exp(-iHΔt)ψ(t)。分裂算符法利用 Strang 分裂，把 exp(-iHΔt) 近似为：

exp(-iHΔt) ≈ exp(-iVΔt/2) · exp(-iTΔt) · exp(-iVΔt/2) + O(Δt³)

其中动能部分 exp(-iTΔt) 在动量空间是对角算符（乘以相位 exp(-i k²Δt/2)），势能部分在位置空间是对角算符（乘以相位 exp(-i V(x)Δt/2)）。因此每一步只需两次 FFT：位置→动量、乘动能相位、动量→位置、乘势能相位。误差为 O(Δt³) 每步，全局 O(Δt²)，且保持幺正性，概率守恒。

完整可运行实现

下面代码自包含：构造网格、初始波函数、分裂算符演化循环，并在若干时刻与解析解比较中心位置、宽度和概率归一化。使用 numpy.fft，无需额外依赖。

import numpy as np

# ---------- 参数 ----------
N = 1024                 # 网格点数
L = 20.0                 # 空间范围 [-L/2, L/2]
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 顺序）

x0 = 2.0                 # 初始相干态中心
dt = 0.01                # 时间步长
n_steps = 628            # 约 2 个周期（T=2π≈6.28）

# ---------- 初始波函数 ----------
psi = (1/np.pi)**0.25 * np.exp(-(x - x0)**2 / 2).astype(complex)

# ---------- 预计算相位因子 ----------
V = 0.5 * x**2
expV_half = np.exp(-1j * V * dt / 2)          # 势能半步相位
expT = np.exp(-1j * k**2 * dt / 2)            # 动能整步相位

# ---------- 解析解 ----------
def analytic(x, t, x0=2.0):
    xc = x0 * np.cos(t)
    pc = -x0 * np.sin(t)
    # 相干态解析波函数（含整体相位，概率密度只用模方）
    norm = (1/np.pi)**0.25
    psi_a = norm * np.exp(-(x - xc)**2 / 2) * np.exp(1j * pc * x)
    return psi_a

# ---------- 演化循环 ----------
def split_step(psi, n_steps, dt):
    psi = psi.copy()
    for _ in range(n_steps):
        psi *= expV_half                       # 势能半步
        psi_k = np.fft.fft(psi)
        psi_k *= expT                          # 动能整步
        psi = np.fft.ifft(psi_k)
        psi *= expV_half                       # 势能半步
    return psi

# ---------- 运行并验证 ----------
print(f"{'t':>6} {'<x>num':>9} {'<x>ana':>9} {'width_num':>10} {'width_ana':>10} {'norm':>8}")
for step in [0, 157, 314, 471, 628]:
    psi_t = split_step(psi, step, dt)
    t = step * dt
    prob = np.abs(psi_t)**2
    norm = np.sum(prob) * dx
    x_mean = np.sum(x * prob) * dx
    x2_mean = np.sum(x**2 * prob) * dx
    width = np.sqrt(x2_mean - x_mean**2)

    psi_a = analytic(x, t)
    prob_a = np.abs(psi_a)**2
    x_mean_a = np.sum(x * prob_a) * dx
    x2_mean_a = np.sum(x**2 * prob_a) * dx
    width_a = np.sqrt(x2_mean_a - x_mean_a**2)

    print(f"{t:6.2f} {x_mean:9.4f} {x_mean_a:9.4f} {width:10.4f} {width_a:10.4f} {norm:8.5f}")

结果判读与收敛性检查

预期输出（示例，实际数值随参数略有差异）：

判读要点：

• 中心轨迹：⟨x⟩ 在 ±2 之间按余弦振荡，与 x₀cos(t) 一致，说明势能相位与动能相位配合正确。
• 宽度不变：数值宽度稳定在 1/√2 ≈ 0.7071，这是谐振子相干态的关键特征——与自由粒子波包会扩散形成鲜明对照。
• 概率守恒：归一化始终为 1.00000，验证分裂算符的幺正性。
• 收敛性：把 dt 从 0.01 减到 0.005，同一时刻的 ⟨x⟩ 误差应约缩小 4 倍（O(Δt²)）。若误差不降，检查网格范围是否足够（波包不能碰到边界）或 FFT 动量网格是否用 fftfreq 正确构造。