NIntegrate積分則

はじめに

積分則は,通常重み付き総和を使って領域上の積分推定値を計算する.NIntegrateが使用されるコンテキストでは,積分則オブジェクトは積分推定値,および積分推定値の確度の測定としての誤差推定の両方を与える.

次の一般的な説明は"GaussKronrodRule""MultidimensionalRule"等の重み付き総和型の積分則に当てはまる."LevinRule"等,他の型の規則には別の解説が必要である.

積分則はサンプル点と呼ばれる点の集合における被積分関数をサンプリングする.文献ではこれらは分点とも呼ばれる.それぞれのサンプル点 に対応して重みの数 がある.積分則は重み付き総和 で積分 を推定する.積分則は関数的であり,区間(あるいはより一般的な領域)上の関数を実数にマップする.

規則が領域 に適用されると,これは は被積分関数)と表される.

以下で考える規則のサンプル点は区間,単位立方,中心化された単位立方 は積分の次元)のいずれかにおける積分の推定値を計算するために選ばれている.よって, がこれらの領域のいずれかであれば, の代りに が使われる.これらの規則が他の領域に適用されると,その横座標と推定値はそれに応じてスケールされる.

積分則 は, ならば関数 に対して厳密であると言われる.

積分則 を関数 に適用することは, の積分といわれる.つまり, で積分されると, を得るのである.

一次元の積分則で,次数 以下の多項式はすべて厳密に積分するが次数 の多項式のうちの少なくとも1つが失敗する場合,その規則は次数 と言われる.

多次元の積分則で,次数 以下の単項式はすべて厳密に積分するが,次数 の単項式の少なくとも1つが失敗する場合,つまり規則が形式 は次元,)の単項式すべてについてすべて厳密である場合,その規則は次数 と言われる.

次数 のヌル規則は,次数 の単項式すべてについてゼロに積分され,次数 の単項式の少なくともひとつについては失敗する.それぞれのヌル規則は,基本的な積分則と低次元の適切な積分則との差分と考えることができる.

次数 の規則 のサンプル点の集合に低次元の )の規則 のサンプル点の集合が含まれる場合,に埋め込まれているという.これは と表される.

一般的な微分およびプロパティを持つが次数が異なる規則族のメンバーである次数 の積分則は, と記述される.ここで は族を識別するように選ばれる(例えば,次数4の台形規則は と参照される).

1つの族のそれぞれの規則が同じ族の別の規則に埋め込まれている場合,その族の規則は進行的であるという(ある に対して が存在しそれについて が成り立つ).

積分則は,被積分関数が区間の終点で評価されないならば「開いたタイプ」と呼ばれる.

NIntegrateの積分則オブジェクトは,積分推定に対して1つの積分則と誤差推定に対して1つあるいは複数のヌル規則を持つ.積分則とヌル規則のサンプル点は一致する.「積分則」もしくは「規則」がNIntegrateの積分則オブジェクトを意味するのか,通常の数学的な意味での積分則を意味するのかは,コンテキストから分かるはずである.

積分則の指定

以下で説明する以外の積分則はすべて,NIntegrate適応型ストラテジーで使用される.NIntegrateでは,モンテカルロ法は単純でも適応型でもを使う.

積分ストラテジーの積分則要素を変更すると,異なる積分法になる.

NIntegrateの適応型ストラテジー(「大域的適応ストラテジー」「局所的適応ストラテジー」を参照)で使用する積分則を指定するためには,Methodサブオプションを使う.

ストラテジーで積分則を使った例である.
In[1]:=
Click for copyable input
Out[1]//InputForm=
同じ積分則を異なるストラテジー()で使ってみる.
In[2]:=
Click for copyable input
Out[2]//InputForm=

NIntegrateに積分則指定(以外)だけを持つメソッドオプションが与えられると,その規則はストラテジーで使われる.以下の2つの入力は同等である.

この積分では,積分則だけが指定されている.デフォルトでストラテジーが使われる.
In[74]:=
Click for copyable input
Out[74]//InputForm=
以下の積分では,積分ストラテジーと積分則が指定されている.
In[75]:=
Click for copyable input
Out[75]//InputForm=

では,適応型モンテカルロストラテジーが自動的に使用される.下の2つのコマンドは等価である.

以下のモンテカルロ積分では,だけが指定されている.
In[5]:=
Click for copyable input
Out[5]//InputForm=
このモンテカルロ積分では,モンテカルロ積分ストラテジーとが指定されている.
In[6]:=
Click for copyable input
Out[6]//InputForm=

"TrapezoidalRule"

積分推定における台形規則は最も簡単で昔からある規則の一つである(古代ギリシャの数学者により使われ,バビロン人にも使われた可能性がある).

複合台形規則は以下の形式

のリーマン総和である.ここで である.

Methodオプションに値が与えられると,複合台形規則は積分則により形成された各部分区間を推定するために使われる.

積分.
In[7]:=
Click for copyable input
Out[7]=
オプション名
デフォルト値
"Points"5粗い台形の点の数
"RombergQuadrature"Trueロンバーグ(Romberg)積分法を使うかどうか
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

台形規則およびその複合(マルチパネル)拡張はあまり正確ではない.複合台形規則は線形関数については厳密であり,被積分関数が連続の第二導関数を持つなら少なくとも の速さで収束する[DavRab84].マルチパネル台形規則の正確さはロンバーグ求積法を使うと増加させることができる.

の横軸は の部分集合であるため,差分は積分推定値 の誤差推定を取ることができ,追加の被積分関数の評価をしないで計算することができる.

オプション"Points"->k はいくつの粗い点を使うかを指定するために使うことができる.で使用される点の合計数はである.

以下で,サンプル点が(1)のようになっているかを検証する.
In[8]:=
Click for copyable input
Out[9]=

は多次元積分にも使える.

これはを使った多次元積分である.厳密解はである.
In[10]:=
Click for copyable input
Out[10]=

注意:NIntegrateには台形規則と台形ストラテジーの両方がある.「積分ストラテジー」チュートリアルの「台形ストラテジー」を参照のこと.NIntegrateの内部的に実装された積分則はすべて-Ruleという接尾辞を持つ.よって,は台形積分則を指定するために使われ,は積分ストラテジーを指定するために使われる.

ロンバーグ積分法

ロンバーグ積分法とは,の打ち切り近似誤差の同じ次数項を除去する の線形結合を使うものである.

オイラー・マクローリン(EulerMaclaurin)公式[DavRab84]から

が得られる.ここで

が成り立つ.よって,以下のように書くことができる.

上の方程式の 項は,2つ目の方程式から最初の方程式を4回引くと除去することができる.結果は下のようになる.

次は,ロンバーグ積分法を使う台形規則の方が,標準の台形規則よりもよいパフォーマンスを見せる例である.前者の結果の方が厳密値 に近い.
In[11]:=
Click for copyable input
Out[11]=
次はロンバーグ積分法を使わない台形規則を使った積分である.
In[10]:=
Click for copyable input
Out[10]=

"TrapezoidalRule"のサンプル点と重み

次は指定された精度での台形サンプル点,重み,誤差の重みを計算する.
In[3]:=
Click for copyable input
Out[4]=
以下はロンバーグ積分法の重みと誤差の重みの導出方法である.
In[5]:=
Click for copyable input
Out[9]=

"NewtonCotesRule"

ニュートン・コーツ(NewtonCotes)積分公式は,等間隔のサンプル点を持つ補間タイプの公式である.

NIntegrateでニュートン・コーツ法を指定するためにはMethodオプション値を使う.
In[20]:=
Click for copyable input
Out[20]=
オプション名
デフォルト値
"Points"3粗いニュートン・コーツ点の数
"Type"Closedニュートン・コーツ規則のタイプ
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

積分区間 を以下の点

で,等間隔の 個の部分区間に分けるとする.すると,補間タイプの積分公式は

で与えられる.ここで

が成り立ち,

である. が大きいとき,ニュートン・コーツの 点係数は大きく,符号も混在する.

In[21]:=
Click for copyable input
Out[21]=

これでは桁落ちにより多数の有効数字が失われる可能性があるため,高次のニュートン・コーツ規則は注意して使う必要がある.

"NewtonCotesRule"のサンプル点と重み

次は指定された精度での台形サンプル点,重み,誤差の重みを計算する.
In[22]:=
Click for copyable input
Out[23]=

"GaussBerntsenEspelidRule"

ガウス積分法では,最適サンプル点を使い(多項式補間による),これらの点上での被積分関数値の重み付き総和を形成する.これらのサンプル点の部分集合上で,低次数の積分則が作られる.2つの規則の差分は誤差を推定するために使うことができる.BerntsenとEspelidは,奇数のサンプル点を使ってガウス規則の中点を削除することにより誤差推定規則を導き出した.

NIntegrateではガウス積分法はMethodオプション値で指定できる.
In[24]:=
Click for copyable input
Out[24]=
オプション名
デフォルト値
"Points"Automaticガウス点の数
"SymbolicProcessing"Automatic記号的前処理を行う秒数

オプション

被積分関数 に対する 点のガウス規則 は,次数の多項式の場合厳密である(つまり, が次数の多項式ならば である).

ガウス規則では被積分関数が区間の終点で評価されないので,開いたタイプである(Lobatto規則クレンショウ・カーチス則台形規則は区間の終点で積分評価を使うので閉じたタイプである).

差分商関数[Ehrich2000]を定義する.

ガウス規則 のサンプル点が であったら,BerntsenとEspelidは以下の誤差推定関数を導いていたことになる([Ehrich2000]を参照).

([Ehrich2000]のもとの公式はのサンプル点のためのものである.上の公式はのサンプル点のためのものである.)

のオプションをさまざまな値にしてNIntegrateで使用されるサンプル点の数を示す.
In[25]:=
Click for copyable input
Out[25]=

"GaussBerntsenEspelidRule"のサンプル点と重み

以下は,指定された数の粗い点と精度でガウスの横座標,重み,BerntsenEspelidの誤差の重みを計算する.
In[26]:=
Click for copyable input
Out[27]=

Berntsen-Espelidの誤差の重みは以下で実装される.

これで差分商を実装する.
In[28]:=
Click for copyable input
の横座標と重みを計算する.
In[30]:=
Click for copyable input
BerntsenEspelidの誤差の重みを計算する.
In[31]:=
Click for copyable input
Out[31]=

"GaussKronrodRule"

ガウス積分法は最適サンプル点(多項式補間による)を使い,これらの点上の被積分関数値の重み付き総和を形成する.ガウス規則のクロンロッド拡張は,ガウス点の間に新しいサンプル点を加え,ガウス規則の被積分関数評価を再利用する高次の規則を形成する.

NIntegrateでガウス・クロンロッド法を指定するためには,Methodオプション値を使う.
In[32]:=
Click for copyable input
Out[32]=
オプション名
デフォルト値
"Points"Automaticクロンロッド点で拡張されるガウス点の数
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

被積分関数 に対する 点のガウス規則 は,次数の多項式の場合厳密である(つまり, が次数の多項式ならば である).

ガウス・クロンロッド規則では被積分関数は区間の終点で評価されないので,開いたタイプである.

点のガウス規則 のクロンロッド拡張 は,個の点を加える.拡張された規則は が偶数の場合次数の多項式に対して, が奇数の場合はの多項式に対して厳密である.ガウス規則に関連付けられた重みはそのクロンロッド拡張で変化する.

の横座標は の部分集合なので,差分 は積分誤差 の誤差推定とされ,余分な被積分関数の評価なしで計算することができる.

のオプションのさまざまな値で,NIntegrateで使われるサンプル点の数を示す.
In[33]:=
Click for copyable input
Out[33]=

ガウス規則のクロンロッド拡張の実装の記述は[PiesBrand74]を参照のこと.

"GaussKronrodRule"のサンプル点と重み

次は,指定された数の粗い点と精度についてのガウス・クロンロッド横座標,重み,誤差の重みを計算する.
In[34]:=
Click for copyable input
Out[35]=

下の計算は,ガウス・クロンロッド積分則(を参照)の次数を示している.

ガウス・クロンロッド積分公式の次数を計算する.
In[36]:=
Click for copyable input
Out[36]=
関数を定義する.
In[37]:=
Click for copyable input

下のコマンドは積分推定値に対する積分則の重み付き総和 および誤差推定 を実装する.ここでは横座標,は重み,は誤差の重みである.

次は規則で計算した に対する積分および誤差推定である.
In[38]:=
Click for copyable input
Out[38]=
上の積分推定は厳密な結果に一致する.
In[39]:=
Click for copyable input
Out[39]=

上の誤差推定は,埋め込まれたガウス規則が次数の多項式に対して厳密であるため,ゼロにはならない.もしその次数の多項式を積分したら,誤差推定はゼロになる.

関数を定義する.
In[40]:=
Click for copyable input
その規則で計算した の積分と誤差推定である.
In[41]:=
Click for copyable input
Out[41]=
Integrateを使った厳密な結果.
In[42]:=
Click for copyable input
Out[42]=

"LobattoKronrodRule"

Lobatto積分則は事前に割り当てられた横座標を持つガウスタイプの規則である.これは積分区間の終点と区間内の最適サンプル点を使い,これらの点上の被積分関数値の重み付き総和を形成する.Lobatto規則のクロンロッド拡張は,Lobatto規則点の間に新しいサンプル点を加え,Lobatto規則の被積分関数の評価を再利用する高次の規則を形成する.

NIntegrateMethodオプションに値が与えられると,Lobatto規則のクロンロッド拡張を使う
In[43]:=
Click for copyable input
Out[43]=
オプション名
デフォルト値
"Points"5クロンロッド点で拡張されるガウス・Lobatto点の数
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

被積分関数 点のLobatto規則 は,次数の多項式について厳密である( が次数の多項式ならば である).

点のLobatto規則 のクロンロッド拡張個の点を加える.拡張された規則は が偶数ならば次数の多項式に対して, が奇数ならばの多項式に対して厳密である.Lobatto規則に関連付けられた重みはそのクロンロッド拡張で変化する.

と同様に,ガウス点の数はオプションで指定される.もし"Points"->n で使われると,規則の点の合計数はになる.

これはのオプションをさまざまな値にしてNIntegrateで使われるサンプル点の数を示したものである.
In[44]:=
Click for copyable input
Out[44]=

Lobatto規則は閉じた規則なので,被積分関数は区間の終点で評価されなければならない.これらの点に特異性がある場合,NIntegrateはそれを無視する.

Lobatto規則のクロンロッド拡張の実装の詳細は[PiesBrand74]を参照のこと.

"LobattoKronrodRule"のサンプル点と重み

次は指定された数の粗い点と精度に対するLobatto・クロンロッド横座標,重み,誤差の重みを計算する.
In[45]:=
Click for copyable input
Out[46]=

下の計算はLobatto・クロンロッド積分則の次数を示す(を参照).

Lobatto・クロンロッド積分則の次数を計算する.
In[47]:=
Click for copyable input
Out[47]=
関数を定義する.
In[48]:=
Click for copyable input

下のコマンドは積分推定に対する積分則の重み付き総和 と誤差推定 を実装する.ここで は横座標,は重み,は誤差の重みである.

この規則で計算された の積分と誤差推定である.
In[49]:=
Click for copyable input
Out[49]=
上の積分推定値は厳密な結果と一致する.
In[50]:=
Click for copyable input
Out[50]=

埋め込まれたLobatto規則は次数の多項式に対して厳密なので,上の誤差推定はゼロにはならない.もしその次数の多項式を積分したならば,誤差推定はゼロになる.

関数を定義する.
In[51]:=
Click for copyable input
その規則で計算した の積分と誤差推定である.
In[52]:=
Click for copyable input
Out[52]=
Integrateを使った厳密な結果.
In[53]:=
Click for copyable input
Out[53]=

"ClenshawCurtisRule"

クレンショウ・カーチス(ClenshawCurtis)則は被積分関数のチェビシェフ(Chebyshev)多項式近似から導かれたサンプル点を使用する.

NIntegrateではクレンショウ・カーチス則はMethodオプション値ので指定することができる.
In[54]:=
Click for copyable input
Out[54]=
オプション名
デフォルト値
"Points"5粗いクレンショウ・カーチスの点
"SymbolicProcessing"Automatic記号的処理を実行する秒数

のオプション

個のサンプル点を持つクレンショウ・カーチス則は理論的に次数 以下の多項式に対して厳密である.しかし,実際にはクレンショウ・カーチス則はガウス則の確度を達成する[Evans93][OHaraSmith68].クレンショウ・カーチス公式の誤差は[OHaraSmith68]で解析されている.

伝統的なクレンショウ・カーチス則のサンプル点は,チェビシェフ多項式の零点である.実践的なクレンショウ・カーチス則のサンプル点はチェビシェフ多項式の極値となるよう選ばれる.伝統的なクレンショウ・カーチス則は進行的でないが実践的なクレンショウ・カーチス則は進行的である[DavRab84][KrUeb98].

が,関数 個のサンプル点の実践的なクレンショウ・カーチス則を表すとする.

進行的属性は, のサンプル点がのサンプル点の部分集合であることを意味する.よって,差分は積分推定の誤差推定とすることができ,余分な被積分関数評価なしで計算することができる.

NIntegrateのオプションMethod->{"ClenshawCurtisRule","Points"->k}個の点を持つ実践的なクレンショウ・カーチス則を使う.
In[55]:=
Click for copyable input
Out[55]=
のオプションにさまざまな値を使い,NIntegrateにより使用されるサンプル点の数を示す.
In[56]:=
Click for copyable input
Out[56]=

"ClenshawCurtisRule"のサンプル点と重み

指定された数の粗い点と精度でのクレンショウ・カーチス則のサンプル点と重みである.
In[57]:=
Click for copyable input
Out[58]=
のサンプル点を計算する別の方法である.
In[59]:=
Click for copyable input
Out[60]=
関数を定義する.
In[61]:=
Click for copyable input
この規則で計算した の積分と誤差推定である.
In[62]:=
Click for copyable input
Out[62]=
Integrateによる厳密値である.
In[63]:=
Click for copyable input
Out[63]=

"MultipanelRule"

は2つ以上の隣接区間上の一次元積分則の適用を1つの規則に結合する.隣接区間のいずれかにもとの規則を適用することをパネルと呼ぶ.

を使った積分の例である.
In[7]:=
Click for copyable input
Out[7]=
オプション名
デフォルト値
Method"NewtonCotesRule"1つのパネルに対する横座標,重み,誤差の重みを提供する積分則指定
"Panels"5パネルの数
"SymbolicProcessing"Automatic記号的処理を実行する秒数

のオプション

単位区間を点 個の部分区間に分割するとする.

今,規則

があるとすると,区間 の規則

に変換される.ここで とすると, に基づく パネルの積分規則は明示的に

と書くことができる.もし が閉じている( がサンプル点としてを持つ)ならば, であり, のサンプル点の数は に還元される(これはの実装で実行される).

マルチパネル則(複合則とも呼ばれる)の理論は[KrUeb98]および[DavRab84]に記述されている.

"MultipanelRule"のサンプル点と重み

のサンプル点と重みは次のコマンドで得ることができる.
In[65]:=
Click for copyable input
Out[66]=
次はガウス・クロンロッド則の横座標と重みである.
In[67]:=
Click for copyable input
Out[67]=
マルチパネル則の横座標はRescaleを使って得ることができる.
In[68]:=
Click for copyable input
Out[68]=
次はもとの重みからマルチパネル則の重みを導出する方法を示している.
In[69]:=
Click for copyable input
Out[69]=

"CartesianRule"

次元の直積型積分則は 個の一次元則のサンプル点の直積であるサンプル積を持つ.直積型積分則のサンプル点に関連付けられた重みは,座標に対応する一次元則の重みの積である.

NIntegrateの直積の積分はMethodオプション値で指定することができる.
In[70]:=
Click for copyable input
Out[70]=
オプション名
デフォルト値
Method"GaussKronrodRule"直積型積分則が形成される規則あるいは規則のリスト
"SymbolicProcessing"Automatic記号的処理を実行する秒数

オプション

例えば,以下の公式があるとする.

これらはそれぞれ次数 の多項式に対して厳密である.よって 個の点を持つ公式

が次数の多項式に対して厳密であることは簡単に分かる.分点に関連付けられた重みは である.

個のサンプル点と重みを持つような, 個の一次元則に対する一般的な直積の公式は,以下の通りである.

明らかに(2)は

と書くことができる.ここで であり,各整数 について であり である.

次は直積型積分則による積分の可視化である. 軸に沿ってが, に沿ってが使ってある.
In[71]:=
Click for copyable input
Out[72]=

直積型積分則は比較的低い次元()に適用できる.それは高次元では,「組合せ的爆発」の対象となるからである.例えば,それぞれが個のサンプル点を持つつの同一の一次元則の次元直積は個のサンプル点を持つ.

NIntegrateでは,積分が多次元であり,Methodオプションに一次元則あるいは一次元則のリストが与えられている場合に直積型積分則が使われる.

で直積型積分則の積分を指定する.
In[73]:=
Click for copyable input
Out[73]=
一次元の積分則のリストで直積型積分則の積分を指定する.
In[74]:=
Click for copyable input
Out[74]=
一次元積分則のリストで直積型積分則の積分を指定する別の例.
In[75]:=
Click for copyable input
Out[75]=

直積型積分則についての詳細は[Stroud71]に記載されている.

"CartesianRule"のサンプル点と重み

則のサンプル点と重みはコマンドで求めることができる.
In[76]:=
Click for copyable input
Out[76]=

はそれぞれの規則の横座標と重みとを別々にしておく.そうでないと(3)から分かるように,高次元の場合に結果が大きすぎることがあるからである.

の結果は以下の関数を使い,(4)の形式にすることができる.
In[77]:=
Click for copyable input
In[78]:=
Click for copyable input
Out[78]=

"MultidimensionalRule"

立方体(ここで )に対する完全対称則は次の2つの属性を持つ点の集合からなる.まず1つ目は,集合内のすべての点は,その集合の任意の固定点の座標の置換および/または符号変更により生成できるというものである.2つ目は,その集合のすべての点は関連付けられた同じ重さの重みを持つということである.

NIntegrateでの完全対称多次元積分(完全対称立体求積法)はMethodオプション値ので指定することができる.
In[6]:=
Click for copyable input
Out[6]=

上の属性を満足する完全対称則の点集合は軌道と呼ばれる.不等式 が成り立つ座標に対する軌道上の点 はジェネレータと呼ばれる([KrUeb98][GenzMalik83]を参照のこと).

オプション名
デフォルト値
"Generators"5完全対称則のジェネレータの数
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

積分則が で表される 個の軌道を持つ場合,その 番目 はそれに関連付けられた重み を持ち,積分推定値は以下の式

で計算される.次数 のヌル規則は,次数 のすべての単項式でゼロに積分され,次数 の少なくとも1つの単項式ではそれに失敗する.それぞれのヌル規則は基礎的な積分則と低次数の適切な積分との間の差分と考えることができる.

NIntegrateオブジェクトは基本的に,1つの積分則と1つあるいはそれ以上のヌル規則とを合せた3つの異なる積分則オブジェクトのインターフェースである.そのジェネレータ数と次数は下の表にまとめてある.6個または9個のジェネレータを持つ規則オブジェクトは3つのヌル規則を使い,それぞれが2つのヌル規則の線形結合である.ヌル規則の線形結合は位相誤差を避けるために使われる.ヌル規則がどのように使われるかは[BerntEspGenz91]を参照のこと.

NIntegrateの完全対称則のジェネレータ数と次数.

ジェネレータ数積分則次数それぞれのヌル規則の次数参照
575[GenzMalik80]
675,3,1[GenzMalik83][BerntEspGenz91]
997,5,3[GenzMalik83][BerntEspGenz91]

(ここで )形式の積分に対する完全対称多次元積分則でNIntegrateにより使用されるサンプル点の数である.
In[80]:=
Click for copyable input
Out[81]=

"MultidimensionalRule"のサンプル点と重み

ここでは完全対称多次元則を使った積分推定値の計算例を示す.

以下はジェネレータ数に対するパラメータである.
In[82]:=
Click for copyable input
この関数はジェネレータ点を取りその軌道を生成する.
In[83]:=
Click for copyable input
指定されたジェネレータ数に対するジェネレータと重み.
In[84]:=
Click for copyable input
各ジェネレータの軌道を計算する.
In[89]:=
Click for copyable input
以下で関数を定義する.
In[90]:=
Click for copyable input
以下で多次元則を適用する.
In[92]:=
Click for copyable input
Out[92]//InputForm=
厳密な結果である.
In[93]:=
Click for copyable input
Out[93]//InputForm=
軌道の点に対するグラフィックスプリミティブを作成する.
In[94]:=
Click for copyable input
異なる軌道がどのように見えるかである.
In[95]:=
Click for copyable input
Out[95]=
すべての規則点を一緒にする.
In[96]:=
Click for copyable input
Out[96]=

"LevinRule"

レビン型規則は,不定積分を被積分関数の多項式と振動部の積として近似することにより振動関数の積分を推定する.

次はを使って計算された,高振動積分である.
In[45]:=
Click for copyable input
Out[45]=

はデフォルトで,自動的に被積分関数の振動部を検知し,下で説明する選点法を適用して積分推定を形成する.振動部を指定するためのオプション,および非振動の別の積分則へと適宜切り替えるかどうか等を含む数値法の細かい動作を制御するオプションを使うことができる.

オプション名
デフォルト値
"Points"Automatic粗い選点の数
"EndpointLimitTerms"Automatic非数値の終点における多項式補間の次数
"TimeConstraint"10積分ステップ毎に選点方程式を解く最大限の時間
"MethodSwitching"True非振動の規則へ自動的に切り替えるかどうか
"OscillationThreshold"10を適用し始める,最小の振動の概数
MethodAutomatic別の非振動規則
"LevinFunctions"Automaticカーネルに含ませる振動関数
"MaxOrder"50カーネルの最大微分階数
ExcludedForms{}カーネルから除外する形式
"AdditiveTerm"Automatic被積分関数の加算的な非振動部
"Amplitude"Automatic被積分関数でのカーネルの係数
"DifferentialMatrix"Automaticカーネルにより満足される線形常微分方程式の行列形式
"Kernel"Automatic明示的な振動カーネル

のオプション

レビン型選点法

スコープ

一次元の数値積分に対する基本的なレビン型のメソッド[Levin96]は次の形式の振動関数に適用される.

ここで振動部 は次の行列形式で書かれる階数 の線形常微分方程式を満足しなければならない.

この の一つである.関数 のベクトルは振動カーネルであり,行列 は微分行列である.

例えば,であるとすると,振動カーネルは次のように

となり,その微分行列は下のようになる.

レビン型のメソッドは次のようなより一般的なケース

を扱うこともできる.上のより単純な場合は,通常 および である.関数 のベクトルは振幅である.

ExpBesselJAiryAi等,数学における振動する特殊関数のほとんどは,線形常微分方程式を満足するため,レビン型メソッドを使って積分することができる.さらに,このような関数の積,総和,整数ベキ,高速でない振動関数との合成関数もすべて線形常微分方程式を満足する.DifferentialRootReduceも参照のこと.

法は次のような,わずかにより一般的なクラスにも適用することができる.

ここで加算項 等の従来の求積法で別に扱われる.

レビン型メソッドは,振幅 ,加算項 ,微分行列 自身が高速で振動しないときに効果的である[Levin97]

選点法

レビン型の選点法では,ある未知の関数に対して被積分関数の不定積分が形式で求められる.

について両辺を微分し,を使うと,すべての が打ち消され,関数 が非振動の微分方程式を満足することが分かる.

これらの方程式は,選点法を使って について解くことができる.選点法では,それぞれの 個の既定された基底関数 の単純集合の線形結合として近似する.

では,基底関数 はチェビシェフ(Chebyshev)の多項式ChebyshevT[k,x]となるよう選ばれる.これを の微分方程式に代入すると, 個の選点の係数 に対する 個の線形方程式ができる.

はすべて既知の関数なので,積分領域内の 個の選点ノードでこれらの線形関数を評価することができる.これにより, 個の係数 に対する 個の線形方程式ができる.これらの方程式を数値的に解くと, に対する近似形式が得られ,次のように微積分学の基本定理を使って直接積分を評価できるようになる.

NIntegrateでは,選点ノードとして対応する"GaussKronrodRule"の分点が使われる.

個の選点ノードを使って積分推定を計算する.ここで メソッドオプションの値である.階数の低い積分推定値もまた,2つおきの選点ノード()を使って計算される.これら2つの積分推定値の差分はアルゴリズムの誤差推定として取られる.

それぞれの積分部分領域にある23個の選点ノードを使って積分推定値を計算し,11個の選点ノードを使って誤差推定に対する階数の低い積分推定値を計算する.
In[88]:=
Click for copyable input
Out[88]=

振動カーネルの指定

デフォルトではは,できるだけ被積分関数を包含する振動カーネルを自動的に選ぶ.この動作はオプションを使って細部に渡り変更することができる.

被積分関数に,積分領域上であまり高振動ではない振動関数が含まれている場合,振動カーネルからそれを除外した方が効率よくなる可能性がある.例えば,Sin[x]は領域上ではあまり高振動ではない.

被積分関数の振動部を指定するためには,メソッドオプションを使うことができる.

カーネルSin[1000x]を使う.因数Sin[x]は振動カーネルには含まれていないが,振幅に含まれている.
In[157]:=
Click for copyable input
Out[157]=

メソッドは,自動的に検出される振動カーネルに含まれるべきではない被積分関数の部分を指定するために使うことができる.

次は,オプションを使って上と同じカーネルを指定する別の方法である.
In[158]:=
Click for copyable input
Out[158]=

メソッドオプションは,自動的に検出される振動カーネルに含まれる関数のタイプを指定するために使うことができる.設定は特定の関数 のみを検出する.

ExpSinにのみ使用する.BesselJは無視される.
In[28]:=
Click for copyable input
Out[28]=

設定は名前付きのグループ から関数を検出する.

三角関数関係および指数関数関係の関数を検出するときだけを使う.BesselJは振動カーネルには含まれない.
In[29]:=
Click for copyable input
Out[29]=

メソッドオプションで認識される名前付きグループ

デフォルトではの設定には以外の名前付きグループがすべて含まれる.この設定は,領域内のいずれかの部分で振動する初等関数および特殊関数に対応する.

NIntegrate`LevinIntegrandオブジェクト

NIntegrateにより検出された振動カーネル,差分行列,振幅,加算項を,内部関数を使って詳細に調べることができる.は,レビン型の選点法を適用することのできる完全に処理された被積分関数を表すオブジェクトを返す.

次は,処理されたが被積分関数sqrt(x)+x ⅇ^(ⅈ x) TemplateBox[{0, {100,  , x}}, BesselJ]を1個の振動カーネルとその他の部分に分解する.
In[4]:=
Click for copyable input
Out[4]=

オブジェクト li のプロパティは,形式 li["prop"] を使ってアクセスすることができる.

次は上の被積分関数に対する振動カーネル である.
In[6]:=
Click for copyable input
Out[6]=
これは同じ被積分関数に対する差分行列 である.多次元の場合,それぞれの次元の積分に対して異なる差分行列があるため,Firstが使われる.
In[9]:=
Click for copyable input
Out[9]//MatrixForm=
この振動カーネルに対する,レビンの微分方程式 を検証する方法である.
In[13]:=
Click for copyable input
Out[13]=
In[14]:=
Click for copyable input
Out[14]=
プロパティすべてのリストは li["Properties"]でアクセスできる.
In[15]:=
Click for copyable input
Out[15]=

プロパティは,レビン型の選点法で使われるの主なプロパティを与える.

上の被積分関数に対する振動カーネルと差分行列に加え,加算項 および振幅 である.
In[16]:=
Click for copyable input
Out[16]=

このプロパティを使うことで,振動カーネル検出に対するオプション等の異なるメソッドオプションの効果をすくに見ることができる.

次は同じ被積分関数に対する2種類の分解である.デフォルトではSin[x]Cos[x]の両方の項が含まれる.
In[19]:=
Click for copyable input
Out[19]=
"AdditiveTerm"->Sin[x]という設定では,Sin[x]の項は検出されるカーネルに含まれない.
In[20]:=
Click for copyable input
Out[20]=

被積分関数に関数のグループを含める,あるいはそれから除外するためにオプションを使う方法が分かる.

次はメソッドオプションの異なる設定の効果を見る方法である.Log[x]は振動関数ではないため,デフォルトでは検出されるカーネルから除外される.
In[43]:=
Click for copyable input
Out[43]=
"LevinFunctions"->Allという設定では,Log[x]は検出されるカーネルに含まれる.
In[44]:=
Click for copyable input
Out[44]=

被積分関数の分解における4つの部分(加算項,振幅,カーネル,差分行列)すべてを完全に指定することができる.これは同様の被積分関数でNIntegrateを複数回効率よく呼び出すときに便利である.

これは被積分関数の分解を完全に手動で指定する場合のの使い方である.
In[24]:=
Click for copyable input
Out[24]=

この場合,NIntegrateの第1引数が完全に無視されることに注意する.

次の積分は上記と等価である.NIntegrateの第1引数は無視される.
In[25]:=
Click for copyable input
Out[25]=

"MonteCarloRule"

モンテカルロ則はランダム(擬似ランダム)なサンプル点上で一様に重みの付いた被積分関数評価の総和を形成することにより,積分を推定する.

を1000個のサンプル点で使用する例である.
In[97]:=
Click for copyable input
Out[97]=
オプション名
デフォルト値
"Points"100サンプル点の数
"PointGenerator"Randomサンプル点座標ジェネレータ
"AxisSelector"Automatic大域的適応的モンテカルロ積分が使用されたときに軸を分割する選択アルゴリズム
"SymbolicProcessing"Automatic記号的前処理を実行する秒数

オプション

モンテカルロ法[KrUeb98]では, 次元積分 は次の期待(平均)値と解釈される.

ここで は, 上の一様分布,つまり確率密度が の分布について,乱数変数として解釈される関数 の平均値である.領域 の特性関数は と表され, の体積は と表される.

期待値 の単純モンテカルロ推定値は, 個の独立した乱数ベクトル を密度 で取る(つまりベクトルは 上で一様に分布するということ)ことで得られ,以下の推定値が求められる.

注意.関数 および 全体で非負であるため,有効な確率密度関数である.

大数の法則に従うと,収束

は確率で起る.大数の法則は誤差 に関する情報は提供しないので,確率的推定が使われる.

と定義されるとする.式(5)は の不偏推定値(のさまざまな集合に対する の期待値が ということ)であり,その分散は

である.ここで の分散を表し, の標準誤差はとなる.

実際には は既知でないため,以下の公式

で推定される.すると, の標準誤差は

となる.モンテカルロ推定値の結果は と書くことができる.

方程式(6)から分かるが,単純モンテカルロ推定値の収束率は積分の次元 に依存しておらず, 個のサンプル点が使われるならば収束率はになる.

NIntegrateの積分則は推定値 を計算する.

推定値は付加的に向上させることができる.つまり,もし推定値がであり,新しい付加的なサンプル関数値の集合があるならば,(7)と(8)を使い,下が求められる.

推定値を求めるためには,推定値を計算するために使われた乱数点が分かっている必要はない.

"AxisSelector"

が多次元の大域的適応型積分に対して使われるとき,それは2つの方法で適用される積分部分領域の分割軸を選ぶ.1つ目の選び方は無作為な選び方である.2つ目は,部分領域が軸に沿って分割される場合に,その部分領域のそれぞれの積分推定値の分散の和が最小になるようにする方法である.

無作為な軸の選択は次のように行われる.が選択に使用する軸の集合 を保持する.初めは にすべての軸が含まれる. の要素が無作為に選択される.選択された軸は から除かれる.次の積分推定の後,軸が から選ばれそこから除外される.これが継続していく. が空になると,すべての軸で補充される.

分散最小化の軸選択は次のように行われる.領域上での積分中,サンプル点の部分集合とその被積分関数値が保存される.次に,すべての軸に対して,軸に沿った分割で生成される2つの部分領域の分散が,保存されているサンプル点と対応する被積分関数値を使って推定される.これらの分散の総和が最小となる軸が分割軸として選ばれる.領域がその軸で分割されると新しい積分誤差推定が最小となるということを意味するからである.ある軸について,すべての保存されている点が半領域のどちらかでクラスタ化されると,その軸が分割に選ばれる.

オプション値
Random無作為な分割軸の選択
MinVariance|{MinVariance, "SubsampleFraction"->frac}新しい領域の分散和を最小にするような分割軸の選択

のオプション値

オプション名
デフォルト値
"SubsampleFraction"1/10分割軸を決定するために使われるサンプル点の分数

のオプション

のオプションを使った例である.
In[98]:=
Click for copyable input
Out[98]=

下の例では,2つの軸選択アルゴリズムを比較する.一般に分散最小化の選択法の方が使用するサンプル点の数が少ない.しかしこの方法を使うとが遅くなる.よって,どちらの軸選択方法でも同数のサンプル点を使う結果になるような積分では,無作為選択法を使った方が速い.また,分散最小化の分割軸を選ぶために多数のサンプル点を使うと積分が遅くなる.

以下の関数を考える.
In[2]:=
Click for copyable input
Out[3]=
次は無作為の軸選択を使った,上の関数に対する適応型モンテカルロ積分のサンプル点である.
In[4]:=
Click for copyable input
Out[5]=
下は分散を最小化する分割軸を選んだときのサンプル点である.上のモンテカルロ積分と比較すると,こちらのサンプル点の方が円 の周囲に密集しており,その数はほぼ半分である.
In[6]:=
Click for copyable input
Out[6]=
無作為軸選択を使用する適応型モンテカルロ積分.
In[104]:=
Click for copyable input
Out[104]=
分散最小化の軸選択を使用する適応型モンテカルロ積分である.これはランダム選択を使うよりも遅い.
In[105]:=
Click for copyable input
Out[105]=
分散最小化の軸選択で多数の保存された点を使うと,積分が遅くなる.
In[106]:=
Click for copyable input
Out[106]=

規則の比較

以外の積分則はすべてNIntegrate適応型ストラテジーで使われる.積分ストラテジーに対して積分則の要素のタイプや点の数を変更することで,異なる積分アルゴリズムとなる.一般に,これらの異なる積分アルゴリズムは,異なる積分に対して実行方法が異なる.そこで次のような問題が生じる.

1.  すべての積分,あるいはあるタイプの積分に対して,他よりもよいという規則のタイプがあるか.

2.  積分ストラテジーが与えられると,どの規則がそれとうまく動作するか.また,どの積分に対してか.

3.  積分,積分ストラテジー,積分則が与えられた場合,規則の中の何番の点が,目標精度を満足する積分推定値に到達するために使用されるサンプル点の合計数を最小化するか.

指定された積分,積分ストラテジーで,サンプル点の最少数で目標精度を満足する結果に達する積分則のことを最適積分則という.最適積分則を決定するいくつかの要因を挙げる.

4.  一般に,規則の次数が高いほど,滑らかな被積分関数および高い目標精度での積分は速くなる.一方,被積分関数に対して規則の次数が高すぎると,適応ストラテジーを使うときに使用されるサンプル点が多すぎる.

5.  規則の誤差推定汎関数は,積分ストラテジーによる仕事の総量に大きく影響を及ぼす.少数の点を持つ規則は,(1)積分の過小評価のために誤った結果となるか (2) 積分の過大評価のために多すぎるサンプル点を適用することになるかのどちらかである(「異常な挙動の例」を参照のこと).さらに,誤差推定汎関数は1つあるいはいくつかの埋め込まれたヌル規則で計算されることがある.一般に,ヌル規則の数が多いほど,誤差推定はよくなる.位相誤差が小さくなるからである.ヌル規則の数と誤差推定を計算する総和においてそのヌル規則に割り当てられた重みは,異常な積分およびその規則では計算しにくい積分の集合を決定する(NIntegrateの多次元規則の中には,誤差推定を計算するためにいくつかの組み込まれたヌル規則を使うものもある.NIntegrateの一次元積分則はすべてヌル規則を1つだけ使う).

6.  局所的適応型ストラテジーは,開規則(等)やサンプル点が非一様に分布する閉規則(等)よりも,サンプル点がより一様に分布する閉規則(等)でより効率的である.

7.  ストラテジーにより再使用される点の割合により何が最適の規則かが決まる.一次元積分では,は閉規則のすべての点を再使用する.は誤差推定値の向上が必要である領域のほぼすべての点を捨てる.

規則の中の点の数

ここでは,被積分関数が滑らかで目標精度が高い場合,規則の次数が高いほど積分が速いという例題を挙げる.また,被積分関数に対して規則の次数が高すぎるため,適応ストラテジーが被積分関数の非連続性を操作する場合に使うサンプル点が多すぎるという例も示す.すべての例題にはガウス則のBerntsenEspelid誤差推定が使われている.

区間のガウス則の誤差は

である([DavRab84]を参照).

Integrateの厳密値を使って,積分 の規則の誤差を計算する関数である.
In[107]:=
Click for copyable input
関数のリストを定義する.
In[108]:=
Click for copyable input
区間で関数をプロットする.
In[109]:=
Click for copyable input
Out[109]=
点の範囲について に対するの誤差を計算する.
In[110]:=
Click for copyable input
それぞれの関数に対して誤差の対数がどのように減少するかのプロットである.点の数が増えると,非連続関数および非連続導関数を持つ関数の積分推定値が向上することが分かる.
In[111]:=
Click for copyable input
Out[117]=

サンプル点の最少数

以下は積分で使用されるサンプル点の数を見付ける関数である.
In[118]:=
Click for copyable input
次は目標精度の範囲および積分則の粗い点の範囲に対して使用されるサンプル点の数を見付ける.
In[19]:=
Click for copyable input
以下は各精度に対するサンプル点の最少合計数を見付ける.これで使用される粗い積分則の点の数も見付かる.
In[121]:=
Click for copyable input
以下は目標精度,および最少数の合計サンプル点を使う積分則の点の数のプロットである.
In[122]:=
Click for copyable input
Out[128]=

規則の比較

Integrateによる厳密値を使って積分 に対する規則の誤差を計算する関数である.
In[129]:=
Click for copyable input
以下で関数のリストを定義する.
In[130]:=
Click for copyable input
区間の関数のプロットである.
In[131]:=
Click for copyable input
Out[131]=
積分 のそれぞれに対しての誤差を計算する.
In[132]:=
Click for copyable input
各規則と関数に対して誤差の対数がどのように減少するかのプロットである.
In[136]:=
Click for copyable input
Out[136]=

異常な挙動の例

誤差推定器をごまかす

ここではNIntegrateが被積分関数の部分を検出することができないため,デフォルト設定で過小評価する積分について説明する.その部分は目標精度を上げると検出される.

誤った推定値

以下の関数を考える.
In[13]:=
Click for copyable input
上での厳密な積分は以下のようになる.
In[138]:=
Click for copyable input
Out[138]=
NIntegrateは,推定値を与える.
In[139]:=
Click for copyable input
Out[139]=
これは厳密値と比較すると不正確である.
In[140]:=
Click for copyable input
Out[140]=
以下は関数のプロットであるが,これも誤りである.
In[141]:=
Click for copyable input
Out[141]=

よりよい結果

NIntegrateのオプションPrecisionGoalを使い,反復の深さを上げるとよりよい結果が得られる.
In[17]:=
Click for copyable input
Out[17]=
次はよい結果が計算されない目標精度を見付ける表である.
In[18]:=
Click for copyable input
Out[18]=
プロット点が増加すると,関数のプロットは異なって見える.
In[144]:=
Click for copyable input
Out[144]=
Plotのデフォルトオプションでは表示されない山形のプロットをズームする.
In[145]:=
Click for copyable input
Out[145]=
関数のこの部分が積分されると,結果はNIntegrateのデフォルトオプション設定では失われている数量にフィットする.
In[146]:=
Click for copyable input
Out[146]=
In[147]:=
Click for copyable input
Out[147]=

推定器が誤解される理由

NIntegrateのデフォルトで使用されるガウス・クロンロッド則の横座標と重みである.
In[146]:=
Click for copyable input
以下は上の規則の適用に対する関数を定義する.
In[147]:=
Click for copyable input
適応型ストラテジーが被積分関数をサンプルする点を見付ける.
In[148]:=
Click for copyable input
以下はサンプル点のプロットである.縦軸は被積分関数を評価するために使用された点の次数である.
In[149]:=
Click for copyable input
Out[149]=

上のプロットではNIntegrate付近の第2の山形の頂上付近で活発な計算を実行しているのが分かる.NIntegrate付近の積分されない山形付近ではあまり計算を行っていない.

領域上のサンプル点の最後の集合のガウス・クロンロッドおよびガウスの横座標である.
In[150]:=
Click for copyable input
Out[150]=
Out[151]=
被積分関数が横座標上に適用される.
In[152]:=
Click for copyable input
以下は横座標上の被積分関数の多項式近似である.
In[154]:=
Click for copyable input
これらの関数は2つの多項式近似がほぼ 上で一致することを示している.
In[156]:=
Click for copyable input
Out[156]=
Out[158]=
多項式がの位置する領域上で積分されると,NIntegrateが誤差推定として使用するその多項式間の差分はかなり小さい.
In[159]:=
Click for copyable input
Out[159]=
Out[160]=
Out[161]//FullForm=

上記の誤差は領域に対して割り当てられた誤差推定なので,デフォルトの目標精度のNIntegrateではそれ以上の積分調整のためにそれを選ぶことはない.

位相誤差

ここでは積分則が,積分推定の実際の誤差を過小評価・過大評価する理由を説明する.同様の説明が[LynKag76]でもされている.

関数を定義する.
In[162]:=
Click for copyable input
領域上でのの積分の数値的および記号的評価を考える.
In[163]:=
Click for copyable input
Out[163]=
In[164]:=
Click for copyable input
Out[164]=
これらは非常に異なる.要求される目標精度は2であるが相対誤差はよりもずっと高い.
In[165]:=
Click for copyable input
Out[165]=

NIntegrateはより高い目標精度に対して正しい結果を与える.)

なぜこれが起るのであろうか.

積分則 が規則 に埋め込まれているとする.偶然積分推定 (ここで )の誤差推定が厳密な誤差と比べて小さすぎることがある.

これを例示するために,5個のサンプル点のガウス則 が埋め込まれた,サンプル点11個のガウスクロンロッド則を考える(これは上の2つの積分で使った規則である).
In[166]:=
Click for copyable input
以下はこの規則を適用する関数である.
In[167]:=
Click for copyable input
次は上で定義されたの積分 である.
In[168]:=
Click for copyable input
における の異なる値に対する推定誤差および実際の誤差を使ってグラフをプロットすることができる.つまり,をプロットするのである.
In[169]:=
Click for copyable input
Out[169]=

上のプロットで,青は推定誤差を表している.実際の誤差は赤で表してある.

パラメータ の値の極小値の1つに非常に近い.

一次元積分則は,規則の横座標およびその被積分関数値にフィットする多項式の積分の結果と見ることができる.さらにの積分に対して実際にフィットする多項式を見ようとすることもできる.

In[170]:=
Click for copyable input

下のプロットでは,関数が赤で,ガウス多項式が青で,ガウス・クロンロッド多項式が紫で表示されており,ガウスサンプル点が黒い点,ガウス・クロンロッドサンプル点が赤い点になっている.

の頂上が2つの横座標のほぼ半分のところにあるため,その近似は過小評価される.

In[172]:=
Click for copyable input
Out[172]=

逆に,の頂上は横座標のほぼどちらかとなっているため,近似は過大評価となる.

In[173]:=
Click for copyable input
Out[173]=

専門用語の索引

横座標

一次元積分則の次数

多次元積分則の次数

厳密規則

埋込み規則

ヌル規則

直積型積分則

進行的規則

サンプル点