積分の方法は色々あるが、その中でQMCという積分の紹介。
概要
10年近く前にlucilleというページで初めて知った。
Prof. Harald Niederreiterさんの「Random Number Generation and Quasi-Monte Carlo Methods (PDF, 11MB)」に詳しく掲載されている。
積分をする上でモンテカルロ法というのがある。積分するポイントを乱数で取得しながら各地点の値を足しあわせていくことで積分を行う方法である。
モンテカルロ法より少ない回数で積分精度をあげようというのがQMCという手法である。積分のオーダーについては上記のPDFを参照。
QMCではlow-discrepancy sequence(LDS)という数列を用いる。その二次元版は「Halton_sequence」という数列と認識している。
サンプルコード
以下のサンプルコードで実際に積分してみる。
Halton_sequence()で二次元のLDSを作っているが、今回は一次元のみ使っている。
なお積分範囲は[0,1]となるので、これ以外の積分の場合は[0,1]で得られるLDSを積分の希望範囲へ拡張すればいい。
#include <stdio.h>
#include <math.h>
static void Halton_sequence(int i0, double *px0, double *py0)
{
// configuration
const int xbase = 2;
const int ybase = 3;
//
double invxbase, invybase;
double facx, facy;
double x0, y0;
int inp;
invxbase = 1.0 / (double)xbase;
facx = 1.0 / (double)xbase;
invybase = 1.0 / (double)ybase;
facy = 1.0 / (double)ybase;
inp = i0;
x0 = 0.0;
while (inp > 0) {
x0 = x0 + (inp % xbase) * invxbase;
inp = inp / xbase;
invxbase = invxbase * facx;
}
inp = i0;
y0 = 0.0;
while (inp > 0) {
y0 = y0 + (inp % ybase) * invybase;
inp = inp / ybase;
invybase = invybase * facy;
}
*px0 = x0;
*py0 = y0;
}
static double myFunc(double x_)
{
// result [1/4]
return (x_ * x_ * x_);
}
const static double kAnswer = 1.0/4.0;
int main(void)
{
double xx, yy;
int idx;
double d0;
int count = 0;
double res;
d0 = 0.0;
for(idx=1; idx<=3000; idx++) {
count++;
Halton_sequence(idx, &xx, &yy);
d0 = d0 + myFunc(xx);
res = d0 / (double)count;
printf("%d %f err=%.2f\%\n", idx, res, fabs(res/kAnswer - 1.0)*100);
}
}
1 0.125000 err=50.00%
2 0.070312 err=71.88%
3 0.187500 err=25.00%
...
300 0.245102 err=1.96%
301 0.245453 err=1.82%
302 0.244952 err=2.02%
...
1000 0.248766 err=0.49%
1001 0.248725 err=0.51%
1002 0.248517 err=0.59%
...
2998 0.249347 err=0.26%
2999 0.249531 err=0.19%
3000 0.249448 err=0.22%
ガウス積分などと比べて同じ精度では計算する点数が多くなってしまう。
QMCの良いと思う所は、点を増やすときは過去の点そのままが使えるという所。
100個で計算したけど収束していない場合、101個目から計算を続ければいい。ガウス積分やシンプソン積分だとこうはいかない。
そのため、収束精度を設定して、その精度が求まるまで計算を続けるという方法が取れそう。