本
稿
目
次
稿
目
次
本
稿
目
次
稿
目
次
前提
ある関数
について区間
の積分を考えます。
区間
を n 等分し、その分点をそれぞれ
とし、
がそれぞれ分かっていて
だとします。つまり、
です。
また、各分点の幅は h とします。
関数をラグランジュ補間公式で近似
関数
をラグランジュ補間公式で近似的に置き換えます。
すると求めたい面積 S は
となります。
ニュートン・コーツの積分公式
前述の積分を行うために
と置いて置換を行います。また、前提で述べた定義により
とおけます。すると
です。また、このとき
となり、積分区間は
となります。以上をまとめると
ここで
とすれば、
ですので、あらかじめこの C を求めておけば面積を求めることができます。
係数 C の各項は展開するとわかるとおり、n 次多項式になっています。従ってすべての項を足しても n 次多項式です。これは n + 1 点からなる点列を n 次式で近似することによって面積を求めていることに他なりません。つまり、2 点からなる点列なら 1 次式（直線）で近似し、3 点からなる点列なら 2 次式で近似させるわけです。
係数 C と台形公式、シンプソンの公式
台形公式
次数 n が低ければ直接求めることができます。
次数 n が 1 のときを特に
同様にして
よって面積は
実際にもっと多くの分点がある場合はこれを繰り返して
となります。
シンプソンの公式
次数 n が 2 のときを特に
同様にして
よって面積は
台形公式のときと同じように、分点が複数あるときはこれを繰り返して
一般的にはこのシンプソンの公式がよく使われます。
注意
シンプソンの公式のように2次以上の n 次式で近似する場合、分割数が n の倍数（分点数が n の倍数プラス 1）でなければなりません。そうでなければ、残った部分は切り捨てたり、そこだけ n 次未満で求めるなどの工夫が必要でしょう。
サンプルプログラム
ニュートン・コーツ積分関数
// ver.041107
#ifndef __NEWTON_COTES_H
#define __NEWTON_COTES_H
#include <vector>
// y データの列
// h データの間隔
// n 次数（n=1 なら台形公式、n=2 ならシンプソンの公式、1 <= n <= 5 に対応）
double NewtonCotes(std::vector<double> &y, double h, long n)
{
if(n < 1) n = 1;
if(n > 5) n = 5;
static bool isFirst = true;
static double cTable[6][6]; //[0][] は未使用
if(isFirst){
cTable[1][0] = 1.0 / 2.0;
cTable[1][1] = 1.0 / 2.0;
cTable[2][0] = 1.0 / 3.0;
cTable[2][1] = 4.0 / 3.0;
cTable[2][2] = 1.0 / 3.0;
cTable[3][0] = 3.0 / 8.0;
cTable[3][1] = 9.0 / 8.0;
cTable[3][2] = 9.0 / 8.0;
cTable[3][3] = 3.0 / 8.0;
cTable[4][0] = 14.0 / 45.0;
cTable[4][1] = 64.0 / 45.0;
cTable[4][2] = 24.0 / 45.0;
cTable[4][3] = 64.0 / 45.0;
cTable[4][4] = 14.0 / 45.0;
cTable[5][0] = 95.0 / 288.0;
cTable[5][1] = 375.0 / 288.0;
cTable[5][2] = 250.0 / 288.0;
cTable[5][3] = 250.0 / 288.0;
cTable[5][4] = 375.0 / 288.0;
cTable[5][5] = 95.0 / 288.0;
isFirst = false;
}
int m = (int)y.size() - 1; //分割数
double s = 0;
long k, i;
for(k = 0; k <= m / n - 1; k++){
for(i = 0; i <= n; i++){
s += cTable[n][i] * y[n * k + i];
}
}
//端数部分の計算
int lastBlockStart = n * k;
int lastBlockEnd = (int)y.size() - 1;
m = lastBlockEnd - lastBlockStart;
n = m;
if(m > 0){
for(i = 0; i <= n; i++){
s += cTable[n][i] * y[lastBlockStart + i];
}
}
return h * s;
}
#endif
メインプログラム
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
#include "../NewtonCotes.h"
main()
{
vector<double> y;
double h = 0.01;
for(int i = 0; i < 314; i++){
y.push_back(sin(i * h));
}
cout << NewtonCotes(y, h, 2) << endl;
return 0;
}
参考
数値計算情報処理入門コース (7)
http://www.amazon.co.jp/exec/obidos/ASIN/4000078577/birdport-22
