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 では 2 桁遅くなります。新しいバージョンでは、エラーをより適切に制御しようとしていると思いますが、この途方もない時間の増加を犠牲にしています. Mathematica 7 と同じ速度で計算を進めるために使用できる設定はありますか?

4

3 に答える 3

16

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」など、2 次元で反復されるいくつかのガウス型 1 次元ルールのいずれよりも低速です。

NIntegrate は、被積分関数の記号解析を実行した後に決定を下します。適用されるシンボリック前処理にはいくつかのタイプがあります。この設定 は、 Method -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}1 つのタイプのシンボリック前処理を無効にします。

基本的に、「OscillatorySelection」を有効にすると、NIntegrate は「LevinRule」を選択します。「OscillatorySelection」を無効にすると、NIntegrate はこの積分に対してより高速な「MultidimensionalRule」を選択しますが、異常に遅い収束が観察されたことを示すメッセージ NIntegrate::slwcon に基づく結果を信用できないかもしれません。

これが Mathematica 8 が Mathematica 7 と異なる部分です。Mathematica 8 は「LevinRule」と関連するメソッド選択ヒューリスティックを「OscillatorySelection」に追加します。

NIntegrate に別のメソッドを選択させるだけでなく、「OscillatorySelection」を無効にすると、実際のシンボリック処理に費やされる時間も節約されます。これは場合によっては重要です。

設定Method -> "GaussKronrodRule"はメソッド選択に関連するシンボリック処理をオーバーライドしてスキップし、代わりに 2-D デカルト積規則 を使用しMethod -> {"CartesianRule", Method -> {"GaussKronrodRule", "GaussKronrodRule"}}ます。これはたまたまこの積分の非常に高速な方法です。

その他のシンボリック処理

ruebenkoMethod -> {"SymbolicPreprocessing", "OscillatorySelection" -> False}と user1091201はどちらも他の形式のMethod -> "GaussKronrodRule"シンボリック処理を無効にしません。これは一般的には良いことです。適用できるシンボリック前処理のタイプのリストについては、NIntegrate の詳細ドキュメントのこの部分を参照してください。特に、"SymbolicPiecewiseSubdivision" は、区分関数が存在するためにいくつかの点で非解析的である被積分関数に対して非常に価値があります。

すべてのシンボリック処理を無効にし、既定のメソッド オプションを使用して既定のメソッドのみを取得するには、 を使用しますMethod -> {Automatic, "SymbolicProcessing" -> 0}。1 次元積分の場合、これは現在Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}、これらのメソッドのすべてのパラメーターの特定のデフォルト設定 (ルール内の補間点の数、グローバル適応戦略の特異点処理のタイプなど) に相当します。多次元積分の場合、現在のところ、これもMethod -> {"GlobalAdaptive", Method -> "MultidimensionalRule"}特定のデフォルト パラメータ値で になります。高次元の積分では、モンテカルロ法が使用されます。

NIntegrate の最初の最適化ステップとして直接切り替えることはお勧めしませんMethod -> {Automatic, "SymbolicProcessing" -> 0}が、場合によっては便利です。

最速の方法

ヒューリスティックに選択された非常に多くのパラメーターがあり、調整することで恩恵を受ける可能性があるため、数値積分を少なくとも少し、場合によっては大幅に高速化する方法はほぼ常にあります。( 「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}

これは非常に高速ではありませんが、かなり高速です。この手順は常にうまく機能するとは限りません (この手順から得られる 2D 積分グリッドが常に最適であるとは限らないため) が、被積分関数が積分を超え​​て十分に「分離」されている場合には十分に機能するはずxですy。

于 2011-12-10T14:40:13.277 に答える