R语言数值积分收敛异常排查与误差降低方法咨询

你遇到的积分结果异常，核心原因首先是被积函数本身的实现/模型设定问题，其次是使用通用数值积分方法没有针对被积函数“短时间收敛慢、长时间拖平缓尾巴”的特性做适配，之前收到的稀疏网格建议是高维积分场景的方案，一维积分完全不需要使用。

第一步：先排查模型与代码实现错误 #

这是结果超出0-1范围的核心诱因：

• 你当前prob1Dbox使用余弦级数，对应Neumann反射边界（粒子碰到边界反弹，永远留在区域内），不考虑突变项时总概率守恒，t→∞时函数值会收敛到常数invL，不会衰减到0。如果你的模型设定是粒子碰到边界就被移除（吸收边界），必须替换为正弦级数形式的Dirichlet热核，否则大t区间的长尾巴必然会在突变率很低时积出远大于1的结果。
• 核对级数的系数：你代码中指数项的系数写为c = pi * (pi/4) * (t * invtau)，和标准扩散热核的指数项exp(-n²π² D t / L²)对比多了1/4缩放因子，需要确认invtau的取值是否已经包含了4倍的系数补偿，否则系数错误会完全改变级数收敛速度。
• 注意t→0处的奇异性：反射边界热核在t→0时是狄拉克δ函数，temp_pbx * temp_pby会趋近于无穷大，被积函数在0附近的行为~1/t，积分本身是发散的。你在高突变率下得到有限值，本质是n_lim=1000的截断给小t区间的奇异性加了虚假的数值截断，结果只是碰巧落在0-1范围内，不代表计算正确。

第二步：数值积分方案调整（确认模型正确后使用） #

针对一维无穷区间积分，不需要稀疏网格，按以下方案调整即可实现收敛到机器精度：

• 拆分积分区间，分区域适配计算方法，避免通用积分的选点误差：
  1. 小t区间(0, t1]：取t1≈0.5（对应c≈0.001，级数收敛极慢），不要用截断级数计算prob1Dbox，直接用无限空间高斯核+镜像法计算：t极小时粒子还未碰到边界，只需要加最近2-3个镜像点的高斯项，精度远高于截断级数，也不会出现奇异性计算错误。
  2. 中间区间[t1, t2]：取t2≈10000（对应c≈20，n≥1的级数项衰减到2e-9以下，可以忽略），用你原本的级数求和计算即可，积分直接调用R的integrate函数做自适应Gauss-Kronrod求积，这段函数平滑，积分误差可以轻松控制到1e-10以下。
  3. 大t尾巴区间[t2, +∞)：这段t足够大，temp_pbx≈invL、temp_pby≈invL，被积函数可以简化为(1/(2*De)) * invL^2 * exp(-2*t*mut_rate)，直接用解析公式计算积分结果：(1/(2*De)) * invL^2 * exp(-2*t2*mut_rate) / (2*mut_rate)，完全没有数值误差。
• 优化级数求和逻辑：去掉固定n_lim=1000的循环，每计算一项判断项的绝对值，小于1e-15时直接终止循环，t较大时通常只需要前10-20项就满足精度要求，既提升速度也避免多余计算引入的误差。

验证方法 #

先用mut_rate=7e-3的参数测试调整后的方案，结果应该和你之前得到的0.001997469偏差小于1e-6，再运行低突变率的场景。如果调整后低突变率结果仍然大于1，直接回到第一步排查模型设定错误即可。

内容的提问来源于stack exchange，提问作者CafféSospeso