数値積分 > ガウス積分 / gaussian quadrature > n点の積分例 [対数積分]

[ガウス積分](http://qiita.com/7of9/items/4aeb1b04dab9ce2fdf1f)の計算を対数積分に対応した。

ついでに[sunil](http://www.physics.iitm.ac.in/~sunil/differentiate-integrate.html)さんのweightsとabscissa計算も取り込んで任意の点で計算できるようにしている。
(コード掲載に問題があれば取り下げます。)

#積分式
以下の式を計算する

```math
\int_{10^{-4}}^{10^0}{\pi dr} = \int_{ln(10^{-4})}^{ln(10^0)}{\pi r dlnr}
```
結果はほぼπの値になるはず

#サンプルコード

```integ_gauss.c
#include <stdio.h>
#include <math.h>

/*
 * to use log scale
 * $gcc integ_gauss.c -DUSE_LOG_SCALE -lm
 * 
 * to use normal scale
 * $gcc integ_gauss.c -lm 
 */

//#define USE_LOG_SCALE

/* define paramter constants */
#define EPS 3.0e-11

static void gauleg(double x1, double x2, double x[], double w[], int n);

#define SIZBUF (100)
const int kNumGauss = 11;
static double s_weight[SIZBUF];
static double s_abscissa[SIZBUF];

const double kRmin =  1.e-4;
const double kRmax =  1.e0;

double myFunc_dr(double r_v)
{
	// d[N(r_v)] / d[r_v]
	double pi = acos(-1.0);
	return pi;
}

double myFunc_dlnr(double r_v)
{
	// d[N(r_v)] / dln[r_v]
	return myFunc_dr(r_v) * r_v;
}

double wrapperFunc(double x_)
{
#ifdef USE_LOG_SCALE
	return myFunc_dlnr( exp(x_) );
#else // USE_LOG_SCALE
	return myFunc_dr(x_);
#endif // USE_LOG_SCALE	
}

int main(void)
{
	int idx;
	double d0;

#ifdef USE_LOG_SCALE
	printf("in log scale\n");
	gauleg(log(kRmin), log(kRmax), &s_abscissa[0], &s_weight[0], kNumGauss);
#else
	printf("in normal scale\n");
	gauleg(kRmin, kRmax, &s_abscissa[0], &s_weight[0], kNumGauss);
#endif	

	d0 = 0.0;
	for(idx=1; idx<=kNumGauss; idx++) {
		printf("%f %f dbg\n", s_abscissa[idx], s_weight[idx]);
		d0 += ( s_weight[idx] * wrapperFunc(s_abscissa[idx]) );
	}
	printf("%f\n", d0);
}

// from sunil -------------------------------------------------

/* the function to genearate the abscissa and weights using 
Gauss-Legandere  qudrature rules */
static void gauleg(double x1, double x2, double x[], double w[], int n)
{
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;

	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
/* find m zeros of the Legandere polynomial in the interval -1 to +1 */
	for (i=1;i<=m;i++) {
		z=cos(3.141592654*(i-0.25)/(n+0.5));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (fabs(z-z1) > EPS);
/* rescale the values to the limits required */
		x[i]=xm-xl*z;
		x[n+1-i]=xm+xl*z;
/* calculate the weight */
		w[i]=2.0*xl/((1.0-z*z)*pp*pp);
		w[n+1-i]=w[i];
	}
}
#undef EPS
```

#実行例
実行方法はソースのヘッダに記載している通り。

normalスケールで積分する場合
```
$ gcc integ_gauss.c -lm
$ ./a.out
in normal scale
0.010985 0.027832 dbg
0.056563 0.062784 dbg
0.135011 0.093136 dbg
0.240528 0.116585 dbg
0.365292 0.131389 dbg
0.500050 0.136449 dbg
0.634808 0.131389 dbg
0.759572 0.116585 dbg
0.865089 0.093136 dbg
0.943537 0.062784 dbg
0.989115 0.027832 dbg
3.141278
```

対数スケールで積分する場合
```
$ gcc integ_gauss.c -DUSE_LOG_SCALE -lm
$ ./a.out
in log scale
-9.110080 0.256363 dbg
-8.690244 0.578319 dbg
-7.967644 0.857898 dbg
-6.995696 1.073897 dbg
-5.846462 1.210260 dbg
-4.605170 1.256866 dbg
-3.363878 1.210260 dbg
-2.214644 1.073897 dbg
-1.242696 0.857898 dbg
-0.520096 0.578319 dbg
-0.100261 0.256363 dbg
3.141278
```

両者とも3.1412...とほぼ同じ結果となった。

kNumGauss (=11)を変更することで分点の数を変更できる。