1. 课程导入与仿真环境搭建

量子动力学仿真研究的是量子态随时间的演化行为。与量子力学课程中大量推导解析解不同，仿真路线关注的是：给定一个哈密顿量 H、一个初始态 |ψ(0)⟩，如何用数值方法求出 |ψ(t)⟩，并从结果中提取可观测量的时间序列。因此，本课程的核心对象可以归纳为三件东西：态矢量（或密度矩阵）、哈密顿量（决定演化的算符）、时间推进算法（把连续演化离散化的数值格式）。学习路径也围绕这三者展开：先建立数值线性代数的基本操作，再掌握含时薛定谔方程的积分方法，最后进入开放量子系统、周期驱动、量子淬火等更复杂的场景。

在动手写第一行仿真代码之前，需要先明确一个贯穿全课程的单位约定：本课程统一取约化普朗克常数 ħ = 1，能量与角频率同量纲。这样含时薛定谔方程写作 i d|ψ⟩/dt = H|ψ⟩，数值实现时不必反复乘除常数。时间单位与能量单位互为倒数，例如哈密顿量矩阵元取 1.0 时，对应的特征角频率就是 1.0，演化周期约为 2π。这个约定在后续所有章节和图示中保持一致。

科学计算栈的安装与验证

Python 生态中承担量子动力学仿真的主力是三个库：NumPy 提供复数数组与矩阵运算，SciPy 提供常微分方程求解器和稀疏矩阵工具，Matplotlib 负责把演化结果画成可读的曲线。推荐使用虚拟环境隔离依赖，避免与系统 Python 冲突：

python -m venv qdyn-env
source qdyn-env/bin/activate        # Windows 下为 qdyn-env\Scripts\activate
pip install numpy scipy matplotlib jupyter

安装完成后不要急着写演化脚本，先做一次最小验证，确认版本与关键功能可用：

import numpy as np, scipy, matplotlib
print(np.__version__, scipy.__version__, matplotlib.__version__)
# 验证复数矩阵指数与特征分解
sigmax = np.array([[0, 1], [1, 0]], dtype=complex)
vals, vecs = np.linalg.eigh(sigmax)
print(vals)          # 期望 [-1.  1.]

若 eigh 返回的特征值顺序为 [-1, 1]，说明线性代数后端工作正常。这里特意用 eigh 而不是 eig，因为量子力学中的哈密顿量是厄米矩阵，eigh 利用厄米性给出实数特征值和正交归一特征向量，数值上更稳定，这也是后续章节反复使用的习惯。

Jupyter Notebook 的配置要点

Jupyter 适合本课程的原因是它能把公式、代码、图形和文字解释放在同一页面，便于边改参数边看结果。启动方式为 jupyter notebook 或 jupyter lab。建议在第一个单元格中固定随机种子并设置绘图风格，保证结果可复现：

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(2024)
plt.rcParams['figure.dpi'] = 110
plt.rcParams['axes.grid'] = True

一个常见问题是 Notebook 中修改了外部模块却看不到更新，此时可用 %load_ext autoreload 与 %autoreload 2 自动重载。另一个问题是中文字体缺失导致图例出现方块，可在绘图前设置 plt.rcParams['font.sans-serif'] = ['SimHei']（或系统已有中文字体），并设置 axes.unicode_minus = False 以正常显示负号。

第一个量子态演化脚本

下面用一个二能级系统完成第一次完整演化。取哈密顿量 H = 0.5·ω·σx，其中 ω = 2.0，初始态为 |0⟩ = [1, 0]ᵀ。由于 σx 的本征态是 |±⟩，初始态在两个本征态上等幅叠加，演化表现为布居数在 |0⟩ 与 |1⟩ 之间以角频率 ω 振荡，即拉比振荡。解析结果为 P₀(t) = cos²(ωt/2)，可作为数值结果的对照基准。

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

omega = 2.0
H = 0.5 * omega * np.array([[0, 1], [1, 0]], dtype=complex)
psi0 = np.array([1.0, 0.0], dtype=complex)

def rhs(t, y):
    psi = y[0::2] + 1j * y[1::2]      # 复数拆成实部/虚部两路
    dpsi = -1j * (H @ psi)
    return np.concatenate([dpsi.real, dpsi.imag])

t_eval = np.linspace(0, 2*np.pi, 400)
sol = solve_ivp(rhs, (0, 2*np.pi), np.concatenate([psi0.real, psi0.imag]),
                t_eval=t_eval, rtol=1e-9, atol=1e-12)

psi_t = sol.y[0::2, :] + 1j * sol.y[1::2, :]
P0 = np.abs(psi_t[0, :])**2
P1 = np.abs(psi_t[1, :])**2
print("最大归一化偏差:", np.max(np.abs(P0 + P1 - 1)))

这里有一个关键实现细节：solve_ivp 的求解器默认处理实数向量，而量子态是复数，因此把复数拆成实部与虚部拼接成实向量送入求解器，回调时再重组。这是量子动力学数值实现中最常见的工程处理方式之一。运行后应看到最大归一化偏差在 1e-10 量级，说明数值积分保持了概率守恒；若偏差明显偏大，通常意味着容差设置过松或哈密顿量不是厄米的。

把 P0 与解析式 cos²(ωt/2) 画在同一张图上，两条曲线应当几乎重合。这个对照步骤不是形式主义：它同时验证了哈密顿量构造、复数拆分、积分器调用和结果提取四个环节，是后续更复杂模型的可信度基线。

从首次运行到后续章节的衔接

这次运行虽然只有二能级、单一哈密顿量，但它已经包含了后续所有仿真的骨架：构造厄米矩阵、准备归一化初态、选择积分器与容差、提取可观测量、与解析或已知极限对照。后续章节会在这个骨架上逐步替换部件——把 σx 换成更复杂的多能级或格点哈密顿量，把闭系演化换成含耗散的林德布拉德方程，把固定参数换成周期驱动或参数扫描。因此，第一次运行是否“跑通”并不是终点，真正要建立的是可复现、可对照、可扩展的工作习惯：参数集中写在单元格顶部、随机种子固定、每次结果都与一个可解释的基准比较。