9.3 数値積分 (3) 実用プログラムのベンチマーク

9.3.1 実用プログラムの選択

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

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

XLPack に収録されているプログラム (有限区間: [a, b])

※それぞれのリバースコミュニケーション版も使用可能(プログラム名に”_r”のサフィックスが付く)

XLPack に収録されているプログラム (有限区間: [a, b], 特殊な被積分関数用)

※それぞれのリバースコミュニケーション版も使用可能(プログラム名に”_r”のサフィックスが付く)

XLPack に収録されているプログラム (半無限区間／無限区間: [a, \(\infty\)]／[\(-\infty\), \(+\infty\)])

※それぞれのリバースコミュニケーション版も使用可能(プログラム名に”_r”のサフィックスが付く)

9.3.2 有限区間積分

それぞれの例題において, 種々の実用プログラムで計算したときの被積分関数評価回数を横軸に, その時の相対誤差を縦軸に対数プロットする.

ここで比較する実用プログラムは (適応) 自動積分プログラムであるから, 要求精度を適当に変化させて計算を行い, そのときの実際の関数評価回数と相対誤差をプロットした.

全体的には汎用の Qag がよい成績を示した. Qag が使用する固定分点プログラムを Qk15 (15 点公式), Qk31 (31 点公式), Qk61 (61 点公式) と 3 種類変えて比較したが, どれがよいかは被積分関数によって異なった. 端点に特異点がある場合には Defin が有効なことが多い. Qawo, Qagp, Qaws などの特定の関数向けのプログラムは, それが使える関数であれば使った方がよい. なめらかな関数であれば, 普通のガウス公式や Qng, Qk15 等の区間分割を行わないプログラムがよいこともあるが, これらは全くだめなこともあるので検算しながら使った方がよい.

数値実験 (1) 解きやすい問題

はじめに比較的解きやすい問題を使って有限区間積分用プログラムの傾向を見ておく. 例題は「9.2.4.1 適応自動積分プログラム QAG」の例題として使用した次の関数である.
\[
f(x) = \frac{1}{(x – 0.3)^2 + 0.01} + \frac{1}{(x – 0.9)^2 + 0.04} – 6
\]

これを Qng, Qag (Qk15, Qk31, Qk61), Defin で計算した結果は次のとおりである.

比較のために 2, 3, …, 31 点ガウス公式単体の結果も表示した. Qng は 21, 43 または 87 点ガウス・クロンロッド公式を自動選択するプログラムである.

ガウス公式と Qng の結果から, この問題は 80 点程度のガウス公式単体を使えば最高精度で計算できそうだと推定できる. 他の自動積分プログラムの結果は滑らかな関数に対する基本性能の傾向を表している.

数値実験 (2) 端点に特異点

数値実験 (2) 以下では, 実用プログラムの比較として比較的解きにくい例題のみを扱う. それぞれの例題の詳細は文献[6]を参照のこと.

次の積分を求める.
\[
\int_0^1 x^c log(1/x) dx = (c + 1)^{-2}
\] 端点 (x = 0) に特異点がある.

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

これを Qng, Qag, Defin で計算した結果は次のとおりである.

比較のためにガウス公式単体 (N = 2 ～ 31) の結果も示した. 端点に特異点があるので, Defin が他の公式よりも速い.

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

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

数値実験 (3) 端点に特異点 (2)

次の積分を求める.
\[
\int_0^{\pi/2} sin(x)^{c-1} dx = \int_0^{\pi/2} x^{c-1}(sin(x)/x)^{c-1} dx = 2^{c-2}\Gamma(\frac{c}{2})^2/\Gamma(c)
\] 端点 (x = 0) に特異点がある.

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

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

数値実験 (4) 鋭いピーク

次の積分を求める.
\[
\int_0^1 \frac{4^{-c}}{(x – \frac{\pi}{4})^2 + 16^{-c}} dx = arctan((4 – \pi)4^{c-1} + arctan(\pi 4^{c-1})
\] c = 1 のときの結果は次のとおりである.

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

C = 15 とすると鋭いピーク (x = \(\pi\)/4) が現れてくる.

この場合には小区間に分割するプログラムでないと計算できなくなった. Qag 以外では解が求められなかった.

数値実験 (5) 区間内に特異点

次の積分を求める.
\[
\int_0^1 |x – \frac{1}{3}|^c dx = \frac{(2/3)^{c+1} + (1/3)^{c+1}}{c + 1}
\] 区間内 (c = 1/3) に特異点がある.

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

この場合は Qag でも解を求めることができたが, 特異点が既知であれば Qagp で特異点を指定するのが有効である.

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

特異点のところでなめらかに x 軸に接するようになったためか, 計算できるプログラムが増えた. ただし, Qagp が有効なことは変わらない.

数値実験 (6) 振動型積分

次の積分を求める.
\[
\int_0^1 e^{20(x – 1)} sin(2^c x) dx = \frac{20sin(2^c) – 2^c cos(2^c) + 2^c exp(-20)}{400 + 4^c}
\] 振動型積分 (重み関数 \(w(x) = cos(\omega x)\) または \(sin(\omega x)\) の有限区間の積分) である. 専用プログラム Qawo を適用することができる.

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

このような場合は Qawo を使って sin の項を分離するとよい. この場合は通常の (小区間に分割しない) ガウス公式でも十分計算できている.

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

今度は通常のガウス公式ではうまく計算できない.

数値実験 (7) 端点に特異点がある振動型積分

次の積分を求める.
\[
\int_0^1 (x(1 – x))^{-1/2} cos(2^c x) dx = \pi cos(2^{c-1}) J_0(2^{c-1})
\] 振動型積分であるが端点に特異点がある.

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

Qawo 以外では \(10^{7}\) 程度の精度までしか求められなかった.

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

状況は c = 8 のときと同様である. Defin は速いが \(10^{-8}\) 付近で精度が頭打ちになった.

9.3.3 無限・半無限区間積分

それぞれの例題において, 種々の実用プログラムで計算したときの被積分関数評価回数を横軸に, その時の相対誤差を縦軸に対数プロットする.

ここで比較する実用プログラムは (適応) 自動積分プログラムであるから, 要求精度を適当に変化させて計算を行い, そのときの実際の関数評価回数と相対誤差をプロットした.

Qagi は, 半無限積分を有限区間 [0, 1] の積分に変数変換しそれに 15 点ガウス・クロンロッド公式を適用することにより必要な積分を求める.
\[
\int_a^{\infty} f(x) dx = \int_0^1 \frac{f(a + (1 – t)/t)}{t^2} dt
\] また, Qagi は無限区間の積分にも適用できるが, 無限積分は 2 つの半無限積分の和として計算される.
\[
\int_{-\infty}^{+\infty} f(x) dx = \int_0^{\infty} (f(x) + f(-x)) dx
\] これに対して, DE 公式 Dehint と Deiint は本文で説明した変換関数を使う (半) 無限区間専用のプログラムである.

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

数値実験 (8) 指数積分 E1(a)

次の半無限積分を求める.
\[
\int_a^{\infty} \frac{e^{-x}}{x} dx = E1(a)
\] E1(a) は特殊関数としてよく出てくる指数積分である. Qagi と Dehint を使った計算結果を示す.

比較のためにガウス・ラゲール公式 (2 ～ 31点) の結果も示した.

数値実験 (9) ガウス積分

次の無限積分を求める.
\[
\int_{-\infty}^{+\infty} x^4 e^{-x^2} dx = \frac{3}{4} \sqrt{\pi}
\] ガウス積分と呼ばれ, 積分に \(x^4\) を含む.

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

数値実験 (10) 非常に早く 0 に収れんする関数

次の半無限積分を求める.
\[
\int_0^{\infty} x^2 exp(-2^{-c}x) dx = 2^{3c+1}
\] 非常に早く 0 に収れんする関数で有限区間用積分プログラムでも計算できる. c = 3 としたとき次のようになった.

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

数値実験 (11) ゆっくりと 0 に収れんする関数

次の半無限積分を求める.
\[
\int_0^{\infty} \frac{x^{c-1}}{(1 + 10x)^2} dx = 10^{-c}(1 – c)\pi/sin(\pi c)
\] ゆっくりと 0 に収れんする関数で, こうなると有限区間用積分プログラムでは計算できない. c = 1.5 としたとき次のようになった.

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

数値実験 (12) フーリエ積分

次の半無限積分 (フーリエ積分) を求める.
\[
\int _0^{\infty} x^{-1/2} exp(-2^{-c}x)cos(x) dx = \sqrt{\pi} (1 + 4^{-c})^{-1/4} cos(\frac{arctan(2^c)}{2})
\] 比較的速く 0 に近づいていく関数である. c = 2 の時の結果は次のとおりである.

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

数値実験 (13) フーリエ積分

次の半無限積分(フーリエ積分)を求める.
\[
\int_0^{\infty} log(x)sin(x) dx = -\gamma
\] 振幅が徐々に大きくなる関数である.

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