数値積分（シンプソン則）

/****************************/
/* シンプソン則             */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>

double snx(double);
double simpson(double, double, int, double(*)(double));

int main()
{
	double y, pi2 = 2.0 * atan(1.0);

	y = simpson(0.0, pi2, 100, snx);

	printf("result %f\n", y);

	return 0;
}

/****************/
/* 関数値の計算 */
/****************/
double snx(double x)
{
	return sin(x);
}

/*******************************************************/
/* 数値積分（シンプソン則）                            */
/*      x1 : 下限                                      */
/*      x2 : 上限                                      */
/*      n : 分割数                                     */
/*      fun : 関数名（f）                              */
/*      return : 積分値                                */
/*******************************************************/
double simpson(double x1, double x2, int n, double(*fun)(double))
{
	double s1 = 0.0;
	double s2 = 0.0;
	double s, x, h, h2;

	h  = (x2 - x1) / n;
	h2 = 2.0 * h;

	for (x = x1+h; x < x2; x += h2)
		s1 += (*fun)(x);

	for (x = x1+h2; x < x2-h; x += h2)
		s2 += (*fun)(x);

	s = h * ((*fun)(x1) + (*fun)(x2) + 4.0 * s1 + 2.0 * s2) / 3.0;

	return s;
}