企业官网建设流程全解析

1. 拆解标题：BFM、Wasserstein 梯度流与“Windows 安装”这三件事的关系

一个能跑出漂亮收敛曲线的数值方法，如果无法在本地工作站运行，就缺少实际验证机会。BFM（The Back-And-Forth Method，前后向法）用于计算 Wasserstein 梯度流，本质是把一个难以直接求解的隐式近端问题拆成前后交替的两半，反复迭代到不动点。这个思路很干净，但参考实现几乎都默认类 Unix 环境：gcc、make、shell 脚本、fork 多进程、正斜杠路径。搬到 Windows 上，方法本身一行没变，环境却要重做一遍。这篇内容讲的是把 BFM 的 Wasserstein 梯度流代码在本机 Windows 上完整跑起来：它解决什么问题、依赖怎么安装、编译链怎么配、参数怎么取、Windows 上有哪些典型问题。无论你刚接触最优传输，还是已经会用数值库写梯度流，都能找到可直接照做的部分。

1.1 用一句话讲清 Wasserstein 梯度流在算什么

把一个概率分布想象成地面上的一堆沙子，用一个“自由能”泛函衡量这堆沙子的状态好坏，比如势能高不高、粒子挤不挤。当沙子自发往低能量状态滚动时，由谁决定沙粒怎么移动？答案是最省力的搬运方式：每一步都选择单位距离成本最小的传输方案，这套成本度量就是 Wasserstein 距离。把概率测度空间配上 W2 距离后，它具备类似弯曲空间的几何结构，泛函沿 W2 测地线方向的最速下降，就是 Wasserstein 梯度流。写成方程是这样的：状态密度 ρ 随时间演化，满足连续性方程 ∂ₜρ = ∇·(ρ ∇(δF/δρ))，其中 F 是自由能，δF/δρ 是它的变分导数。用生活化方式解释：每时每刻，每个粒子都朝“降低整体能量”的方向移动，但移动必须保持总质量守恒，所以只能形成流的形状，不能在原地凭空产生或消失。

实际计算时，连续时间无法直接处理，标准做法是做时间离散，得到所谓最小化运动格式：给定上一步的分布 ρₖ，下一步的解是 ρₖ₊₁ = argmin_σ [ F(σ) + W2²(σ, ρₖ) / (2τ) ]。这个式子的含义很直白：既要让能量 F 降下去，又不能离上一步太远，否则要付一笔与 W2 距离平方成正比的运费，τ 是时间步长，相当于运费折扣。τ 越小越接近连续流，但迭代次数越多；τ 越大跑得越快，但隐式方程的近似误差会明显。问题在于，这个 argmin 是在整个概率测度空间上取的，还带一个非光滑的 W2 平方项，直接求解几乎不可能。这就引出了 BFM 的动机。

1.2 BFM 到底省掉了哪一步

普通的近端迭代（也就是把上面那个 argmin 用各种方式近似）通常要正面处理“能量下降 + 传输约束”这两个耦合在一起的目标，耦合带来的麻烦是：你让能量下降的方向未必是 W2 意义下最近的方向，两者会互相冲突，导致迭代振荡甚至发散。BFM 的处理思路是把这一步拆开：先按能量下降的方向做一个前向（forward）提议，得到中间状态；再用后向（backward）步骤把这个提议按输运意义拉回去，检验它是否真的让整体目标（能量 + 运费）变小了，如果变小了就接受，没变小就缩短步子重试。两个方向的子问题各自都是显式的、可并行的，交替几次到残差足够小，就得到了这一步的近端点。

这个拆法并非数学技巧，工程上有三个实实在在的好处。第一，前向步可以做得很便宜，就是一次梯度下降，甚至可以在 GPU 上一次性算完所有粒子；第二，后向步本质上是一次回溯线搜索加收缩，天然保证了能量单调不上升，避免显式格式在步长偏大时数值爆炸；第三，前后向的步长比例本身就是一个诊断量——如果后向步平均只肯接受一半以下的步长，说明你的 τ 取大了，或者交互核的尺度选错了。我在 Windows 上第一次跑通之后，最先盯的就是这个接受率曲线，它比看损失曲线直观得多。要注意的是，不同文献里“前向”“后向”的具体定义并不完全统一，有的把它写成两个交替优化的映射，有的写成两个半隐式更新；我下面给的实现是工程上最常见的那种松弛版本，思路一致，但并非对某篇论文的逐行复现。

1.3 这套东西的适用面与影响范围

从应用角度看，Wasserstein 梯度流这一类方法真正有价值的地方，是把“分布层面的优化”变成了可以被数值方法处理的动力系统。传统做法里，做贝叶斯后验采样用的是马尔可夫链，收敛慢且难以并行；用梯度流视角的粒子方法，一次迭代能同时更新上万个粒子，而且天然支持多峰目标——因为在输运意义下，粒子是被整体推着走的，不会像某些采样器那样陷在一个峰里出不来。这也是它在变分推断、生成模型训练、平均场神经网络、点云与图像配准、人群与粒子系统仿真这些方向被反复使用的原因。

BFM 在其中的定位是“把隐式格式做得能动”。纯显式格式简单，但步长被稳定性条件卡死；纯隐式格式理论好，但每步都要解一个不可解的非线性问题。前后向交替恰好卡在中间：既保留隐式格式的单调性，又只需要显式格式的计算量。影响范围上，它适合做中大规模粒子的长时间演化，比如粒子数几千到几万、迭代几千步的规模；如果只是几十个粒子、几十步的演示，用什么格式差别不大。我在实际使用中发现，规模一上去，前后向交替省下来的时间很可观，因为后向步的回溯在大多数步只需试一两次就能通过，成本几乎可以忽略。

1.4 为什么“Windows 安装”值得单独写

严格来说，装个 Python 环境跑个脚本不算难事。麻烦在于这类数值方法的参考代码通常带一堆隐式假设：它假设你能make，假设路径分隔符是正斜杠，假设多进程fork可用，假设某个 C 扩展在安装时能自动找到编译器。Windows 上这些假设一条都不成立。更麻烦的是报错信息往往指向错误的地方——比如提示cl.exe找不到，真实原因可能是环境变量里 MSVC 的路径被另一个工具链覆盖了；提示某个dll加载失败，真实原因可能是两套 OpenMP 运行时同时被拉进来打架。

再加上前置的数学库（比如最优传输的计算库）在 Windows 上并不总是有现成的二进制包，逼着你去配编译链；而一旦决定配编译链，Visual Studio Build Tools 这个大块头的勾选项又多到让人怀疑人生——装多了占几十 GB，装少了某个头文件找不到。这些坑与 BFM 本身无关，但每一个都足以让人在跑通之前放弃。所以这篇内容的重点会放在环境与方法之间的衔接上：先让你知道代码在算什么，再让你知道为什么 Windows 上要这么装。

2. Windows 侧的准备工作：把不可能变成可编译

2.1 机器自检：先算一遍内存和算力账

很多人第一反应是看显卡，其实这类粒子方法最先卡住的是内存，因为最常见的实现方式要算粒子两两之间的交互核，规模是 N² 的。假设粒子数 N = 4096，用双精度存一个 N×N 的矩阵是 4096² × 8 字节 = 134 MB，还只是单个数组；如果你的代码里同时存在x[:, None] - x[None, :]的差值矩阵、核矩阵、梯度矩阵三个副本，瞬时占用就接近 400 MB。N 提到 16384，单数组就是 2.1 GB，三个副本直接超过 6 GB，加上 Python 本身和其他缓存，16 GB 内存的机器很容易触发内存交换，表现就是越跑越慢、最后卡死。

所以动手之前先做一道除法：用你打算用的粒子数平方乘以 8，再乘以 3 到 4，估算峰值内存，对比一下自己机器的可用内存。如果不够，有三条路：一是降低粒子数、增加迭代步数（经验上粒子数减少带来的误差是 O(1/√N)，比时间步长的 O(τ) 好控制一些）；二是把交互核换成分箱或直方图近似，把 O(N²) 降到 O(N)；三是把交互项矩阵分块计算，用循环分块把峰值压到 O(N × B)。GPU 这边，显存同理，而且 Windows 上显存的共享机制更容易让你误判——任务管理器显示的“GPU 内存”包含共享内存，实际能用的独显容量要看专用显存那一栏。

CPU 核数决定了你能并行多少条粒子轨迹。纯 NumPy 实现基本吃满单核，想做并行要么上多线程 BLAS，要么自己拆任务。这一点在 Windows 上要特别注意，因为多进程的启动方式和 Linux 完全不同，后文会专门讲。至于 CUDA，只有在你确认要跑 GPU 后端时才需要，装之前先确认驱动版本与新版本 CUDA 的对应关系，别装完了发现驱动太旧不兼容。

2.2 Python 环境选型：为什么我最终选了 conda 而不是 venv

venv 轻、干净，是纯 Python 项目的首选。但 BFM 这类数值代码会踩到一堆非纯 Python 依赖：BLAS 后端、OpenMP 运行时、稀疏矩阵库、某些需要编译的 C/C++ 扩展。venv 只管 Python 包，遇到这些二进制依赖时要么自己编，要么等别人发 wheel，在 Windows 上等 wheel 的过程很折磨。

conda（我一般用精简发行版）的好处是它把二进制依赖也当成包来管，BLAS、OpenMP、FFT 库都有预编译版本，装上就能用，路径也不需要你手动配置。代价是环境稍微臃肿一些，包体积偏大。我的取舍标准很简单：如果你只是跑一个纯 Python 的小脚本，venv 足够；如果你要装科学计算栈、要跨多个数值库、还要在这个环境里编译扩展，那 conda 省下来的时间远大于它占的磁盘。Python 版本建议钉在 3.10 或 3.11。这两个版本的轮子覆盖最全，3.12 之后虽然新，但某些偏门的科学计算包还没有对应的二进制包，你得自己编，而这恰恰是把人卡住的地方。

还有一个细节：装环境时别把 base 环境搞乱。所有东西都装进命名环境，不要习惯性pip install到 base，否则几个月后你会面对一个完全无法复现的环境。我给自己定了个规矩，只要涉及 BFM 相关实验，一律conda activate到专用环境，退出时直接关掉终端，不留在 base 里瞎操作。

2.3 编译工具链的最小安装清单

只有当某个依赖必须以源码方式安装时，你才需要 MSVC。但“什么时候需要”往往事先不知道，所以我的建议是先装好一个最小可用集合，别等报错了再补。去下载安装 Visual Studio Build Tools（社区版即可），在组件勾选页面上只装三样：

• “使用 C++ 的桌面开发”工作负载下的MSVC v143 - VS 2022 C++ x64/x86 生成工具（如果你要编译 32 位就再加 x86，一般不需要）
• Windows 10 SDK或Windows 11 SDK（选一个与你系统匹配的即可，两个都装会多占几个 GB）
• C++ CMake 工具（某些包用 CMake 作为构建后端，装上省事）

其他所有可选项，包括 .NET、Web 开发、移动开发、Unity 相关，一律不勾。这样装完大约几个 GB，而不是几十 GB。装完之后不要指望直接在普通命令行里cl就能用，Windows 上 MSVC 的环境变量必须通过专门的“x64 Native Tools Command Prompt for VS 2022”来初始化，或者你在普通 PowerShell 里手动调用它的环境初始化脚本。我习惯用前者，因为路径是自动配好的，不容易漏。

验证方式很直接：在那个命令行窗口里敲cl，如果输出一行编译器的版本信息，说明工具链就位；如果提示不是内部或外部命令，那就是路径没生效，回到安装器里确认有没有装“生成工具”这一项。另外建议顺手把长路径支持打开，因为这类项目在 conda 环境里嵌套几层之后，路径长度轻松突破 260 字符，Windows 默认的限制会让解压和编译在莫名其妙的地方失败。这个开关在组策略里，用管理员权限改完之后要重启才生效。

2.4 三条落地路线横向对比

到这一步你其实有三个选择，各自的成本和风险差别很大，我把它们放在一张表里方便对比：

我的建议是分两步走：先用第三条路线把方法跑通、确认结果对，再决定要不要投入时间走第二条路线。第一和第二条路线并不冲突，你完全可以在原生 Windows 上做日常调参和画图，把大规模实验放到 WSL2 里跑，两边共享同一份代码目录（注意不要在子系统里对 Windows 路径做高强度 I/O，会很慢）。有一点需要提前说明：某些主打即时编译和自动微分的数值库在 Windows 原生环境下没有官方发布的轮子，装的时候会直接卡住，这种情况下要么切到 WSL2，要么换用别的后端，别硬刚。

3. 依赖安装与代码落地实操

3.1 建环境与版本钉死

我一般用下面这几条命令起环境。之所以把 Python 版本写死在创建命令里，是因为后续所有包的版本选择都依赖它，一旦版本漂移，之前踩过的坑要重新踩一遍：

conda create -n bfm python=3.11 -y conda activate bfm python -c "import sys; print(sys.version)"

装包之前先把下载源指向一个在你所在网络环境下响应最快的节点，这一步能省掉大量“下载到一半卡住”的时间。具体地址以你自己的实测为准，不必照搬别人的。设置完可以用pip config list确认生效。

接着按下面的顺序装，顺序是有讲究的：先装数值基础库，再装科学计算扩展，最后装可视化和业务相关的库。原因是某些库在安装时会去探测已装的 NumPy 版本并据此决定编译参数，顺序反了容易装出一个 ABI 不匹配的环境：

conda install -c conda-forge numpy scipy matplotlib scikit-learn tqdm pyyaml -y

这五个是几乎所有实现都绕不开的：NumPy 提供数组和向量化，SciPy 提供距离计算、插值和稀疏矩阵，matplotlib 负责出图，scikit-learn 提供一些现成的密度估计和聚类工具，tqdm 负责进度条（迭代几千步的实验没有进度条会让人焦虑），pyyaml 用来读配置文件。把它们装在 conda 层而不是 pip 层，是为了拿到带 MKL 或 OpenBLAS 加速的二进制版本，矩阵运算能快出一截。

3.2 科学计算与最优传输依赖的安装顺序

真正可能卡住的是最优传输计算库。它内部有 C++ 实现，Windows 上的轮子覆盖不如 Linux 全，直接pip install有时会走到源码编译，然后因为找不到编译器或缺少某个头文件而失败。我的做法是优先走 conda 通道，实在没有再回退：

conda install -c conda-forge pot -y

装完之后立刻验证，而不是等到跑实验时才验证：

import ot import numpy as np a = np.ones(100) / 100 b = np.ones(100) / 100 x = np.linspace(0, 1, 100).reshape(-1, 1) y = np.linspace(0, 1, 100).reshape(-1, 1) M = ot.dist(x, y) print(ot.emd2(a, b, M))

这一小段在做的是一维均匀分布到自身的传输，理论最优值应当是 0。如果输出接近 0，说明库和底层编译都没问题；如果直接抛异常，多半是 C++ 运行时缺失或者 DLL 加载失败，这时候装一个通用 C++ 运行库通常能解决。这一步别跳过，我见过太多人在写了三百行代码之后才发现是库没装好，然后去改算法，方向完全错了。

GPU 后端是可选的。如果你的实现提供 CuPy 或 PyTorch 后端，直接pip install cupy-cuda12x（其中 12x 对应你装的 CUDA 大版本）或从官方渠道装 PyTorch 的 CUDA 版本。这里有个 Windows 特有的注意点：不要在同一个环境里同时装 CuPy 和 PyTorch 的 CUDA 版本，两者可能链接到不同版本的 CUDA 运行时，表现是运行到一半随机崩溃且没有任何报错信息。真要同时用，就分两个环境。

3.3 源码目录结构与文件职责

把代码拿到本地之后，先别急着运行主脚本，花五分钟看一眼目录结构。这类项目通常长这样：

bfm-wasserstein/ ├── configs/ # 参数配置文件，一维和高维实验各一份 ├── data/ # 初始分布、预计算的目标测度 ├── src/ │ ├── energy.py # 自由能泛函与变分导数 │ ├── operators.py # 前向步、后向步、交替调度 │ ├── solver.py # 外层时间循环与收敛判据 │ └── utils.py # 可视化、统计量、日志 ├── scripts/ # 命令行入口 └── tests/ # 单元测试，先跑这个

要重点看的是energy.py和operators.py：前者决定了你的自由能里有几项（势能、交互能、熵项），后者决定了前后向是怎么交替的。多数实现里的能量函数是带参数化系数的，比如F = ∫V dρ + β·交互能 + ε·熵，这三个系数直接决定最终分布的形状，比时间步长还重要。配置文件里如果有beta、sigma、eps这类名字，先弄清楚它们对应哪一项，别拿别人的默认值直接跑自己的问题。

还有一个容易忽略的地方：tests/目录。很多项目会在这里放一两个极小规模的正确性检查，比如“一维高斯到高斯，跑 100 步之后均值应当接近理论值”。先跑测试，快则几十秒，慢则几分钟，但它能帮你把环境问题和算法问题彻底分开。我在 Windows 上第一次跑的时候，测试挂了三次，前两次都是路径相关的，第三次才发现是核函数的尺度参数单位不一致——如果我直接跑主实验，会以为是算法不收敛，然后去调一堆无关的参数。

3.4 后端切换与配置文件写法

配置文件用 YAML 是最省心的，因为它的键值对读起来清楚，改起来也不会破坏缩进结构（只要你用空格不用 Tab）。下面是我常用的一份一维实验配置，注释里写清了每一项的作用：

seed: 20240117 particles: 1024 # 粒子数，先小规模验证再放大 dim: 1 # 空间维度 tau: 0.002 # 时间步长，见参数推导 steps: 400 # 外层时间迭代次数 inner: 4 # 前后向交替的最大轮数 tol: 1.0e-6 # 不动点残差阈值 energy: beta: 0.05 # 交互项系数 sigma: 0.25 # 交互核的空间尺度 eps: 0.0 # 熵项系数，一维演示先置零 backend: numpy # numpy 或 gpu，GPU 需额外装后端 output: results/run01

这里有一个在实践中被反复验证的习惯：backend单独做成一个开关，而不是写死在代码里。因为同一个实验在 CPU 和 GPU 上的数值结果应当一致（浮点误差范围内），这个开关正好用来做交叉校验。如果两个后端结果差得离谱，那一定是某一边的实现有 bug，而不是“GPU 就是不准”这种想当然的解释。

还有output这一项，务必让每次运行写到独立目录。这类实验动辄几百上千步，配上绘图和日志，一次输出就是几百 MB，混在一起之后你根本分不清哪个是哪个。我一般用“日期 + 关键词 + 序号”的命名，跑之前先看目录里有没有重名，有就顺延。

4. 跑通第一个可复现实验

4.1 环境自检脚本：一个一维双阱的 BF 迭代

下面这段代码是我在每台新机器上都会跑一遍的“环境自检 + 方法验证”脚本。它不依赖任何编译型扩展，只用 NumPy 和 matplotlib，所以只要环境能 import 这两个库，它就能跑。它做的事是：在一维双阱势能加上高斯排斥核的自由能下，用前向提议 + 后向回溯的方式迭代粒子位置，最后画出一张分布演化的对比图。核心结构与 BFM 一致——前向步按能量下降方向走，后向步检查整体目标有没有下降并据此收缩步长，交替到不动点。

import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams["font.sans-serif"] = ["Microsoft YaHei"] plt.rcParams["axes.unicode_minus"] = False def potential(x): """双阱势 V(x) = (x^2 - 1)^2，最小值在 x = ±1""" return (x ** 2 - 1.0) ** 2 def grad_potential(x): """V 对 x 的梯度""" return 4.0 * x * (x ** 2 - 1.0) def interaction_force(x, sigma): """核 phi(r)=exp(-r^2/2sigma^2) 产生的排斥力，均值场近似""" d = x[:, None] - x[None, :] w = np.exp(-(d ** 2) / (2.0 * sigma ** 2)) / (sigma ** 2) return (d * w).sum(axis=1) / x.size def free_energy(x, beta, sigma): """自由能：势能项 + 交互项""" v = potential(x).mean() d = x[:, None] - x[None, :] k = np.exp(-(d ** 2) / (2.0 * sigma ** 2)) return v + 0.5 * beta * k.sum() / (x.size ** 2) def forward_op(x, tau, beta, sigma): """前向：沿能量下降方向走一步""" return x - tau * (grad_potential(x) - beta * interaction_force(x, sigma)) def backward_op(x_old, proposal, beta, sigma, shrink=0.5, max_try=20): """后向：回溯收缩，保证自由能不上升""" f_old = free_energy(x_old, beta, sigma) alpha = 1.0 for _ in range(max_try): cand = x_old + alpha * (proposal - x_old) if free_energy(cand, beta, sigma) <= f_old: return cand, alpha alpha *= shrink return x_old.copy(), 0.0 def run(x0, tau, beta, sigma, steps, inner=4, tol=1e-6): x = x0.copy() log = [] for k in range(steps): alpha = 0.0 for _ in range(inner): proposal = forward_op(x, tau, beta, sigma) x_new, alpha = backward_op(x, proposal, beta, sigma) res = np.linalg.norm(x_new - x) / max(np.linalg.norm(x), 1e-12) x = x_new if res < tol: break log.append((k, free_energy(x, beta, sigma), alpha)) return x, np.array(log) if __name__ == "__main__": rng = np.random.default_rng(20240117) n, tau, beta, sigma, steps = 1024, 2e-3, 0.05, 0.25, 400 x0 = rng.normal(loc=3.0, scale=0.4, size=n) x, log = run(x0, tau, beta, sigma, steps) fig, axes = plt.subplots(1, 2, figsize=(10, 3.6)) axes[0].hist(x0, bins=60, density=True, alpha=0.55, label="初始分布") axes[0].hist(x, bins=60, density=True, alpha=0.55, label="演化后分布") axes[0].legend() axes[0].set_title("粒子分布对比") axes[1].plot(log[:, 0], log[:, 1], lw=1.2) axes[1].set_title("自由能随迭代变化") axes[1].set_xlabel("迭代步") axes[1].set_ylabel("F") fig.tight_layout() fig.savefig("bfm_check.png", dpi=150) print("平均接受步长比例:", log[50:, 2].mean().round(3))

这段脚本刻意写得直白，没有做任何向量化优化，因为它的用途是验证环境和确认方法逻辑，不是跑大规模实验。实际跑下来，1024 个粒子、400 步、内层最多 4 轮，在普通笔记本上十几秒就能出图，非常适合作为“这台机器能不能干活”的判断依据。

4.2 参数取值的计算过程

参数不是拍脑袋定的，每一个都有推导或经验依据，我把关键几项的取法写清楚，方便你换问题的时候按同样的逻辑重算。

时间步长的上限由势能项的曲率决定。显式格式的局部误差放大因子近似为 |1 - τ·V''(x)|，要求它落在 (-1, 1) 之间，也就是 τ < 2 / max|V''|。上面的双阱势 V''= 12x² - 4，粒子初始在 x = 3 附近，V'' ≈ 104，所以 τ < 2/104 ≈ 0.019。我取 2e-3，只有上限的十分之一，留足了安全余量。这不是保守，而是因为交互项的曲率不好估，一旦加上排斥核，实际的有效曲率会更大，步长太激进会在中后期突然发散。

交互核的尺度 σ 要和粒子间距匹配。1024 个粒子分布在一个宽度大约 4 到 5 的区间上，平均间距约 0.005，σ 取 0.25 意味着每个粒子能“看到”周围半径 0.25 范围内的约 50 个粒子，核矩阵不会退化成一个近乎对角阵（那样排斥力几乎为零），也不会宽到把整个分布当成一个点（那样变成全局平移）。经验上取平均间距的 20 到 50 倍比较合适。

交互项系数 β 决定排斥强度和势能约束的平衡。β 太大，粒子互相推得太远，分布摊得很平；β 太小，粒子全部塌到 ±1 两个点上，退化成两个 δ 峰。0.05 这个值是我在小规模实验里扫出来的：在 0.01 到 0.2 之间取若干个值各跑一次，看最终分布是“两个有宽度的峰”而不是“两条竖线”或者“一个大平台”，就停下。

粒子数 N 与统计误差的关系是 O(1/√N)。N = 1024 时相对误差约 3%，用来做方法验证足够；要出论文级别的图，至少 2¹⁴ 以上，此时内存开销按 2.1 节算出来的账要重新核对。内层交替轮数inner不用取太大，一般 3 到 5 轮就够，因为前后向交替的收敛速度在前几轮最快，后面是线性收敛的尾巴，收益递减。

4.3 结果怎么判断是对的

跑出图不等于做对了。至少要过三道检查。第一道是能量单调性：自由能曲线必须单调不增，如果中间出现上升，说明后向回溯逻辑有 bug（或者free_energy和forward_op用的不是同一个泛函，这是最常见的错误，尤其是梯度推导里漏了个系数）。第二道是终态形状：双阱势加排斥核的平衡态应当是两个峰，峰位在 ±1 附近，峰之间有一定宽度。如果只出现一个峰，说明粒子全跑到了势能的同一侧，通常是初始分布太集中或者时间步太长，局部能量地形被跨过去了。

第三道是和解析解或参考实现对比。一维情况下有些特例可以精确算，比如自由能只含势能项且势能是二次型时，解是高斯，均值和方差都有闭式表达式，可以直接算差值。哪怕不做严格对比，至少做一次收敛性检查：把 τ 减半、粒子数翻倍，看结果变化是否在预期量级内（时间误差 O(τ)，统计误差 O(1/√N)）。如果减半之后结果剧烈变化，说明之前的结果只是数值假象。

我在 Windows 上第一次跑通时，能量曲线在 200 步左右有个极小的回升，只有