3. 薛定谔方程数值求解基础

含时薛定谔方程是量子动力学仿真的出发点。对二能级系统，它可写成两个耦合的一阶复微分方程，形式简单却包含了量子演化最核心的相位与跃迁机制。本节以二能级系统为对象，把“连续方程”变成“可执行的离散迭代”，并比较欧拉法与Crank-Nicolson格式在稳定性与精度上的差异，最后用拉比振荡作为验证场景。

一、从物理模型到可计算形式

设二能级系统的基态为 |0⟩、激发态为 |1⟩，态矢量写作 ψ(t) = c₀(t)|0⟩ + c₁(t)|1⟩。在旋转波近似下，驱动场频率与能级差共振时，含时薛定谔方程化为：

iℏ d/dt [c₀]   =   ℏ/2 [ 0    Ω ] [c₀]
        [c₁]           [ Ω   0  ] [c₁]

其中 Ω 为拉比频率。取 ℏ = 1 的自然单位，方程简化为 i dc/dt = H c，H = (Ω/2)·σₓ。解析解为 c₀(t) = cos(Ωt/2)、c₁(t) = −i sin(Ωt/2)，因此激发态概率 P₁(t) = sin²(Ωt/2)，这就是拉比振荡：概率在 0 与 1 之间以角频率 Ω 周期往返。这个解析解正是后面检验数值格式的“标准答案”。

数值求解的关键一步是把时间导数离散化。设时间步长为 Δt，记第 n 步的态矢量为 cⁿ。最直接的做法是用前向差分近似导数：

i (cⁿ⁺¹ − cⁿ) / Δt = H cⁿ
⇒ cⁿ⁺¹ = cⁿ − i Δt H cⁿ = (I − i Δt H) cⁿ

这就是显式欧拉法。它每一步只需一次矩阵-向量乘法，实现极简，但它是条件稳定的：对纯虚特征值的薛定谔方程，|1 − iΔtE| = √(1 + (ΔtE)²) > 1，模长会随步数单调增长，概率不守恒。步长稍大，激发态概率就会“溢出”到大于 1，振荡幅度被虚假放大。

二、Crank-Nicolson：用隐式平均换稳定性

为克服欧拉法的模长漂移，Crank-Nicolson 格式对时间导数在 n 与 n+1 之间取平均：

i (cⁿ⁺¹ − cⁿ) / Δt = H (cⁿ⁺¹ + cⁿ) / 2
⇒ (I + i Δt H / 2) cⁿ⁺¹ = (I − i Δt H / 2) cⁿ

该格式是酉的：演化算符 U = (I + iΔtH/2)⁻¹(I − iΔtH/2) 满足 U†U = I，因此每一步都严格保持概率守恒，且对任意步长无条件稳定。代价是每步需要解一个线性方程组。对二能级这种 2×2 小矩阵，可以直接写出闭式更新，无需通用求解器：

import numpy as np

Omega = 1.0          # 拉比频率
dt    = 0.05         # 时间步长
T     = 20.0         # 总时长
N     = int(T / dt)

H = 0.5 * Omega * np.array([[0, 1], [1, 0]], dtype=complex)
I = np.eye(2, dtype=complex)

# 欧拉法
c_eu = np.array([1.0, 0.0], dtype=complex)
P1_eu = []
for _ in range(N):
    c_eu = (I - 1j * dt * H) @ c_eu
    P1_eu.append(abs(c_eu[1])**2)

# Crank-Nicolson
A = I + 1j * dt * H / 2
B = I - 1j * dt * H / 2
c_cn = np.array([1.0, 0.0], dtype=complex)
P1_cn = []
for _ in range(N):
    c_cn = np.linalg.solve(A, B @ c_cn)
    P1_cn.append(abs(c_cn[1])**2)

这段代码把两种格式放在同一组参数下对照。注意初态取 c₀ = 1、c₁ = 0，即系统从基态出发；Ω = 1、dt = 0.05 是教学示例参数，实际仿真中应根据振荡周期 T_rabi = 2π/Ω 选取步长，通常要求每周期至少几十到上百步。

三、结果判读与格式对比

运行后应重点观察三件事：第一，P₁(t) 是否呈现周期约 2π/Ω 的正弦平方振荡；第二，|c₀|² + |c₁|² 是否恒为 1；第三，数值曲线与解析解 sin²(Ωt/2) 的偏差是否随步长减小而收敛。欧拉法在 dt = 0.05 时，若干周期后概率和会明显偏离 1，振荡峰值超过 1；Crank-Nicolson 则在整个时间窗内保持概率和接近机器精度，峰值稳定在 1 附近。若把 dt 减半，欧拉法的漂移会变慢，但不会消失——这正是“条件稳定”与“无条件酉”的本质区别。

上图把两种格式的更新公式与 P₁(t) 曲线放在同一视野中：绿色实线是解析解，红色虚线是欧拉法在较大步长下的典型失真——峰值越冲越高、谷值越陷越深，概率和不再为 1。Crank-Nicolson 的曲线则与解析解几乎重合。需要说明的是，图中曲线为教学示意，用于表现格式差异的趋势，具体数值以实际运行代码为准。

四、实现要点与常见陷阱

第一，矩阵与态矢量必须用复数类型。若用 float 数组，虚部会被静默截断，振荡直接消失。第二，H 的构造要保证厄米性，二能级情形下 σₓ 的系数取 Ω/2 才能让解析解与 P₁ = sin²(Ωt/2) 对应；若误写成 Ω，振荡周期会差一倍。第三，Crank-Nicolson 中 A = I + iΔtH/2 与 B = I − iΔtH/2 的符号不能对调，否则得到的是反向时间演化。第四，判断收敛时不要只看某一时刻的 P₁，而应扫描多个 Δt，观察最大概率偏差随 Δt 的下降阶数：欧拉法一阶收敛，Crank-Nicolson 二阶收敛。

把二能级拉比振荡跑通，就掌握了含时薛定谔方程数值求解的最小闭环：写出哈密顿量、选择时间离散格式、迭代更新态矢量、用守恒量与解析解双重校验。后续更复杂的多能级、含时驱动或空间格点问题，都可以沿用这套“离散—迭代—校验”的骨架，只是把 2×2 矩阵换成更大规模、把闭式更新换成稀疏线性求解。