積分式

以下の式を計算する

\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)を変更することで分点の数を変更できる。

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

積分式

サンプルコード

実行例