reference

方法:离散化、抽样与收敛

一份带编号方程的方法学处理:主平衡、代表性切片约化、双重积分的求积离散化、带 1/√N 误差的蒙特卡洛估计量,以及相对方差的传播律,并讨论收敛。

最后更新: 2026-07-26

本章以方法学章节的风格,阐述从 Schubart Master Formula 中提取一个数字所依凭的计算程序。它从连续的主平衡出发,分五个带编号的步骤走向其实际求值:精确积分、其向代表性切片的约化、其作为求积和的离散化、其由蒙特卡洛的估计,以及输入不确定度的传播。每条关键关系都加以编号,以便上下文和后续各章能无歧义地将其称为”式 (n)“。为这些步骤提供动机的完整解析推导见完整推导;此处的重点是积分究竟如何被计算,以及每种方法所带的误差。

方程 (1) —— 扩展的平衡

出发点是 Schubart Master Formula 的扩展的物理表示——以其谱–输运形式写出的、关于瞬时转换功率的核算恒等式,

P(t) = \iint \Phi(E,\theta,t)\,\sigma(E)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;dE\,d\theta \tag{1}

其中积分遍及入射能量 $E$ 和角度坐标 $\theta$ ,五个因子分别是通量密度、散射截面、有效动量转移、声子速度和转换效率。式 (1) 是在二维域上的双重积分;本章随后的每种方法都是求值它的一种方式。表明右端化简为瓦特的量纲核算在量纲分析中进行,此处假定已知。

方程 (2) —— 代表性切片

当被积函数在一个精选窗口 $\Delta E \times \Delta\theta$ 上缓变,或当需要一个初步估计时,式 (1) 约化为单个代表性切片——以在代表值处求得的被积函数乘以窗口大小来取代积分,

P \approx \Phi\,\sigma\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E\,\Delta\theta . \tag{2}

式 (2) 是式 (1) 的中值形式:若各因子在窗口上为常数则精确,否则为一个可控的近似。它的有效性恰是假设与局限中假设 C 的代表性平均条件——窗口必须宽到足以具代表性,又窄到使被积函数不在其上剧烈摆动。式 (2) 对数量级估计以及核查更完整的计算是否出偏极为有用,但在被积函数尖锐成峰之处——例如横跨一个截面共振——绝不可信任它。

方程 (3) —— 求积离散化

为在不假设单个代表性切片的情况下求值式 (1),将域划分为以 $j$ (能量)和 $k$ (角度)为指标的单元网格,并以在这些单元上的黎曼或求积和来逼近积分,

P \approx \sum_{j,k} \Phi_{jk}\,\sigma_j\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E_j\,\Delta\theta_k , \tag{3}

其中 $\Phi_{jk}$ 和 $\sigma_j$ 是在单元中心抽样的因子值, $\Delta E_j$ 、 $\Delta\theta_k$ 是单元宽度。式 (3) 不过是式 (2) 对许多切片求和;当单元收缩时它收敛于式 (1)。如所写的均匀中点法则,每单元的离散化误差在被积函数的二阶导数上为 $(\Delta E)^3$ 和 $(\Delta\theta)^3$ 的量级,故中点或梯形法则的全局误差按单元尺寸 $h$ 标度为 $O(h^2)$ ;辛普森法则在光滑被积函数上将其改进为 $O(h^4)$ 。把单元更密地置于被积函数变化最快之处——譬如围绕一个谱峰——的非均匀网格,远比均匀加密更廉价地恢复精度。这种基于网格的求值是确定性数值模拟的骨干。

方程 (4) —— 蒙特卡洛估计量

当被积函数高维、尖锐成峰或形状别扭时,网格求积变得昂贵,因为单元数随维数呈幂律增长。蒙特卡洛的替代办法把积分改写为一个期望,并以随机抽样来估计它。从域中抽取 $N$ 个点 $x_n = (E_n,\theta_n)$ ,并记 $f(x)$ 为式 (1) 的完整被积函数乘以域体积,则估计量是样本均值

\hat P = \frac{1}{N}\sum_{n=1}^{N} f(x_n), \qquad \mathrm{SE}\big[\hat P\big] \propto \frac{1}{\sqrt{N}} . \tag{4}

式 (4) 是无偏的——其期望等于真实积分——且其标准误差按 $1/\sqrt{N}$ 下降,与域的维数无关。正是这种维数无关性,使蒙特卡洛在复杂被积函数上超越张量网格求积,也正是该方法自然推广到不确定度传播的原因:通过在一个外层循环中让输入因子本身从其不确定度分布中抽取,每一趟蒙特卡洛就给出 $P$ 的一个实现,而整个系综再现了输出分布。完整的嵌套程序见蒙特卡洛验证。

方程 (5) —— 不确定度传播

最后一种方法关心的不是积分的值,而是它的不确定度。因为式 (1) 是各因子之积,各因子的相对方差相加。对于具统计独立因子的 $P = \prod_i x_i$ ,一阶传播给出

\left(\frac{\sigma_P}{P}\right)^{2} = \sum_i \left(\frac{\sigma_{x_i}}{x_i}\right)^{2} . \tag{5}

式 (5) 由对数微分 $\mathrm{d}\ln P = \sum_i \mathrm{d}\ln x_i$ 得出,并连同相关输入的协方差项在不确定度分析中完整导出。它在诸方法中的作用,是在任何昂贵抽样之前告诉我们哪个输入主导输出误差:因为各项以平方相加,最大的相对不确定度控制结果,精化的力气应投向那里。在小不确定度极限下,式 (4) 的蒙特卡洛散布再现式 (5),故二者互为一致性检验。

收敛

三种求值方法——式 (2)、(3) 和 (4)——带有性质不同的收敛行为,在它们之间抉择,是把方法与被积函数相匹配的问题。

代表性切片,式 (2),在任何精化意义上都不收敛;它是一个固定的近似,其误差完全由所选窗口有多具代表性来设定。它是一个初步估计,而非一个收敛的方案。

求积,式 (3),按单元尺寸的幂收敛:在光滑被积函数上,中点与梯形法则为 $O(h^2)$ ,辛普森法则为 $O(h^4)$ 。因此单元尺寸减半会使梯形误差减为四分之一。这种快速的高阶收敛,使求积成为光滑、低维被积函数的首选方法——正是二维式 (1) 在各因子行为良好时的情形。它的弱点在于单元数随维数呈幂律增长,且尖锐的峰迫使用极细的网格。

蒙特卡洛,式 (4),按 $1/\sqrt{N}$ 收敛:将误差降到十分之一需要样本增加百倍。这在低维下比求积慢,但指数 $-\tfrac12$ 与维数无关,故对于求积网格爆炸的高维或尖锐结构被积函数,蒙特卡洛胜出。它缓慢的速率通常靠方差缩减来改善:重要性抽样,它把抽样集中于被积函数最大之处并以 $f/g$ 重新加权;分层抽样,它保证覆盖重要的子区域;以及控制变量,它减去一个积分已知的相关量。这每一种都在不改变 $1/\sqrt{N}$ 标度的情况下降低标准误差分子中的有效方差 $\mathrm{Var}[f]$ ,使一个适中的 $N$ 就能达成否则需要多出数量级的均匀样本才能达成的效果。

实践中,对式 (1) 的稳健求值会用到三者:式 (2) 用于做合理性核查的数量级,式 (3) 用于被积函数光滑时的确定性值,式 (4) 加方差缩减用于不确定度传播以及网格变得不切实际的任何情形。随后,传播律式 (5) 提供抽样预期应当再现的解析误差预算。