【C言語で数値計算】定積分の近似計算（シンプソン法）
シンプソン法 > テーブルデータでの積分

計算について

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

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

対数スケールにする時
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)
なる分布関数を [10^-4, 10^0]の範囲で積分するとする。

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となる。

計算コード

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

prep_table.c : テーブル作成
calc_integ.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

#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);
}

使い方

動作環境

CentOS 6.5
gcc v4.4.7

実行例

$ 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にて積分の分割数を指定する。

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

計算について

計算コード

prep_table.c

calc_integ.c

使い方