Integral1 - c++で自己流数値計算

積分ルーチン

はじめに

信頼性と存在意義

ここでは、LibMK4に含まれる積分ルーチンを紹介する。
ただし、信頼性や計算速度といった観点から、できることなら以下で紹介するものは使わないことを強くお勧めする。c言語から呼び出せる積分ルーチンはGSLなど無料で入手できるものが存在し、直接比較したことこそないものの、計算速度という観点からも信頼性という観点からもLibMK4はそれらのルーチンに遠く及ばないと思われるためである。

それではなぜ、筆者自身は自作の積分ルーチンを使い続けるのか。第一の理由はweb上で見つけたルーチンで関数オブジェクトの積分が行えるものを見たことがないためである。関数オブジェクトがもたらす柔軟性は慣れてしまうと手放せないため、この条件は筆者にとっては必須条件である。
また、既存の積分ルーチンのインターフェースで好きになれるものがなかったことも挙げられる。積分ルーチンが関数として提供される場合に(積分ルーチンがクラスとして提供される可能性に関する考察は後述)、テンプレートおよび関数のオーバーロードの機能を用いたインターフェースの簡略化はなされて然るべきだと思うのだが、ほとんどの積分ルーチンは汎用性のためかcで書かれているため、これらの機能が使われていない。

逆に、LibMK4の積分ルーチンはこれらの要求を満たすために書かれたものであり、それこそが以下で紹介するルーチン達の存在意義でもある。

積分アルゴリズム

LibMK4が提供する有限閉区間の積分アルゴリズムは、大きく分けて以下の二つがある。
(1)IntegR
スタンダードなRomberg積分。アルゴリズムとしてはNumerical RecepiesのRomberg積分と同じである。ただし、明らかに非効率な部分には一通りの最適化を施してあるのでNumerical Recepiesのものよりは速い。通常の積分ならばこのルーチンで問題ない。
(2)IntegA
積分領域の一部でのみ急激な変化をする、もしくは滑らかでない(不連続もしくは高次の微係数に不連続がある)被積分関数を想定したルーチン。Romberg積分がベースだが、積分実行中に積分区間を分割し、収束性の悪い部分のみを重点的にサンプリングすることで被積分関数の計算回数を減らし、効率を上げる。被積分関数によっては、Romberg積分よりこちらの方が圧倒的に効率が良い。

これら双方にそれぞれ、一変数および多変数関数用のルーチン、また積分区間の端点での発散を許すものと許さないもの、の組み合わせで合計四通りが存在する。

これら有限閉区間用のルーチンに加え、以下のようなものも実装されている。
(1)積分区間の上端が無限大の積分。
(2)被積分関数が、積分区間内の一点で発散する場合の積分。ただし当然ながら積分結果は有限の値でなくてはならない。
(3)主値積分( P\int dx f(x)/(a-x) )を行うための積分ルーチン。

関数リスト

IntegR

端点を含む積分区間に発散を含まない被積分関数用の、スタンダードなRomberg積分。ただし簡略化のため、Rombergの次数は５に固定してある。LibMK4/integral.hppに実装されている。

使用例：

#include