10

我有一个 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 中相同的速度进行？

4

3 回答 3

16

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}禁用一种类型的符号预处理。

本质上，启用“OscillatorySelection”后，NIntegrate 会选择“LevinRule”。在禁用“OscillatorySelection”的情况下，NIntegrate 选择“MultidimensionalRule”，这对于这个积分来说更快，尽管我们可能不相信基于消息 NIntegrate::slwcon 的结果，这表明观察到异常缓慢的收敛。

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

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

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

其他符号处理

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

要禁用所有符号处理并仅获取具有默认方法选项的默认方法，请使用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 将被积函数压缩为有限范围，以几何间隔的方式对其进行有效采样。这通常比用于有限积分范围的线性（均匀）间隔样本更有效。

于 2011-12-10T19:58:32.897 回答
8

这是一种解决方法：

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
于 2011-12-10T14:04:30.707 回答
2

对于这个特殊的积分，罪魁祸首似乎是积分超过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}

这不是很快，但相当快。这个过程并不总是那么好（因为这个过程产生的二维积分网格并不总是最优的），但是当被积函数使得积分结束x并且y充分“解耦”时应该工作得很好。

于 2011-12-10T14:40:13.277 回答