有限差分的比热看上去合理，其实差 1.3%

主线位置：TRG 主线的最后两步「分段求导得内能 → 对照精确解」：从内能再求一次导得比热。

版本：JAX 0.12.0.1（CPU，x64）、tenferro-rs 0.8.0（CPU，native，单线程）

flowchart TB
  jx["JAX 用 jvp grad"] --> cj["c 精确到 AD"]
  port["tenferro 没有 二阶"] --> fd["对 AD 能量差分"]
  fd --> h5["h 1e-5"]
  h5 --> cross["跨过奇异值交叉"]
  cross --> bias["c 偏 1.3%"]
  fd --> h6["h 1e-6"]
  h6 --> ok["和 jvp grad 一致"]

实验代码：JAX 一侧是动手页的 ising_trg.py（jit(value_and_grad)，比热同时给 jvp(grad) 和对梯度的中心差分），tenferro 一侧是动手页的 ising_trg（eager 分段，逐步 vjp，release，CpuEngine::serial()）。χ = 16，20 步。M2 Max、macOS 27、rustc 1.99.0。

现象

比热 c = β² ∂²(lnZ/N)/∂β²。JAX 里它是一行 jax.jvp(jax.grad(f), (b,), (1.0,))，前向套反向，对整张图也不贵。tenferro 的分段程序只给一阶导数；二阶导数要么对整张图做（指数增长，见高阶导数），要么把每步的 vjp 再做 jvp。最省事的是对 AD 给出的能量做中心差分：

let (_, gp) = value_and_grad(&rt, &mut eng, b + h, &ks)?;
let (_, gm) = value_and_grad(&rt, &mut eng, b - h, &ks)?;
let c = b * b * (gp - gm) / (2.0 * h);

h = 1e-5，临界点附近三个 β：

# JAX，c 是对 grad 的差分，c_jvp 是 jvp(grad)
beta=0.4300 ... c=1.524985849603 c_jvp=1.524985776407 | ...
beta=0.4400 ... c=3.270066566479 c_jvp=3.311956762689 | ...
beta=0.4500 ... c=1.607084600399 c_jvp=1.607084461505 | ...
# tenferro
beta=0.4300 ... c=1.524985849596 | value_and_grad 562.3 ms
beta=0.4400 ... c=3.270066566430 | value_and_grad 570.1 ms
beta=0.4500 ... c=1.607084600412 | value_and_grad 559.9 ms
Onsager c [(0.43, ...(1.5268483322243886)), (0.44, ...(2.8896870923571734)), (0.45, ...(1.607952937640869))]

tenferro 的差分和 JAX 的差分对得上（差 1e-11 量级），看起来没有问题。
但 β = 0.44 时差分 3.2701 和 jvp(grad) 3.3120 差 1.3%；β = 0.43、0.45 只差 1e-7 量级。
和 Onsager 精确解比，β = 0.44 本来就差 +0.38 到 +0.42（χ = 16 的截断误差，临界点附近最大），1.3% 的偏差混在里面，单看 tenferro 的输出发现不了。

出处

截断 TRG 的 lnZ/N(β) 不是处处光滑的：最后一步保留 16 个奇异值，第 16 和第 17 个在 β ≈ 0.44000296 交叉。交叉时留下的那个奇异向量换了一个，lnZ/N 跳一下（相邻两个 f64 的 β 上差 1.5e-9），能量也随之跳。β ± 1e-5 正好跨过这个点，差分把跳变当成了斜率。AD 给的是跳变两侧各自光滑那一段的导数，jvp(grad) 不受影响。二分出交叉点、打印截断处奇异值的过程见 Ising 模型的 TRG 的比热一节。这是截断 TRG 本身的性质，和框架无关：JAX 的差分同样偏。

修法

用两个步长各算一次，对不上就说明中间有交叉，取小的那个，或者换一个 β：

# h = 1e-6
# JAX
beta=0.4300 ... c=1.524985783521 c_jvp=1.524985776407 | ...
beta=0.4400 ... c=3.311956764748 c_jvp=3.311956762689 | ...
beta=0.4500 ... c=1.607084498843 c_jvp=1.607084461505 | ...
# tenferro
beta=0.4300 ... c=1.524985783398 | value_and_grad 564.5 ms
beta=0.4400 ... c=3.311956765457 | value_and_grad 569.9 ms
beta=0.4500 ... c=1.607084498821 | value_and_grad 563.6 ms

h = 1e-6 时 β = 0.44 的差分和 jvp(grad) 差 2.8e-9，β = 0.43、0.45 的差也从 7e-8、1.4e-7 降到 7e-9、3.7e-8。h 也不能一直减小：梯度末位的噪声（换线程数就差 1e-12 量级）除以 2h，h = 1e-7 时已经开始主导（动手页的步长扫描）。

怎么确认没中招

每个 β 用 h 和 h/10 各算一次 c：正常时两者只差 1e-7 量级（上面 β = 0.43、0.45 是 6.6e-8、1.0e-7）；差到 1e-2 量级（β = 0.44 是 4.2e-2），就是跨过了交叉。
有 JAX 参考实现时，用 jvp(grad) 对照 tenferro 的差分，不要用 JAX 的差分对照 tenferro 的差分：两边的差分会一起错。

小结

tenferro 的分段程序没有便宜的二阶导数，比热靠对 AD 能量做差分；JAX 用 jvp(grad)。
β = 0.44、h = 1e-5 跨过了截断处第 16、17 个奇异值的交叉，c 偏 1.3%，两个框架的差分一起偏，和精确解的误差混在一起。
用两个步长交叉检查；对照时拿 jvp(grad) 当参考。

有限差分的比热看上去合理，其实差 1.3%

现象

出处

修法

怎么确认没中招

小结

💬 关联讨论

登录