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 的编辑 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 在对被积函数执行一些符号分析后做出决定。应用了多种类型的符号预处理；设置方法-> {“SymbolicPreprocessing”，“OscillatorySelection”-> False} 禁用一种类型的符号预处理。

本质上，启用“OscillatorySelection”后，NIntegrate 选择“LevinRule”。禁用“OscillatorySelection”后，NIntegrate 选择“MultiDimensionalRule”，该积分速度更快，尽管我们可能不信任基于消息 NIntegrate::slwcon 的结果，该消息表明观察到异常缓慢的收敛。

这是 Mathematica 8 与 Mathematica 7 不同的部分：Mathematica 8 在“OscillatorySelection”中添加了“LevinRule”和相关的方法选择启发式。

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

设置方法-> “GaussKronrodRule”覆盖并跳过与方法选择相关的符号处理，而是使用二维笛卡尔积规则Method -> {“笛卡尔规则”，方法-> {“GaussKronrodRule”，“GaussKronrodRule”}}。对于这个积分来说，这恰好是一种非常快速的方法。

其他符号处理

两个ruebenko的方法 -> {“SymbolicPreprocessing”，“OscillatorySelection”-> False} 和 user1091201 的 Method -> “GaussKronrodRule” 不会禁用其他形式的符号处理，这通常是一件好事。有关符号类型的列表，请参阅NIntegrate 高级文档的这一部分可以应用的预处理。特别是，“SymbolicPiecewiseSubdivision”对于由于分段函数的存在而在多个点上不可解析的被积函数非常有价值。

要禁用所有符号处理并仅获取具有默认方法选项的默认方法，请使用Method ->; {自动，“SymbolicProcessing”-> 0}。对于一维积分，目前相当于方法 -> {“GlobalAdaptive”，方法-> "GaussKronrodRule"} 具有这些方法的所有参数的某些默认设置（规则中的插值点数量、全局自适应策略的奇点处理类型等）。对于多维积分，目前相当于Method -> {“GlobalAdaptive”，方法-> “MultiDimensionalRule”}，同样具有某些默认参数值。对于高维积分，将使用蒙特卡罗方法。

我不建议直接切换到 Method -> {自动，“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

Answer 3

对于这个特定的积分，罪魁祸首似乎是对 x 的积分，这可能是由于快速衰减项和高振荡项的存在。另外，对于这种情况，可以通过分析方式对 x 进行积分：（

In[92]:= 
-Integrate[(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]/.x->0//FullSimplify

Out[92]= (-\[Pi] (1+y) (2+Cos[(\[Pi] y)/(1+y)])+(2 \[Pi]^2+(1+y)^2) Sin[(\[Pi] y)/(1+y)])/
(4 \[Pi]^2+(1+y)^2)

我丢弃了上限的值，因为它对于 y 来说都非常小）。然后可以对 y 进行数字积分以获得正确的结果：

In[94]:= NIntegrate[%,{y,0,100}]
Out[94]= 1.17525

一种更通用的解决方法是拆分 x 和 y 积分，如下所示：

ClearAll[f];
f[y_?NumericQ, xrange : {_?NumericQ, _?NumericQ}] :=
  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, First@xrange, Last@xrange}];

然后我们就得到了：

In[105]:= NIntegrate[f[y,{0,100}],{y,0,100}]//Timing
Out[105]= {2.578,1.17525}

这不是极快的速度，但相当快。此过程并不总是那么有效（因为由此过程产生的 2D 积分网格并不总是最优的），但当被积函数使得 x 和 < code>y 已充分“解耦”。