第一次在 Windows 上跑 Wasserstein 梯度流的数值实现，最容易让人放弃的往往不是那套测度论，而是pip install之后满屏的红色报错。我在三台不同配置的机器上装过这套东西，最夸张的一次卡在一个DLL load failed上耗掉整个下午，最后发现跟算法毫无关系——只是解释器版本和预编译轮子对不上。这篇就把 BFM（The Back-And-Forth Method）在 Windows 上的落地过程完整走一遍：它到底在“来回”什么、环境该怎么搭才不返工、核心循环怎么写、以及跑起来之后密度变负、质量漂移、时间步爆炸这些坑怎么排。适合已经会写 Python、想动手复现 Wasserstein 梯度流数值格式的人，也适合被 Windows 编译环境折腾怕了、想找一条确定路径的读者。

1. 先把BFM讲成人话：它到底在"来回"什么

1.1 从能量泛函出发，梯度流想让什么东西下降

先把这个概念落到地面上。给你一个密度 ρ 和一个能量泛函 E(ρ)，Wasserstein 梯度流要回答的问题是：如果我允许质量以"搬运"的方式移动，代价按搬多远、搬多少来算，那么让 E 下降最快的演化方向是什么？

答案是速度场 v = -∇(δE/δρ)，其中 δE/δρ 是能量对密度的第一变分。代进连续性方程，就得到

∂_t ρ + ∇·(ρ v) = 0, v = -∇(δE/δρ)

不同 E 对应完全不同的物理过程：E = ∫ρ log ρ 给的是热方程；再加上外势 ∫Vρ 就是 Fokker-Planck 方程；换成 E = (1/(m-1))∫ρ^m 就变成多孔介质方程；把吸引势和扩散项混在一起，就是聚集-扩散型方程。也就是说，同一个数值框架能覆盖一整类演化方程，这正是这套东西值得花时间的原因。

麻烦在于，这个方程高度非线性，而且经常退化——密度一旦在某处趋近于零，扩散系数就跟着趋近于零，显式格式立刻炸掉。直接对 PDE 硬上差分，基本是自找苦吃。

1.2 拆成两步：一次点态求值 + 一次椭圆方程求解

BFM 的核心思路是把这个非线性问题拆成两个线性子问题反复交替，这也正是"来回"这个名字的由来。

第一步是点态的，几乎不花计算量：给定当前的 ρ，直接算出 μ = δE/δρ，再取 w = -∇μ。这一步没有任何全局耦合，就是在网格上做一次差分。

但这里有个问题——w 不一定是个梯度场。它可能带着一个很大的旋度分量，而这个分量对密度的实际变化毫无贡献，只会在显式输运里白白消耗时间步、放大噪声。所以第二步要做一次加权亥姆霍兹投影：在 ρ 加权的最小二乘意义下，找最接近 w 的梯度场 ∇φ，也就是解

-∇·(ρ ∇φ) = -∇·(ρ w)

这是个标准的变系数椭圆方程，系数由 ρ 给出。解出来之后令 v = ∇φ，再拿这个 v 去推 ρ——推的那一步是线性的输运方程，保守格式写起来很成熟。

于是整个循环就是：算速度 → 投影成梯度场 → 推密度 → 回到算速度。两个子问题各自线性、各自有成熟的求解器，非线性和退化性全部被塞进了泊松方程的系数 ρ 里。这是它比直接硬解 PDE 稳得多的根本原因，不是玄学调参调出来的。

需要说明的是，不同文献里这套交替迭代的具体写法有出入，有的在拉格朗日视角下用映射复合，有的在欧拉视角下做通量更新。我这里给的这一版是最容易写、也最容易改的：全程欧拉网格，砍掉任何粒子重采样环节，Windows 上不会遇到额外的编译依赖。

1.3 为什么偏偏在Windows上格外容易翻车

这套代码的依赖横跨好几层：底层的 BLAS/LAPACK、稀疏矩阵求解器（SUPERLU/UMFPACK）、FFT 后端，再加上可选的加速扩展。每一层在 Windows 上的轮子覆盖率都不如 Linux 完整，而且失败的时候报出来的往往是一行 DLL 名字，跟算法一点关系没有。

更现实的问题是环境本身：Windows 没有系统级包管理器，Python 解释器来源五花八门，Microsoft Store 版和 python.org 版的路径规则还不一样；虚拟环境路径太深会撞上长路径限制；目录名带中文在某些老版本工具链上会直接崩。我见过最多的场景是：算法代码完全正确，但在别人机器上跑不起来，原因只是对方装了 3.13 而某个依赖还没有对应的轮子。

所以下面这一节虽然在讲"安装"，但它其实是整个项目能否推进的前置条件，值得认真对待。

2. Windows上的环境搭法：版本选择比装什么更重要

2.1 Python解释器：别追新，锁定3.11

这是我最想说的一条经验。科学计算生态的轮子发布永远慢半拍，新版本 Python 刚出来那几个月，稀疏求解、FFT 这些底层包大概率还没跟上，pip 会退化成源码编译，然后你就需要一整套 MSVC 工具链——而这本来是完全可以避免的。

我的建议是锁Python 3.11，3.10 也可以。3.12 大部分包已经齐了但个别加速库还在补，3.13 目前不建议碰。解释器来源优先选 python.org 的官方安装包，安装时务必勾上 "Add python.exe to PATH"，同时不要勾选 "Install for all users" 之外的花哨选项，也不要装 Microsoft Store 版本——那个版本的sys.executable指向一个带版本号的重定向路径，某些需要定位自身位置的包会在这上面出错。

装完先验证：

python -c "import sys, platform; print(sys.version); print(platform.architecture())"

输出里要看到64bit。32 位解释器在稀疏矩阵规模一上去就会因为内存寻址撞墙，而且报错信息很不直观。

2.2 venv 还是 conda：两条路的取舍我踩过

结论很直接：先走 venv + pip，除非你确实需要 MKL 那点性能。如果非要用 conda，记住一条铁律——在 conda 环境里绝对不要用 pip 去装 numpy/scipy，混装导致的符号冲突排查起来极其痛苦。

venv 的创建和激活：

cd D:\work\bfm python -m venv .venv Set-ExecutionPolicy -Scope CurrentUser RemoteSigned .\.venv\Scripts\Activate.ps1

第二行是 PowerShell 用户最容易漏的。默认执行策略是 Restricted，直接运行激活脚本会被拒绝，报错信息还带一堆安全提示，很容易让人误以为是环境坏了。设成RemoteSigned就够了，只影响当前用户，不需要管理员权限。

另外强烈建议把工作目录放在盘符根下一两层，比如D:\work\bfm，别放在C:\Users\某某某\Documents\我的项目\毕业设计\代码\第二版\...这种地方。路径一长，加上中文目录名，后面总有一处会出问题。

2.3 依赖安装顺序与镜像配置

装包的顺序是有讲究的，不是随便排的：

python -m pip install --upgrade pip python -m pip config set global.index-url https://<你的镜像地址>/simple python -m pip install --only-binary=:all: numpy python -m pip install --only-binary=:all: scipy python -m pip install --only-binary=:all: matplotlib

核心依赖其实就三个：numpy 负责数组和基础线性代数，scipy 提供稀疏矩阵和迭代求解器，matplotlib 用来画收敛曲线。numpy 必须排在 scipy 前面，因为 scipy 的轮子对 numpy 的 ABI 有版本约束，pip 在解析依赖时会按已装版本去挑对应的 scipy 轮子。

--only-binary=:all:这个参数很关键。它的作用是：如果没有可用的预编译轮子，就直接报错，而不是偷偷转去源码编译。听起来像是自找麻烦，但实际是帮你把问题提前暴露——你宁愿在装包阶段看到一句清晰的 "no matching distribution"，也不愿意在半小时编译之后看到一个 LNK2019。

这个项目不需要pyfftw。虽然它理论上有更快的 FFT，但 Windows 上要手动处理 FFTW 的 DLL，性价比极低。scipy 自带的 pocketfft 后端（scipy.fft）足够快，而且完全免配置。同理，numba 可以装（有现成轮子），但先别装，等确认性能确实是瓶颈再说。

验证一下装完的状态：

python -c "import numpy, scipy; print(numpy.__version__, scipy.__version__); import scipy.sparse.linalg as sla; print('ok')"

scipy.sparse.linalg单独 import 一次是有意义的，这个模块下挂着 SUPERLU 的绑定，如果轮子有问题，异常会在这一行暴露，而不是等你跑到一半才炸。

最后设两个环境变量，省掉后面一堆编码麻烦：

$env:PYTHONUTF8 = "1" chcp 65001

matplotlib 如果在脚本里画图，记得在 import 之前加一行matplotlib.use("Agg")，否则在某些终端里会因为没有交互后端直接卡住。

3. 把BFM主循环写出来：网格、投影、输运三件套

3.1 网格与初值：为什么用单元中心而不是格点

先定网格。取周期边界、二维方形区域 [0,1)²，N×N 个单元中心：

import numpy as np N = 128 L = 1.0 h = L / N x = (np.arange(N) + 0.5) * h X, Y = np.meshgrid(x, x, indexing="ij") def make_gaussian(s0): rho = np.exp(-(X - 0.5)**2 / (2 * s0) - (Y - 0.5)**2 / (2 * s0)) return rho / rho.sum() / (h * h) # 总质量归一化为 1

用单元中心而不是格点，原因是后面所有格式都是通量形式：密度存在单元里，速度定义在单元之间的面上，质量的进出正好抵消。用格点的话边界处理会多出一堆特例，很容易在质量守恒上出偏差。

indexing="ij"也要留意。默认是"xy"，会让第一个下标对应 y 方向，跟后面np.roll(..., axis=0)的心智模型对不上，是非常隐蔽的 bug 来源。我建议全程统一用"ij"。

外势取 V = |x-0.5|²/2 + |y-0.5|²/2，这样整个系统有解析解可以对拍，后面第 5 节要用。

3.2 投影步：加权泊松方程的离散与零空间处理

先把组装变系数算子的函数写出来。目标算子 A 满足 (Aφ)_i ≈ -∇·(ρ∇φ)：

from scipy.sparse import coo_matrix from scipy.sparse.linalg import cg def weighted_laplacian(rho, h): N, M = rho.shape idx = np.arange(N * M).reshape(N, M) nb = [ (np.roll(rho, -1, axis=0), np.roll(idx, -1, axis=0)), # +x 面 (np.roll(rho, +1, axis=0), np.roll(idx, +1, axis=0)), # -x 面 (np.roll(rho, -1, axis=1), np.roll(idx, -1, axis=1)), # +y 面 (np.roll(rho, +1, axis=1), np.roll(idx, +1, axis=1)), # -y 面 ] rows, cols, vals = [], [], [] diag = np.zeros((N, M)) for r_f, i_f in nb: w = 0.5 * (rho + r_f) / h**2 # 面上系数取算术平均 rows.append(idx.ravel()); cols.append(i_f.ravel()); vals.append(-w.ravel()) diag += w rows.append(idx.ravel()); cols.append(idx.ravel()); vals.append(diag.ravel()) return coo_matrix( (np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(N * M, N * M) ).tocsr()

注意符号：对角是加，非对角是减。这样组装出来 A 是半正定的，特征值都非负，正好满足 CG 的要求。如果符号搞反了，CG 会给出完全无意义的结果而且不报错——这是最难查的一类错误，建议写完先用常系数、常数右端手算两三个点对一下。

零空间的处理是关键细节。周期边界下 A 的零空间是常数向量，也就是说 φ 只定义到一个常数。这时候如果右端 b 的均值不为零，方程组根本无解，CG 会在残差里震荡最后返回一个垃圾。做法是：

def project_velocity(rho, wx, wy, h, rtol=1e-10): N, M = rho.shape A = weighted_laplacian(rho, h) px, py = rho * wx, rho * wy b = -((np.roll(px, -1, axis=0) - np.roll(px, 1, axis=0)) / (2 * h) + (np.roll(py, -1, axis=1) - np.roll(py, 1, axis=1)) / (2 * h)) b = b.ravel() b -= b.mean() # 关键：去掉均值 A = A + 1e-12 * coo_matrix((np.ones(N * M), (np.arange(N * M), np.arange(N * M)))).tocsr() phi, info = cg(A, b, rtol=rtol, maxiter=500) if info != 0: print(f"[warn] CG 未收敛, info={info}, 残差={np.linalg.norm(A @ phi - b):.3e}") phi = phi.reshape(N, M) phi -= phi.mean() vx = (np.roll(phi, -1, axis=0) - np.roll(phi, 1, axis=0)) / (2 * h) vy = (np.roll(phi, -1, axis=1) - np.roll(phi, 1, axis=1)) / (2 * h) return vx, vy, np.linalg.norm(A @ phi.ravel() - b)

那个1e-12的单位阵是刻意的——给零空间加一个极小的正位移，让矩阵严格正定，纯数值技巧，对结果的影响在双精度下看不出来。真正起作用的是b -= b.mean()那一行。

cg的参数名要小心：scipy 1.12 开始tol被改名为rtol，1.14 之后旧名字直接移除。老代码搬过来跑会直接 TypeError。用rtol是安全的。

3.3 输运步：通量形式迎风，保住质量

输运步要解的是一阶守恒律。用一阶迎风 + 显式欧拉，保守形式写出来是这样的：

def transport(rho, vx, vy, tau, h): vx_f = 0.5 * (vx + np.roll(vx, 1, axis=0)) # 面 (i-1/2, j) 上的速度 Fx = np.where(vx_f > 0, vx_f * np.roll(rho, 1, axis=0), vx_f * rho) vy_f = 0.5 * (vy + np.roll(vy, 1, axis=1)) Fy = np.where(vy_f > 0, vy_f * np.roll(rho, 1, axis=1), vy_f * rho) div = (np.roll(Fx, -1, axis=0) - Fx) / h + (np.roll(Fy, -1, axis=1) - Fy) / h return rho - tau * div

np.where那一行就是迎风的核心：面上速度为正时取左侧单元的值，为负时取右侧单元的值。写成这种形式之后，每个面的通量只在一个方向上取上游值，全局求和时内部通量两两抵消，剩下只有边界的净流出——周期边界下净流出为零，所以质量严格守恒到机器精度。

这里没有任何归一化操作。如果你发现质量在漂移，说明格式写错了，而不是"精度不够"，不要用重新归一化的方式去掩盖它。

另一条路是半拉格朗日格式，无条件稳定、可以取大时间步，但它不严格守恒，每次都要补质量，而且插值会引入抹平效应。这套方法里我更推荐迎风，简单、守恒、行为可预测。

3.4 主循环：把三步串起来并记录诊断量

def energy(rho, V, h): r = np.maximum(rho, 1e-14) return (r * np.log(r)).sum() * h * h + (V * r).sum() * h * h def run(rho, V, tau, n_steps, h, inner=1): for k in range(n_steps): # 1) 点态速度 mu = np.log(np.maximum(rho, 1e-14)) + 1.0 + V wx = -((np.roll(mu, -1, axis=0) - np.roll(mu, 1, axis=0)) / (2 * h)) wy = -((np.roll(mu, -1, axis=1) - np.roll(mu, 1, axis=1)) / (2 * h)) # 2) 投影：取梯度场分量 vx, vy, res = project_velocity(rho, wx, wy, h) # 3) 输运 rho = transport(rho, vx, vy, tau, h) neg = (rho < 0).sum() rho = np.maximum(rho, 0.0) rho /= rho.sum() * h * h if k % 50 == 0: print(f"step={k:5d} E={energy(rho, V, h):+.6f} " f"mass={rho.sum()*h*h:.12f} neg={neg} res={res:.2e}") return rho

neg这个计数器一定要加，它是你判断时间步是否合适的第一手指标。只要它非零，说明当前步长已经越界，后面所有结果都不可信。用1e-14而不是0去兜 log，是为了避免-inf直接污染整条链路。

inner参数在这版里没用上，但留个口子——如果你想在单个时间步里多迭代几次投影（把时间步内的分裂误差压下去），把它包一层循环就行。

4. 数值坑：负密度、质量漂移、时间步爆炸

4.1 负密度出现时，不要直接用clip糊弄过去

这是新手最容易养成坏习惯的地方。看到负密度，第一反应是rho = np.maximum(rho, 0)，跑通了，收工。但截断会在密度接近零的区域凭空造出或抹掉质量，而这类区域恰恰就是梯度流最活跃的地方——你等于在最关心的地方引入了最大的误差。

正确顺序是：先看neg的数量级。如果是零星几个单元，且截断掉的总质量占比低于 1e-10，那属于舍入误差范畴，截断后重新归一化无妨。如果neg占比超过 0.1%，或者截断的质量超过总量的 1e-6，那不叫误差，那叫格式失效，必须回头减小时间步。

更稳的做法有几条，按改造成本从低到高排：

第一条是减 τ。这是最直接的，成本只是多跑几步。第二条是把能量泛函从 ∫ρ log ρ 换成 ∫ρ log(ρ + ε)，用一个小 ε 把对数奇点正则化掉，代价是收敛到的方程不再是原来那个。第三条是改用通量限制器（比如在迎风通量上加一道 minmod 限制），能显著抑制振荡，但实现复杂度上一个台阶。

我的实际建议是：先用减 τ 把手上的问题跑通，再考虑换格式。绝大多数"必须改进格式"的场景，追下去都是时间步选太大了。

4.2 时间步与收敛性的双重约束

τ 到底怎么选？有两个约束同时起作用。

第一个是迎风格式的 CFL 条件：τ · max|v| / h ≤ 1，工程上留两倍余量，取 τ ≤ 0.5h / max|v|。麻烦的是 max|v| 不是常数——它近似等于 |∇log ρ|，密度分布越尖锐，这个值越大。一个初始方差 s0 = 0.01 的高斯，在 t=0 时刻的 max|v| 可以到几十，比 t 大到 1 的时候高一个量级。所以要么用自适应步长，要么按最坏情况取一个保守值。

第二个是投影步本身的平滑性质。这一版格式在数学上接近一种混合有限元离散，对高频模态有很强的阻尼，但显式推进仍然受限于网格。经验规则是 τ 不超过 h²/4 的量级，N=128 时 h≈0.0078，h²/4 ≈ 1.5e-5，比 CFL 给出的限制还紧。

这两个数摆一起就很明白了：真正卡住你的往往是 h² 那一项。N 翻倍，τ 要降到四分之一，总计算量涨十六倍。这也是为什么这类方法在大网格上必须上隐式或者半隐式，显式只适合验证和调参。

实践里我会先写一个自适应逻辑：

def pick_tau(vx, vy, h, safety=0.4): vmax = max(np.abs(vx).max(), np.abs(vy).max()) + 1e-12 return min(safety * h / vmax, 0.25 * h * h)

跑起来之后把每一步实际的 τ 打印出来，你会看到它在初期很小、后期逐渐放大，这本身就是一个很好的健康度指标。如果 τ 一直卡在极小值上不去，说明密度里存在持续尖锐的结构，需要回头看看是不是初值或者边界设置有物理上不合理的地方。

4.3 求解器静默失败：不报错但结果全错的那一类

比崩溃更可怕的是"跑完了、有数、但数是错的"。这类问题我遇到过三次，原因各不相同，但都属于同一类：迭代求解器返回了非零 info，代码却没检查。

第一次是忘了减右端均值，CG 在零空间上一直震荡，info返回 1（未收敛），但代码没看。第二次是密度在某片区域已经降到 1e-16 以下，加权泊松方程的系数矩阵在那里几乎退化，条件数爆到 1e18。第三次是maxiter设得太小（默认是 10*N，看着够用其实不够），残差还差三个量级就退出了。

防守手段很朴素但有效：每次调用 CG 都把info和最终残差打到日志里，并且设一个硬阈值——残差相对值超过 1e-6 就raise，让它当场爆出来，不要让它污染下游。宁可程序在第一百步挂掉，也不要它在第一千步给你一张看起来很漂亮的错误图。

另外提一句预条件。默认 CG 没有预条件，收敛速度对条件数很敏感。这个算子加个对角预条件（scipy.sparse.linalg.spilu或者简单用对角逆）通常能砍掉三到五成迭代次数，实现成本很低，值得做。

4.4 Windows上常见报错与真实原因的对照

顺带说一个不怎么像报错的问题：OpenBLAS 在某些 Windows 机器上会打印OpenBLAS : Program is Terminated. Because you tried to allocate too many memory regions。这通常是嵌套并行导致的——numpy 内部多线程和你的进程池撞在一起。设环境变量OMP_NUM_THREADS=1就能压下去。如果你用multiprocessing，记得把入口代码包在if __name__ == "__main__":里，Windows 用的是 spawn 而不是 fork，不加就是无限递归创建进程。

5. 怎么确认算出来的数是对的：三个自检手段

5.1 用高斯解析解做回归测试

没有解析解的数值代码是没法调试的，所以第一步永远是把问题收缩到一个有闭式解的特例上。

取外势 V = |x-c|²/2，能量为 E(ρ) = ∫ρ log ρ + ∫Vρ。对应的方程是 Fokker-Planck，若初值是各坐标方差为 s0 的高斯，那么解永远是高斯，方差按 s(t) = 1 + (s0 - 1