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

Ответ 1

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

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

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

Это та часть, где Mathematica 8 отличается от Mathematica 7: Mathematica 8 добавляет "LevinRule" и связанные с ней методы выбора эвристики в "Осцилляторный выбор".

Помимо того, что 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"}, снова с определенными значениями параметров по умолчанию. Для высокоразмерных интегралов будет использоваться метод monte-carlo.

Я не рекомендую переключаться прямо на 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 компактифицирует подынтегральное выражение в конечном диапазоне, эффективно сэмплирует его геометрически разнесенным способом. Обычно это намного эффективнее линейных (равномерно) интервальных выборок, которые используются для конечных диапазонов интеграции.

Ответ 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

Ответ 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, [email protected], [email protected]}];

а затем имеем:

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

который не пылает быстро, но достаточно быстро. Эта процедура не всегда будет работать так хорошо (поскольку сетка интегрирования 2D, полученная в результате этой процедуры, не всегда будет оптимальной), но она должна работать достаточно хорошо, если подынтегральное выражение таково, что интегрирование по x и y развязаны".