// ==================================================================
//  GaussLegendre.sci: ガウスの求積法による数値積分を行う関数
// ------------------------------------------------------------------
//  ＜引数＞   f: 被積分関数
//             a: 積分区間の始まり
//             b: 積分区間の終わり
//             N: 補間点数 (但し，N は 2,3,4,5 の場合にのみ対応)
//  ＜返り値＞ I: 積分値
// ==================================================================
function I = GaussLegendre(f, a, b, N)
xi = [-sqrt(1/3)  -sqrt(3/5)  -0.8611363116  -0.9061798459  // p.53 の表3.3の補間点
       sqrt(1/3)          0    0.3399810436  -0.5384693101
              0    sqrt(3/5)   0.3399810436              0
              0           0    0.8611363116   0.5384693101
              0           0               0   0.9061798459];
 
wi = [        1         5/9    0.3478548451   0.2369268851  // p.53 の表3.3の重み
              1         8/9    0.6521451549   0.4786286705
              0         5/9    0.6521451549   0.5688888889
              0           0    0.3478548451   0.4786286705
              0           0               0   0.2369268851];

x = xi(1:N, N-1);                  // [1,-1]の区間での「補間点」の行列 xi から必要な部分を抽出
x = ((b - a) * x + b + a) / 2;     // [a,b]の区間での「補間点」に変換．p.51のξの逆変数変換
y = f(x);                          // y = [f(x_1), f(x_2),..., f(x_N)]
w = wi(1:N, N-1);                  // 「重み」の行列 wi から必要な部分を抽出
I = (b - a) / 2 * (w' * y);        // 式(3.20)を (b-a)/2 倍して [a,b]の区間に変換
endfunction