数値積分(シンプソン則)

<?php

/****************************/
/* シンプソン則             */
/*      coded by Y.Suganuma */
/****************************/

	$pi2 = 2.0 * atan(1.0);

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

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

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

/*******************************************************/
/* 数値積分(シンプソン則)                            */
/*      x1 : 下限                                      */
/*      x2 : 上限                                      */
/*      n : 分割数                                     */
/*      fun : 関数名(f)                              */
/*      return : 積分値                                */
/*******************************************************/
function simpson($x1, $x2, $n, $fun)
{
	$s1 = 0.0;
	$s2 = 0.0;
	$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;
}

?>