NIntegrate - почему в данном случае Mathematica 8 намного медленнее?

Ответ ruebenko и комментарии от user1091201 и Leonid вместе дают правильные ответы.

Ответ Редактировать 1 от ruebenko является первым правильным ответом для общих ситуаций, подобных этой, то есть добавьте параметр 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» работает медленнее, чем «MultidiversityRule» (хотя «LevinRule» дает более удовлетворительную оценку ошибки). "LevinRule" также медленнее, чем любое из нескольких одномерных правил типа Гаусса, повторяющихся в двух измерениях, таких как "GaussKronrodRule", которое нашел user1091201.

NIntegrate принимает решение после выполнения некоторого символьного анализа подынтегральной функции. Применяется несколько типов символьной предварительной обработки; настройка Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False} отключает один тип символической предварительной обработки.

По сути, при включенном «OscillatorySelection» NIntegrate выбирает «LevinRule». При отключенном «OscillatorySelection» NIntegrate выбирает «MultidirectionalRule», который быстрее для этого интеграла, хотя мы можем не доверять результату, основанному на сообщении NIntegrate::slwcon, которое указывает на то, что наблюдалась необычно медленная сходимость.

Это та часть, в которой Mathematica 8 отличается от Mathematica 7: Mathematica 8 добавляет «LevinRule» и соответствующую эвристику выбора метода в «OscillatorySelection».

Помимо того, что NIntegrate выбирает другой метод, отключение «OscillatorySelection» также экономит время, затрачиваемое на фактическую обработку символов, что в некоторых случаях может быть значительным.

Параметр Method -> "GaussKronrodRule" отменяет и пропускает символьную обработку, связанную с выбором метода, и вместо этого использует правило 2-D декартова произведения Method -> {"CartesianRule", Method -> {"GaussKronrodRule", "GaussKronrodRule"}}. Это оказывается очень быстрым методом для этого интеграла.

Другая обработка символов

И Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False} ruebenko, и Method -> "GaussKronrodRule" user1091201 не отключают другие формы символьной обработки, и в целом это хорошо. Список типов символьных предварительная обработка, которая может быть применена. В частности, «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 компактизирует подынтегральную функцию до конечного диапазона, эффективно выполняя выборку с геометрическим интервалом. Обычно это намного эффективнее, чем линейные (равномерно расположенные) выборки, которые используются для конечных диапазонов интегрирования.

Вот обходной путь:

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

Для этого конкретного интеграла главным виновником, по-видимому, является интегрирование по 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 достаточно "развязаны".