NIntegrate - 在给定的情况下,为什么Mathematica 8中的速度要慢得多？

Question

我有一个Mathematica代码,我必须在数值上评估数千个与此类似的积分

NIntegrate[
    (Pi*Cos[(Pi*(-2*x + y))/(1 + y)] + (1 + y)*(-Sin[(2*Pi*x)/(1 + y)] + 
    Sin[(Pi*(-2*x + y))/(1 + y)]))/(E^x*(1 + y)), 
    {x, 0, 100}, {y, 0, 100}
] //AbsoluteTiming

被积函数是一个很好的绝对可积函数,没有奇点,它在一个方向上呈指数衰减,在另一个方向上以1/y ^ 3衰减.

该NIntegrate命令在Mathematica 7中工作正常,但在最新版本8.0.4中,它减慢了两个数量级.我假设在新版本中它试图更好地控制错误,但代价是时间的巨大增加.是否有一些我可以使用的设置,以便计算以与Mathematica 7相同的速度进行？

Answer 1

ruebenko的回答以及user1091201和Leonid的评论结合在一起,给出了正确的答案.

对于像这样的一般情况,ruebenko的Edit 1回答是正确的第一个答案,即添加选项:Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}

expr = (Pi*
      Cos[(Pi*(-2*x + y))/(1 + y)] + (1 + y)*(-Sin[(2*Pi*x)/(1 + y)] +
         Sin[(Pi*(-2*x + y))/(1 + y)]))/(E^x*(1 + y));

NIntegrate[expr, {x, 0, 100}, {y, 0, 100}, 
  Method -> {"SymbolicPreprocessing", 
    "OscillatorySelection" -> False}] // AbsoluteTiming

而user1091201的评论建议Method -> "GaussKronrodRule"接近这个特定问题的最快答案.

我将在这个具体的例子中描述NIntegrate中正在发生的事情,并且一路上给出了处理明显类似情况的一些提示.

方法选择

在这个例子中,NIntegrate检验expr,得出的结论是多维"LevinRule"是这个被积函数的最佳方法,并应用它.然而,对于这个特定的例子,"LevinRule"比"MultidimensionalRule"慢(尽管"LevinRule"得到更令人满意的误差估计)."LevinRule"也比在两个维度上迭代的几个高斯型一维规则中的任何一个慢,例如user1091201找到的"GaussKronrodRule" .

NIntegrate在对被积函数进行一些符号分析后做出决定.应用了几种类型的符号预处理; 该设置 Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}禁用一种符号预处理.

基本上,启用"振荡选择"后,NIntegrate选择"LevinRule".在禁用"振荡选择"的情况下,NIntegrate选择"MultidimensionalRule",这对于此积分更快,尽管我们可能会根据消息NIntegrate :: slwcon不信任结果,这表明观察到异常缓慢的收敛.

这是Mathematica 8与Mathematica 7不同的部分:Mathematica 8将"LevinRule"和相关的方法选择启发法添加到"OscillatorySelection"中.

除了使NIntegrate选择不同的方法之外,禁用"OscillatorySelection"还可以节省执行实际符号处理所花费的时间,这在某些情况下可能很重要.

设置Method -> "GaussKronrodRule"覆盖并跳过与方法选择关联的符号处理,而是使用二维笛卡尔积规则Method -> {"CartesianRule", Method -> {"GaussKronrodRule", "GaussKronrodRule"}}.这恰好是这种积分的一种非常快速的方法.

其他符号处理

无论ruebenko的Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}和 user1091201的Method -> "GaussKronrodRule"不禁止其他形式的象征性的处理,而这通常是一件好事.有关可能应用的符号预处理类型的列表,请参阅NIntegrate高级文档的这一部分.特别地,"SymbolicPiecewiseSubdivision"对于由于存在分段函数而在若干点处非分析的被积分非常有价值.

要禁用所有符号处理并仅使用默认方法选项获取默认方法,请使用Method -> {Automatic, "SymbolicProcessing" -> 0}.对于一维积分,这当前相当于Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}这些方法的所有参数的某些默认设置(规则中的插值点的数量,全局自适应策略的奇点处理的类型等).对于多维积分,它目前相当于Method -> {"GlobalAdaptive", Method -> "MultidimensionalRule"}某些默认参数值.对于高维积分,将使用蒙特卡罗方法.

我不建议直接切换Method -> {Automatic, "SymbolicProcessing" -> 0}为NIntegrate的第一个优化步骤,但在某些情况下它可能很有用.

最快的方法

几乎总有一些方法可以加速数值积分至少一点,有时很多,因为有很多参数可以启发式选择,你可以从调整中受益.(查看"LevinRule"方法或"GlobalAdaptive"策略的不同选项和参数,包括所有子方法等)

也就是说,这是我为这个特殊积分找到的最快的方法:

NIntegrate[(Pi*
      Cos[(Pi*(-2*x + y))/(1 + y)] + (1 + y)*(-Sin[(2*Pi*x)/(1 + y)] +
         Sin[(Pi*(-2*x + y))/(1 + y)]))/(E^x*(1 + y)), {x, 0, 
   100}, {y, 0, 100}, 
  Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", 
    "SingularityDepth" -> Infinity}] // AbsoluteTiming

(该设置"SingularityDepth" -> Infinity禁用奇点处理转换.)

整合范围

顺便说一句,是你想要的积分范围真的{x, 0, 100}, {y, 0, 100},或者是{x, 0, Infinity}, {y, 0, Infinity}真正所需的积分范围为您的应用程序？

如果你真的需要{x, 0, Infinity}, {y, 0, Infinity},我建议改用它.对于无限长度积分范围,NIntegrate将被积函数压缩到有限范围,有效地以几何间隔方式对其进行采样.这通常比用于有限积分范围的线性(均匀)间隔样本更有效.

Answer 2

这是一个解决方法:

NIntegrate[(Pi*
      Cos[(Pi*(-2*x + y))/(1 + y)] + (1 + y)*(-Sin[(2*Pi*x)/(1 + y)] +
         Sin[(Pi*(-2*x + y))/(1 + y)]))/(E^x*(1 + y)), {x, 0, 
   100}, {y, 0, 100}, 
  Method -> "AdaptiveMonteCarlo"] // AbsoluteTiming

您还可以使用ParallelTry并行测试各种方法.

当实现新方法或修改启发式时,可能会发生特定参数的减速.这些可能有助于解决一类新问题,而其他一些问题则变得更慢.人们必须仔细调查这里发生了什么.

您可能想要更改问题的主题 - 它表示所有积分在V8中评估较慢,这是不正确的.

编辑1:我认为它被卡在LevinRule中(V8中的新振荡积分)所以,我认为,这应该关闭它.

NIntegrate[(Pi*
      Cos[(Pi*(-2*x + y))/(1 + y)] + (1 + y)*(-Sin[(2*Pi*x)/(1 + y)] +
         Sin[(Pi*(-2*x + y))/(1 + y)]))/(E^x*(1 + y)), {x, 0, 
   100}, {y, 0, 100}, 
  Method -> {"SymbolicPreprocessing", 
    "OscillatorySelection" -> False}] // AbsoluteTiming