validation

数值模拟

数值求解 P(t) 主公式的方法论:将对能量和角度的积分离散化、提供输入,并检验收敛性。

最后更新:

由于 Schubart Master Formula 是一个对能量和角度的二重积分,针对任何具体案例求解它都需要数值方法。本章描述做这件事的一般方法论——如何将积分离散化、如何提供输入,以及如何确认计算结果已经收敛。全章的重点在于方法而非具体数字;凡提到的任何数值都是示意性的占位符,并以此身份呈现,而非测量结果。

将积分离散化

出发点是扩展的物理表示,

P(t)=∬Φ(E,θ,t) σ(E) Δp vph ηconv  dE dθ.P(t) = \iint \Phi(E,\theta,t)\,\sigma(E)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;dE\,d\theta .

最直接的数值方法是用网格上的求和代替连续积分。能量轴被划分为 NEN_E 个区间、角度轴被划分为 NθN_\theta 个区间,被积函数在每个单元的代表点处求值。将被积函数写作 f(E,θ,t)f(E,\theta,t),一个简单的求积式为

P(t)≈∑i=1NE∑j=1Nθwij f(Ei,θj,t) ΔEi Δθj,P(t) \approx \sum_{i=1}^{N_E}\sum_{j=1}^{N_\theta} w_{ij}\,f(E_i,\theta_j,t)\,\Delta E_i\,\Delta\theta_j ,

其中 ΔEi\Delta E_i 和 Δθj\Delta\theta_j 是单元的宽度,wijw_{ij} 是求积权重。对于中点法则,权重为 1,采样点位于单元中心。更精确的格式——梯形、辛普森或高斯-勒让德求积——用加权求和代替中点求值,从而在给定点数下减小离散化误差。格式的选择是精度与计算代价之间的标准权衡,而框架本身不偏好任何特定的求积方法。

离散化策略与收敛阶

由于定义域是能量轴与角度轴的乘积,二维法则最简单的构造方式是把两个一维法则作张量积,一个用于 EE,一个用于 θ\theta。这种因子化对代价和精度都有影响。在代价方面,张量网格使用 NE×NθN_E \times N_\theta 次被积函数求值,因此同时细化两个轴的代价是二次增长的;当所有结构都集中在一个轴上时——比如 σ(E)\sigma(E) 沿能量轴有一个尖锐的共振——单独细化那个轴往往便宜得多。在精度方面,每个一维法则都有自己的收敛阶,即渐近误差标度中的指数 pp

∣Pexact−PN∣  ∼  C N−p,\left| P_{\text{exact}} - P_{N} \right| \;\sim\; C\,N^{-p} ,

其中 NN 是沿该轴的点数,CC 取决于被积函数的光滑程度。对于梯形法则 p=2p = 2,因此把步长减半大致可将误差减小到四分之一;对于辛普森法则,在光滑被积函数上 p=4p = 4;而对于高斯-勒让德求积,收敛更快,对定义域内无奇点的解析被积函数实际上是指数收敛。这些阶数正是 Press 等人《Numerical Recipes》等标准参考书中列出的经典结果,它们通过张量构造延续到二重积分。实际后果是:决定可达精度的是被积函数的光滑程度,而非格式的名义阶数:一个四阶法则若跨越 σ(E)\sigma(E) 中近乎不连续处,会在局部退化为一阶行为,全局误差便由那一个未被妥善解析的区域所主导。

提供输入

一次模拟的意义不会超过其输入的意义。必须提供四类输入。第一,通量模型 Φ(E,θ,t)\Phi(E,\theta,t):对环境谱及其角分布的描述,可以是实测的,也可以取自假定的分布。第二,截面 σ(E)\sigma(E):与能量相关的耦合,通常以表格化或参数化的函数形式提供;其解析行为在 sigma 分析中讨论。第三,动量转移、输运和效率参数 Δp\Delta p、vphv_{\text{ph}} 和 ηconv\eta_{\text{conv}},它们可以是常数,也可以本身依赖于 EE 和 θ\theta。所有这些量的定义和单位固定于变量定义中。由于这些输入中有几个带有不确定性,一次负责任的模拟应传播该不确定性,而不是报告单一的确定性数字——这一点在蒙特卡洛验证中展开。

检验收敛性

确定性求积的核心质量控制问题,是结果是否已就网格分辨率而言收敛。标准做法是在一系列逐步细化的网格上计算 P(t)P(t)——每一步都把 NEN_E 和 NθN_\theta 加倍——并观察估计值如何变化。当连续的细化使结果变化小于选定的容差时,网格就被认为足够。这一检验的更定量版本直接使用假定的收敛阶:如果误差按 N−pN^{-p} 标度,那么在步长为 hh、h/2h/2 和 h/4h/4 的三个网格上计算积分,并比较相继的差值,就能揭示观测到的阶数,并通过理查森外推给出改进的估计。将三个结果写作 PhP_h、Ph/2P_{h/2}、Ph/4P_{h/4},比值

Ph−Ph/2Ph/2−Ph/4  ≈  2p\frac{P_{h} - P_{h/2}}{P_{h/2} - P_{h/4}} \;\approx\; 2^{p}

应当复现预期的阶数;若某值远离 2p2^p,便是一个警告,说明被积函数不如假设的那样光滑,名义误差估计不可信赖。对于光滑的被积函数,一个选得好的求积会快速收敛;σ(E)\sigma(E) 或通量中的尖锐特征,例如一个窄共振,要求在该特征区域更细的采样,是收敛缓慢的常见原因。自适应格式会在局部误差估计超过阈值处细分单元,从而在被积函数变化快的地方布置更多点,是自然的补救办法。

仅作为示意性示例,人们可能发现把能量网格从粗设置加密到中等精细会使某个占位输出估计值移动百分之几,此后进一步细化使其变化远小于百分之一——这是计算已收敛的标志。这类数字纯属概念性,对任何真实案例都需重新确立;此处给出它们只是为了展示收敛研究的形状,而非断言某个结果。

误差估计与预算编制

一个站得住脚的数值结果应配有一份明确的误差预算,将各个不同的贡献分开。离散化误差就是上文的 C N−pC\,N^{-p} 项,实践中由相继细化的差值 ∣Ph/2−Ph/4∣|P_{h/2} - P_{h/4}| 估计,对一个收敛的序列而言,这将剩余误差限制在一个小因子之内。截断误差来自在有限上限处切断能量定义域;它由被忽略尾部的积分所界定,定义域应一直延伸到该界限落到离散化误差之下为止。输入误差由 Φ\Phi、σ(E)\sigma(E)、Δp\Delta p、vphv_{\text{ph}} 和 ηconv\eta_{\text{conv}} 的不确定性传播而来,通常是三者中最大的,应由蒙特卡洛验证的随机方法处理,而非靠网格细化。一个常见且有教益的错误,是把网格细化到机器精度,却忽略一个高达数十个百分点的输入不确定性——打磨一个其真实误差由求积无法触及的因子所主导的数字。纪律在于报告全部三项贡献,并让最大的那一项决定每次进一步细化值得投入多少精力。

数值误差的来源

除离散化误差之外,谨慎的实践者还会留意几个反复出现的陷阱。过早截断能量积分会丢弃尾部的通量;定义域应延伸到被积函数真正可忽略为止,如数学架构中所述。使用一个忽略了正确测度的角度网格——例如在立体角积分需要 sin⁡θ\sin\theta 因子时将其遗漏——会引入任何细化都无法消除的系统性偏差,因为它是被积函数本身的误差,而非其采样的误差。而在跨越尖锐特征处对粗糙表格化的 σ(E)\sigma(E) 进行插值,可能低估或高估它与通量的重叠。上述每一项都可通过标准的数值卫生加以避免,每一项都应与任何报告的结果一并记录。

从确定性求值到统计求值

确定性求积回答的问题是”给定这些精确的输入,P(t)P(t) 是多少?“它本身并不回答”给定这些输入的不确定性,我们对 P(t)P(t) 应有多大信心?“第二个问题需要随机方法:从输入各自的分布中抽样,并以统计方式估计积分。这些方法及其产生的置信区间,是蒙特卡洛验证和 sigma 分析的主题。一个收敛的确定性求积在那里还承担第二个角色:在单一固定的输入集上运行时,它是蒙特卡洛估计器据以交叉检验的参照,因此两种方法是互补的而非彼此替代的。

延伸阅读

  • Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press——求积格式、收敛阶与自适应积分的标准参考书。
  • Stoer, J., & Bulirsch, R. (2002). Introduction to Numerical Analysis (3rd ed.). Springer——对求积误差分析与理查森外推的严格论述。
  • Davis, P. J., & Rabinowitz, P. (1984). Methods of Numerical Integration (2nd ed.). Academic Press——关于高斯求积与多维积分的全面论述。
↑