メインコンテンツに移動

やっぱ積分ですよね

GSLで数値積分

GSLでは, いろいろな数値積分スキームが準備されています.

Gauss-Legendre積分

基本中の基本, ガウス求積法である. 古いものが好きな人が使うらしい.

Gauss-Kronrod積分

ガウス求積法で, 区分点を入れ子型に選ぶ工夫を行うと, 関数の評価点数を小さくできます. これがGauss-Kronrod積分だそうです. 

QNG 基本の積分法

QAG 誤差評価を行って，関数の変化に追従する適応型積分.

QAGS QAGにおいて特異点の自動検出を試みるQAG

QAGP 特異点リストを与えるQAG

QAGI 変数変換 x=(1-t)/t などを用いて無限積分, 半無限積分を行うQAG

QAWC コーシーの主値積分を行うQAG

QAWS 端点においてべきおよび対数形特異点を持つ関数のQAG

QAWO 周期関数cos, sinを含む関数のQAG

QAWF フーリエ変換型の積分に対するQAG

二重適応型積分

CQUAD なんだか知らんが, 精度がいいらしい.

Romberg積分

  • なんだかよくわからんが, 積分するらしい

なぜか, スーパー性能という京都大学数理研の二重指数型積分は無いんですね・・・あ？変数変換したら，台形則が超高性能になるという理論だから，ライブラリーで自動積分とかいう話ではないんですね．これは・・・お？お弟子さんのWebページになんか寂れているので, そいつも合わせて説明しよう.

京風数値積分

添付のIntDEクラスにより, 二重指数型積分が可能です. GSLとは関係がありません.

#include "intde.hpp"
struct my_type {
double K=1.0;
} myParameter;
IntDE IDE;
IDE.Parameter=&myParameter; //パラメータのポインタを指定してもよい
IDE.Initialize(1e-15);//ご希望の精度で初期化
double xmin=0.0,xmax=0.5, err; 
double v=IDE.Integrate([](double x, void *p){
	my_type& P=(my_type *)p;//パラメータのポインタを指定したなら，復元できる
    return P.K/std::sqrt(x*(1.-x));
}, xmin, xmax, err)*2;

格別説明は不要だと思います. 上の例では被積分関数をラムダ式で書いていますが, もちろん, グローバルな関数の関数ポインタでもOKです. 関数の他にもいくつかあって,

  • class IntDEI(err_rel)
    • Integrate([](double x){...},double a, double & err)   [a,∞)半無限積分を行う
  • class IntDEO(err_rel)
    • Integrate([](double x){...},double a, double omega, double & err)  f(x)は振動項を含み,xの大きいところで漸近的に(単調な関数)* sin(omega*x+theta)となっている

格子点(デフォルト8000)が不足する場合, コンストラクタでIntDE(n) と指定して増やせば良い.

errは推定絶対誤差で収束に失敗する(長さnの表で許容誤差に 達しない)と負で返す. このときの原因と対策は

  1. f(x)に積分端点以外で微分不可能や不連続な点,あるいは それに近い点(鋭いピークなど)がある. この場合そのような点で積分区間を分割すれば計算可能.
  2. f(x)の計算に大きな誤差が含まれる. この場合桁落ちなどが起きないようにf(x)を変形すればよい. ただし端点で発散するような積分の場合はさらに特異点を原点にずらしておく必要がある.
    • これは意外と重要です．上の例では 2/sqrt(x*(1-x)) を(0,1/2]で積分しているので結果はπですが, 実際に実行すると3.14159265358979が得られます. ところが, 1/sqrt(x*(1-x) を(0,1)で積分すると, NANになってしまい積分できません. 
    • 当然ですが, x=1で特異性を持つ場合, g(x)-g(1) という形が含まれてはならないんで, そうなるように式変形するんですよ？変数の定義をずらして特異点をx=0にして, 第0項は解析でキャンセルdevilキャンセルdevil のこりを h(x)*x という形にできたらグッドだにゃー
  3. intdeo以外でf(x)が激しく振動している. この場合はintdeoが使える形に変形できれば計算可能. 変形できない複雑な振動積分は別のルーチンが必要. 積分の和の計算で桁落ちが起きると, 実際の相対誤差は epsで指定した精度よりも悪くなる. 例えばeps=1.0e-15と 指定したときに5桁の桁落ちが起きた場合10桁の精度になる. 誤差の確認はerrを参照すること.
  4. 推定誤差errはあくまで推定であり厳密な誤差の上限ではない. 計算が遅くても高い信頼性を要求する場合にはパッケージ中の efsの値を小さめ(例えばefs=sqrt(eps)程度)に選ぶとよい.

QAGS

GSLのQAGS積分を行ってみましょう.

#include 
struct my_type {
double K=1.0;
} myParameter;
gsl_function QAGS_f;
QAGS_f.Parameter=&myParameter; //パラメータのポインタを指定してもよいunc
QAGS_f.function=[](double x, void *p){
	my_type& P=(my_type *)p;//パラメータのポインタを指定したなら，復元できる
    return P.K/std::sqrt(x*(1.-x));
};
gsl_integration_workspace *QAGS_ws=gsl_integration_workspace_alloc (1000);
double xmin=0.0,xmax=0.5, err; 
double v;
gsl_integration_qags(&QAGS_f,xmin,xmax,0.0,1.e-12,1000,QAGS_ws,&v,&err);

格別説明は不要ですよね？

QAGP

GSLのQAGP積分を行ってみましょう. これは, 特異点が既知の場合に有用です. 特異点には, 無条件で積分範囲の両端を入れなければなりません.

#include 
struct my_type {
double K=1.0;
} myParameter;
gsl_function QAGP_f;
QAGP_f.Parameter=&myParameter; //パラメータのポインタを指定してもよいunc
QAGP_f.function=[](double x, void *p){
	my_type& P=(my_type *)p;//パラメータのポインタを指定したなら，復元できる
    return P.K/std::sqrt(x*(1.-x));
};
std::vector<double> QAGP_pts;
QAGP_pts.resize(2);
QAGP_pts.front()=0.0;QAGP_pts.back()=0.5;
double err;
gsl_integration_workspace *QAGP_ws=gsl_integration_workspace_alloc (1000);
  double v;
gsl_integration_qagp(&QAGS_f,QAGP_pts.data(),QAGP_pts.size(),0.0,1.e-12,1000,QAGS_ws,&v,&err);

格別説明は不要ですよね？

 

スピード競走

いろいろ公式があって，どれにしようかな，というところですね．そこで，スピード競走例を示しましょう．

例題：2つのLennard-Jones分子が2体衝突を行うとき, その散乱角θは, 以下の積分で与えられる.

f(x) = 1 - a x + b x^3 (x^3-1),     x>0であり, a, bは非負定数である. 

xc=min{ x | f(x)=0, x>0}

θ=[sqrt(a)/2] int_0^xc  dx/sqrt(x f(x))

被積分関数は x=0, xc で特異点を持つ. 運が悪いと(b>5で, aの値が，とある関数a(b)の値と一致すると), 積分が発散する(惑星の周りでスイングばいwink). まあ普通に選んでたら当たらないけど, ギリギリまずいケースもあるかも，ですheart

これを2268000点の(a,b)で行った場合のCPU時間は,

積分ライブラリ

CPU(sec)
MacbookPro
debug/release

 結果θの平均 err NAN
/
SIGABRT
QAGS 41.897 / 27.662 0.589963 1.19e-13 0 / 1
QAGP 41.553 /  27.526 0.589963 1.19e-13 0 / 0
IntDE 12.350 /  8.363 0.589925 1.0e-13 16  / 0
QAGS 31.450 / 27.144 0589963 1.5e-10 0 /  0
QAGP 29.886 / 18.538 0.589963 1.5e-10 0 / 0
IntDE 10.169 /  7.396 0.589952 1.0e-9 4 / 0

なお, errは, ライブラリーが自分で推定した誤差の言い値の平均ですから, どこまで信じられるかは不明. SIGABRTはテスト中に計算が止まってしまった総数(出たケースでは, パラメータをずらして回避). 

  • IntDEが3倍以上早いが, 特異点にちょいと弱い感じですね. isnan(結果変数) で検出可能ですが．
    • まあ被積分関数を漸近展開するとかして，テキトーに回避してください
  • IntDEでは破局ケース(SIGABORTによる計算停止)は出なかった．

ついでにLinuxでも計算時間を見てみると(ic++ -Ofast --gsl --boost -std=c++11)

積分ライブラリ

is2015e
Corei7-5930K

 ig2016b
XeonE5-269

ig2017a 
XeonE5-2695

 ib2014b 
Corei7-4790K
QAGS *1  *1 *1 *1
QAGP *1 *1 *1 *1
IntDE 4.0069 sec 3.547 sec 3.483 sec 3.253 sec

*1 このやろう, エラー吐きまくり動かねえ・・・roundoffエラーに引っかかるってことらしいが

さすがDevil's Canyonで鳴らした4790K, 老いたりとはいえ負けてねえな．
Linuxにコピーする際に, 最初isnan()だったら特異点前/中/後の3領域に分けて積分し直すように改造したらNANが消えた

やっぱり積分は和風の味が似合うようだな．

添付 サイズ
intde.zip 23.02 KB
intde++.zip 6.39 KB