シンプソン法 > テーブルデータでの積分 [対数積分]

[【C言語で数値計算】 定積分の近似計算（シンプソン法）](http://qiita.com/HirofumiYashima/items/6636bd3a8ce0460e5a4f)
[シンプソン法 > テーブルデータでの積分](http://qiita.com/7of9/items/5b6f4e7509368efa9b2d)

#計算について
色々な物事は対数スケールで変化していく。

上記のリンクを元に対数スケールでの積分をするようにコードを変更してみた。

対数スケールにする時
t = ln(r_v) -- (1)と置き
(1)より、r_v = exp(t) -- (2)
(1)を微分して、dt/dr_v = 1/r_v -- (3)
(3)よりdr_v=r_v * dtとなる

例として
1/r_v/r_v * pi * r_v * r_v -- (4)
なる分布関数を <font color=red>[10^-4, 10^0]の範囲</font>で積分するとする。

normalスケールで積分する場合、上記の分布関数は (4) = pi をdr_vで積分をすればいいので、[pi * r_v]となり
pi(1 - 10^(-4)) >> およそ 3.14となる。

logスケール積分する場合、(4)は(3)より [pi * r_v]をdln(r_v)で積分となる。
(2)よりr_v = exp(t)となるので、[pi * exp(t)]のdln(r_v)での積分となる。
pi * [exp(ln(10^0)) - exp(ln(10^-4))] となり、こちらもおよそ3.14となる。

#計算コード

シンプソン法の計算をするためのテーブルを作成するコードとテーブルを元に積分するコードに分けている。

1. prep_table.c : テーブル作成
2. calc_integ.c : テーブルを元に積分する

## prep_table.c
```c:prep_table.c
#include <stdio.h>
#include <math.h>

#define USE_LOG_SCALE

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

//  (result is 3.14...)
	res = (pi * r_v * r_v) / r_v / r_v;
	return res;
}

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

int main(void) 
{
	//----- configuration -----
	// for simpson integral
	const double kRmin = 1.e-4;
	const double kRmax = 1.e0;
	const int kDivnum = 120;
	const int kLastIdx = 2*kDivnum;
	// other
	double pi = acos(-1.0);

	//----- local variables -----
	double r_v, log_rv;
	int idx;
	double xs, ys;
	double hsimp; // width of simpson integral

	printf("%d :kLastIdx\n", kLastIdx);

	//
#ifdef USE_LOG_SCALE	
	// for integartion by log scale
	hsimp = (log(kRmax) - log(kRmin)) / (2.0*(double)kDivnum);
	for(idx=0; idx<=kLastIdx; idx++) {
		log_rv = log(kRmin) + (double)(idx)*hsimp;
		r_v = exp(log_rv);

		xs = log_rv;
		ys = myFunc_dlnr(r_v);
		printf("%le %le :xs,ys%d\n", xs, ys, idx);
	}
#else // USE_LOG_SCALE
	// for integraty by normal scale
	hsimp = (kRmax - kRmin) / (2.0*(double)kDivnum);
	for(idx=0; idx<=kLastIdx; idx++) {
		r_v = kRmin + (double)(idx)*hsimp;

		xs = r_v;
		ys = myFunc_dr(r_v);
		printf("%le %le :xs,ys%d\n", xs, ys, idx);
	}
#endif // USE_LOG_SCALE

}
```

## calc_integ.c
```c:calc_integ.c
#include <stdio.h>
#include <stdbool.h>
#include <math.h>

char *readFilename = "simp.tbl";

static double simpson_integ_table(double *pxs, double *pys, int lastIdx_)
{
	double integral;
	double hsimp; // width (when divided with n intervals)
	int idx;

	hsimp = (pxs[lastIdx_] - pxs[0]) / (double)lastIdx_;
	integral = pys[0];
	for (idx=1; idx<lastIdx_-2; idx+=2) {
		integral += 4.0*pys[idx] + 2.0*pys[idx+1];
	}
	integral += ( 4.0*pys[idx] + pys[idx+1] );
	integral *= hsimp/3.0; 

	return integral;
}

#define BUF_SIZ (5000)
static double s_xs[BUF_SIZ] = {0.0};
static double s_ys[BUF_SIZ] = {0.0};
int s_lastIdx;

static void readFile_dlnr(void)
{
	FILE *fp;
	int idx;
	char dummy[200];
	int ret;
	float fx, fy;
	double r_v;
	double pi = acos(-1.e0);

	fp = fopen(readFilename, "r");
	ret = fscanf(fp, "%d %s", &s_lastIdx, dummy);

	for(idx=0; idx<=s_lastIdx; idx++) {
		ret = fscanf(fp, "%le %le %s", &s_xs[idx], &s_ys[idx], dummy);
	}
	fclose(fp);
}

int main(void)
{
	double denom;

	readFile_dlnr();
	denom = simpson_integ_table(s_xs, s_ys, s_lastIdx);
	printf("%lf\n", denom);

	return(0);
}
```

## 使い方

```txt:動作環境
CentOS 6.5
gcc v4.4.7
```

```txt:実行例
$ gcc prep_table.c -lm -o prep_table
$ gcc calc_integ.c -lm -o calc_integ
$ ./prep_table > simp.tbl
$ ./calc_integ
3.141278
```

prep_table.c内にて`USE_LOG_SCALE`デバッグオプションにより[対数積分]と[normal積分]を切り替える。
どちらも3.1412...程度となる。

`kDivnum`にて積分の分割数を指定する。