9. 数値積分 (3) 実用ルーチンのベンチマーク

実用ルーチンの選択

XLPackに収録されている数値積分用のプログラムは, QUADPACK(文献[1])のプログラムと二重指数関数(DE)公式のプログラムに分けられる. 前者はSLATECに収録されたQUADPACKのプログラムを元にしており, 基本的には通常の関数はガウス・クロンロッド公式, 特殊な関数はクレンショウ・カーチス公式を使用している. 後者は主に文献[2]および[3]を元に作成したもので, 二重指数関数(DE)公式を使用している.

表において「適応自動積分」は要求精度を満足するように積分区間を適応的に自動分割して計算するものを指す. 適応的とは, 関数の動きをみて場所により分割を細かくしたり大きくしたり調節することをいう.　「自動積分」は, 要求精度を満足するように分点数(きざみ幅)を(場所によらず一律に)自動選択するものを指す.

XLPackに収録されているルーチン (有限区間: [a, b])

XLPackに収録されているルーチン (有限区間: [a, b], 特殊な被積分関数用)

XLPackに収録されているルーチン (半無限区間／無限区間: [a, +∞]／[-∞, +∞])

有限区間積分

それぞれの例題において, 種々の実用ルーチンで計算したときの被積分関数評価回数を横軸に, その時の相対誤差を縦軸に対数プロットする.
ここで比較する実用ルーチンは(適応)自動積分ルーチンなので, 要求精度を適当に変化させて計算を行い, そのときの実際の関数評価回数と相対誤差をプロットした.

全体的には汎用の Qag がよい成績を示した. Key = 1 (15点公式), 3 (31点公式), 6 (61点公式) と変えて比較したが, どれがよいかは被積分関数によって異なった. 端点に特異点がある場合には Defin が有効な可能性が高い. Qawo, Qagp, Qaws などの特定の関数向けのルーチンは, それが使える関数であれば使った方がよさそうである. なめらかな関数であれば, 普通のガウス公式や Qng, Qk15等の区間分割を行わないルーチンがよいこともあるが, これらは全くだめなこともあるので検算しながら使った方がよさそうである.

数値実験 (0)

はじめに比較的解きやすい問題を使って有限区間積分用ルーチンの傾向を見ておく. 例題は適応自動積分の説明に使用した次の関数である.

  
  f(x) = 1/((x - 0.3)2 + 0.01) + 1/((x - 0.9)2 + 0.04) - 6

これを Qng, Qag (Key = 1, 3, 6), Nc8, Defin で計算した結果は次のとおりである.

ここで, Nc8 は8次のニュートン・コーツ公式を使った適応自動積分ルーチンである(詳細は文献[2]を参照). ニュートン・コーツ公式は分点を等間隔に取れるため適応自動積分ルーチンを作るのに向いており, 以前はよく使われた. Nc8 はXLPackには含まれていないが, ガウス・クロンロッド公式以外の適応自動積分ルーチンと比較するために表示した. また, 同じく比較のために31次までのガウス公式単体の結果も表示した.

ガウス公式とQngの結果から, この問題は80次程度のガウス公式単体を使えば最高精度で計算できそうである. 他の自動積分ルーチンの結果は滑らかな関数に対するその基本的な性能の傾向を表しているものと思われる.

数値実験 (1)

数値実験(1)以下では, 実用ルーチンの比較なので比較的解きにくい例題のみを扱う. それぞれの例題の詳細は文献[1]を参照のこと.

次の積分を求める.

  
  ∫xclog(1/x)dx [0, 1] = (c + 1)-2

 
端点 (x = 0) に特異点がある.

c = 2.6 のときの関数の形は次のとおりである.

これを Qng, Qag (Key = 1, 3, 6), Qags, Defin で計算した結果は次のとおりである.

比較のために通常の(自動積分でない)ガウス公式(N = 2～31)の結果も示した. この関数は小区間に分割しなくても高精度のガウス公式で計算できるようである. 端点に特異点があるので, Defin が他の公式よりも速い.

同じ関数で c = 0 とすると次のようになった.

この場合は逆に小区間に分割しないとうまく求めることができないようである. Qag も Key = 1 (15分点)の場合が速い. Defin はわすか51回の関数計算で 10^-16 のオーダーの精度に達した.

数値実験 (2)

次の積分を求める.

  
  ∫4-c/((x - π/4)2 + 16-c)dx [0, 1] = arctan((4 - π)4c-1 + arctan(π4c-1)

 
c = 1 のときの結果は次のとおりである.

これも小区間に分割しなくても高精度のガウス公式で計算できるようで, Qng は43分点の公式を使用して 10^-16 のオーダーの精度で解を求めた.

C = 15 とするとピーク(x = π/4)が尖ってくる.

こうなると小区間に分割するルーチンでないと計算できなくなった. Qag 以外では解が求められなかった.

数値実験 (3)

次の積分を求める.

  
  ∫|x - 1/3|cdx [0, 1] = ((2/3)c+1 + (1/3)c+1)/(c + 1)

 
区間内(c = 1/3)に特異点がある.

c = 0.5 のときの結果は次のとおりである.

Qag でも解を求めることができるが, こういう場合には Qagp で特異点を指定してやるのが有効である. Defin は区間内に特異点がある場合には効果的ではない.

c = 2.5 とすると次のようになった.

特異点のところでなめらかにx軸に接するような形になっているためかどの公式でも計算できるようになった. ただし, Qagp はあいかわらず有効なようである.

数値実験 (4)

次の積分を求める.

  
  ∫sin(x)^(c-1)dx [0, π/2] = ∫xc-1(sin(x)/x)^(c-1)dx [0, π/2] = 2c-2(Γ(c/2)2)/Γ(c)

 
端点 (x = 0) に特異点がある.

c = 1.5 のときの結果は次のとおりである.

被積分関数を変形して上の2番目の形で Qaws を適用すると効果的であった. なお, 端点に特異点があるため Defin が非常に効果的である.

数値実験 (5)

次の積分を求める.

  
  ∫exp(20(x - 1))sin((2c)x)dx [0, 1] = (20sin(2c) - 2ccos(2c) + 2cexp(-20))/(400 + 4c)

 
振動型積分である.

c = 6 のときは次のようになった.

このような場合は Qawo を使ってsin()の項を分離するとよい. この場合は普通の(小区間に分割しない)ガウス公式でも十分計算できている. Qag (Key = 6) も1回で計算が終了している.

c = 8 とすると様子が全く異なり次のようになった.

今度は普通のガウス公式ではうまくいかない. しかし, Qag (Key = 6) はよい成績を示している.

数値実験 (6)

次の積分を求める.

  
  ∫(x(1 - x))-1/2cos((2^c)x)dx [0, 1] = π cos(2c-1)J0(2c-1)

 
これも振動型積分であるが端点に特異点がある. c = 8 のときは次のようになった.

Qawo 以外では 10^-7 くらいの精度までしか求められなかった.

c = 1 としたときは次のようになった.

状況は似ている. Defin は速いが 10^-8 付近で精度が頭打ちになった. Qag も同様である. この積分は cos() の項を分離してやらなければ精度が出ないのかもしれない.

無限・半無限区間積分

それぞれの例題において, 種々の実用ルーチンで計算したときの被積分関数評価回数を横軸に, その時の相対誤差を縦軸に対数プロットする.
ここで比較する実用ルーチンは(適応)自動積分ルーチンなので, 要求精度を適当に変化させて計算を行い, そのときの実際の関数評価回数と相対誤差をプロットした.

Qagi は, 半無限積分を有限区間[0, 1]の積分に変数変換しそれに15点ガウス・クロンロッド公式を適用することにより必要な積分を求める.

  
  ∫ f(x)dx [a, +∞] = ∫ f(a + (1 - t)/t) / t^2 dt [0, 1]

 
また, Qagi は無限区間の積分にも適用できるが, 無限積分は2つの半無限積分の和として計算する.

  
  ∫ f(x)dx [-∞, +∞] = ∫ (f(x) + f(-x)) dx [0, +∞]

 
これに対してDE公式を使う Dehint と Deiint は本文で説明した変換関数を使う(半)無限区間専用のルーチンである.

フーリエ積分用の Qawf は, 積分区間を小区間に分割しそれを順次 Qawo で積分する. 収束を速くするために加速法を使う. Deoint は本文で説明した変換関数のDE公式を使うフーリエ積分用のルーチンである.

数値実験 (1)

次の半無限積分を求める.

  
  ∫ exp(-x)/x dx [a, +∞] = E1(a)

 
E1(a) は特殊関数としてよく出てくる指数積分である.

比較のためにガウス・ラゲール公式(2～31点)の結果も示した. Dehint は 10^-10 の精度で頭打ちになった. Qagi は計算量は多いが 10^-16 の精度で計算できた.

数値実験 (2)

次の無限積分を求める.

  
  ∫ x4exp(-x2) dx [-∞, +∞] = 3/4 √π

 
ガウス積分と呼ばれる積分に x⁴ を含む積分である.

比較のためにガウス・エルミート公式(2～31点)の結果も示した. この積分は指数の項を除くと4次の多項式なので3次以上のガウス・エルミート公式は(丸め誤差を除いて)正しい値を与える. そのため図のように1点を除いて 10^-15 の精度付近にプロットされている. Qagi は半分づつに分けて計算するので, 手間はおおむね2倍になっていると思われその分計算量的には不利である.

数値実験 (3)

次の半無限積分を求める.

  
  ∫ x2exp(-2-cx) dx [0, +∞] = 23c+1

 
非常に早く0に収れんする関数で有限区間用積分ルーチンでも計算できる. c = 3 としたとき次のようになった.

ガウス公式は積分区間を[0, 300]とした. それでも 10^-13 以上の精度で計算できた.

数値実験 (4)

次の半無限積分を求める.

  
  ∫ xc-1/(1 + 10x)2 dx [0, +∞] = 10-c(1 - c)π/sin(πc)

 
ゆっくりと0に収れんする関数で, こうなると有限区間用積分ルーチンでは計算できない. c = 1.5 としたとき次のようになった.

端点(x = 0)が特異点なので有限区間の場合と同じようにDE公式の Dehint が有効である.

数値実験 (5)

次の半無限積分(フーリエ積分)を求める.

  
  ∫ x-1/2exp(-2-cx)cos(x) dx [0, +∞] = √π　(1 + 4-c)-1/4cos(arctan(2c)/2)

 
比較的速くゼロに近づいていく関数である. c = 2 の時の結果は次のとおりである.

Deoint では2種類の変換関数を使うことができる. 数値積分(2)で説明した文献[7]の関数を Deoint, 文献[6]の関数を Deoint2 と表示した. 後から提案された Deoint の方が少し速いことがわかる.

数値実験 (6)

次の半無限積分(フーリエ積分)を求める.

  
  ∫ log(x)sin(x) dx [0, +∞] = -γ

 
振幅が徐々に大きくなる関数である.

この積分は Qawf ではエラーになり計算できなかった.