数値積分 > ガウス積分 / gaussian quadrature > 5点の積分例

ガウス積分という積分方法について。

概要

積分のポイントは以下のテーブルとして用意されている。
Gaussian Quadrature Weights and Abscissae

計算でもWeights and Abscissaeを求めることができるはず。
sunilさんのページにgaus_leg.cとしてweightsとabscissaeを計算するコードがある。
但し、gaus_leg.c内の関数gauleg()で得られるweight w[]とabscissa x[]はインデックスが1から使う点に注意。

サンプルコード

例として5点のガウス積分を実装してみた。
積分式は以下のものとした。

\int_{-1}^{1}{\pi}dr = [\pi r]^{1}_{-1} = 2\pi

sample-gauss.c

# include <stdio.h>
# include <math.h>

# define IDX_WEI (0)
# define IDX_ABS (1)
# define NUM_ITEMS (5)

double cf[NUM_ITEMS][2] =
{
// weight, abscissa
{0.5688888888888889,	0.0000000000000000},
{0.4786286704993665,	-0.5384693101056831},
{0.4786286704993665,	0.5384693101056831},
{0.2369268850561891,	-0.9061798459386640},
{0.2369268850561891,	0.9061798459386640	}
};

static double myFunc(double x_)
{
	const double pi=acos(-1.0);
	return pi;
}

int main(void)
{
	double d0, xx, yy, weight;
	int idx;

	d0 = 0.0;
	for(idx=0; idx<NUM_ITEMS; idx++) {
		xx = cf[idx][IDX_ABS];
		weight = cf[idx][IDX_WEI];
		yy = myFunc(xx);
		d0 += (weight * yy);
	}

	printf("%lf\n", d0); // 2pi
}

実行結果

$ gcc sample-gauss.c -lm
$ ./a.out
6.283185