補間法(ラグランジュ)

<?php

/****************************/
/* 補間法(ラグランジュ)   */
/*      coded by Y.Suganuma */
/****************************/

	$x    = array(3);
	$y    = array(3);
					// 線形補間
	$n    = 1;
	$x[0] = 0.0;
	$x[1] = 2.0;
	$y[0] = 0.0;
	$y[1] = 4.0;

	printf("線形補間\n");
	$x1 = 0.0;
	for ($i1 = 0; $i1 <= 8; $i1++) {
		$y1 = Lagrange($n, $x, $y, $x1);
		printf("   x %f y %f\n", $x1, $y1);
		$x1 += 0.25;
	}

					// 2次補間
	$n    = 2;
	$x[0] = 0.0;
	$x[1] = 1.0;
	$x[2] = 2.0;
	$y[0] = 0.0;
	$y[1] = 1.0;
	$y[2] = 4.0;

	printf("2次補間\n");
	$x1 = 0.0;
	for ($i1 = 0; $i1 <= 8; $i1++) {
		$y1 = Lagrange($n, $x, $y, $x1);
		printf("   x %f y %f\n", $x1, $y1);
		$x1 += 0.25;
	}

/**********************************/
/* ラグランジュ補間法             */
/*      n : 次数(n=1は線形補間) */
/*      x, y : (n+1)個の点の座標  */
/*      x1 : 補間したい点のx座標 */
/*      return : 補間結果         */
/**********************************/
function Lagrange($n, $x, $y, $x1)
{
	$y1 = 0.0;

	for ($i1 = 0; $i1 <= $n; $i1++) {
		$s1 = 1.0;
		$s2 = 1.0;
		for ($i2 = 0; $i2 <= $n; $i2++) {
			if ($i2 != $i1) {
				$s1 *= ($x1 - $x[$i2]);
				$s2 *= ($x[$i1] - $x[$i2]);
			}
		}
		$y1 += $y[$i1] * $s1 / $s2;
	}

	return $y1;
}

?>