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

【C言語で数値計算】 定積分の近似計算（シンプソン法）
をもとに、テーブルのデータを積分するコードに変更した。

実際にシンプソン積分をしたいものとして、時間のかかる計算結果を用いることがある。その際に先に結果をテーブルとして出力しておき、それを後で積分するということに使う。

simpson_integ_table.c

#include <stdio.h>

// 積分対象関数.
float function (float x) 
{
    float result;
    result = x; // 一次直線
    return result;
}

static double simpson_integ_table(double *xs, double *ys, int lastIdx)
{
    double integral; // 積分結果
    double h; // 積分範囲を n 個に分割したときの幅
    double dA;
    int idx;

    h = (xs[lastIdx] - xs[0]) / lastIdx; // 分割幅

    integral = ys[0]; // 積分結果の変数の初期値

    for (idx=1; idx<lastIdx-2; idx+=2) {
        dA = 4.0*ys[idx] + 2.0*ys[idx+1];
        integral += dA; // Σ
    }

    integral += ( 4.0*ys[idx] + ys[idx+1] );

    integral *= h/3.0; 

    return integral;
}

#define BUF_SIZ (500)

int main()
{
    double integral;
    double xs[BUF_SIZ]; // x
    double ys[BUF_SIZ]; // y
    double h;
    int count;
    int idx;
    double dval;

    // configuration
    double xmin = 0.0;
    double xmax = 100.0;
    int divnum = 120;
    int lastIdx = 2*divnum;

    // 1. preparation of tables
    h = (xmax - xmin) / (2.0*(double)divnum); // 分割幅
    for(idx=0; idx<=lastIdx; idx++){
        dval = xmin + (double)(idx) * h;
        xs[idx] = dval;
        ys[idx] = function(dval);
    }

    // 2. integration
    integral = simpson_integ_table(xs, ys, lastIdx);

    // === Simpson 法による積分 (終了） ===
    printf ("積分結果= %lf\n", integral);
    printf ("解析的に求めた結果 (底辺 100，高さ 100 の直角三角形)=%lf\n",
        (xmax - xmin)*function(xmax) / 2.0 );

    return(0);
}