補間法（スプライン）

<?php

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/

					// 設定
	$n = 5;
	$h = array($n);
	$b = array($n);
	$d = array($n);
	$g = array($n);
	$u = array($n);
	$r = array($n+1);
	$x = array($n+1);
	$y = array($n+1);

	$x[0] = 0.0;
	$x[1] = 1.0;
	$x[2] = 2.0;
	$x[3] = 3.0;
	$x[4] = 4.0;
	$x[5] = 5.0;
	$y[0] = 0.0;
	$y[1] = 1.0;
	$y[2] = 4.0;
	$y[3] = 9.0;
	$y[4] = 16.0;
	$y[5] = 25.0;
					// 実行と出力
	$x1 = 0.0;
	$sp = 0.1;
	while ($x1 <= 5.01) {
		$y1 = spline($n, $x, $y, $x1, $h, $b, $d, $g, $u, $r);
		printf("%f %f\n", $x1, $y1);
		$x1 += $sp;
	}

/**********************************/
/* ３次スプライン関数による補間   */
/*      n : 区間の数              */
/*      x, y : 点の座標           */
/*      x1 : 補間値を求める値     */
/*      h, b, d, g, u, r : 作業域 */
/*      return : 補間値           */
/*      coded by Y.Suganuma       */
/**********************************/
function spline($n, $x, $y, $x1, $h, $b, $d, $g, $u, $r)
{
	$i = -1;
					// 区間の決定
	for ($i1 = 1; $i1 < $n && $i < 0; $i1++) {
		if ($x1 < $x[$i1])
			$i = $i1 - 1;
	}
	if ($i < 0)
		$i = $n - 1;
					// ステップ１
	for ($i1 = 0; $i1 < $n; $i1++)
		$h[$i1] = $x[$i1+1] - $x[$i1];
	for ($i1 = 1; $i1 < $n; $i1++) {
		$b[$i1] = 2.0 * ($h[$i1] + $h[$i1-1]);
		$d[$i1] = 3.0 * (($y[$i1+1] - $y[$i1]) / $h[$i1] - ($y[$i1] - $y[$i1-1]) / $h[$i1-1]);
	}
					// ステップ２
	$g[1] = $h[1] / $b[1];
	for ($i1 = 2; $i1 < $n-1; $i1++)
		$g[$i1] = $h[$i1] / ($b[$i1] - $h[$i1-1] * $g[$i1-1]);
	$u[1] = $d[1] / $b[1];
	for ($i1 = 2; $i1 < $n; $i1++)
		$u[$i1] = ($d[$i1] - $h[$i1-1] * $u[$i1-1]) / ($b[$i1] - $h[$i1-1] * $g[$i1-1]);
					// ステップ３
	$k       = ($i > 1) ? $i : 1;
	$r[0]    = 0.0;
	$r[$n]   = 0.0;
	$r[$n-1] = $u[$n-1];
	for ($i1 = $n-2; $i1 >= $k; $i1--)
		$r[$i1] = $u[$i1] - $g[$i1] * $r[$i1+1];
					// ステップ４
	$xx = $x1 - $x[$i];
	$qi = ($y[$i+1] - $y[$i]) / $h[$i] - $h[$i] * ($r[$i+1] + 2.0 * $r[$i]) / 3.0;
	$si = ($r[$i+1] - $r[$i]) / (3.0 * $h[$i]);
	$y1 = $y[$i] + $xx * ($qi + $xx * ($r[$i]  + $si * $xx));

	return $y1;
}

?>
		

以下は，クラスを使用して書いた例です．

<?php

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/

/****************/
/* クラスの定義 */
/****************/
class Spline
{
	private $r = array();
	private $q = array();
	private $s = array();
	private $x = array();
	private $y = array();
	private $n;

	/**************************/
	/* コンストラクタ         */
	/*      n1 : 区間の数     */
	/*      x1, y1 : 点の座標 */
	/**************************/
	function Spline($n1, $x1, $y1)
	{
					// 領域の確保
		$this->n = $n1;

		$h = array();
		$b = array();
		$d = array();
		$g = array();
		$u = array();

		for ($i1 = 0; $i1 <= $this->n; $i1++) {
			$this->x[$i1] = $x1[$i1];
			$this->y[$i1] = $y1[$i1];
		}
					// ステップ１
		for ($i1 = 0; $i1 < $this->n; $i1++)
			$h[$i1] = $this->x[$i1+1] - $this->x[$i1];
		for ($i1 = 1; $i1 < $this->n; $i1++) {
			$b[$i1] = 2.0 * ($h[$i1] + $h[$i1-1]);
			$d[$i1] = 3.0 * (($this->y[$i1+1] - $this->y[$i1]) / $h[$i1] - ($this->y[$i1] - $this->y[$i1-1]) / $h[$i1-1]);
		}
					// ステップ２
		$g[1] = $h[1] / $b[1];
		for ($i1 = 2; $i1 < $this->n-1; $i1++)
			$g[$i1] = $h[$i1] / ($b[$i1] - $h[$i1-1] * $g[$i1-1]);
		$u[1] = $d[1] / $b[1];
		for ($i1 = 2; $i1 < $this->n; $i1++)
			$u[$i1] = ($d[$i1] - $h[$i1-1] * $u[$i1-1]) / ($b[$i1] - $h[$i1-1] * $g[$i1-1]);
					// ステップ３
		$this->r[0]          = 0.0;
		$this->r[$this->n]   = 0.0;
		$this->r[$this->n-1] = $u[$this->n-1];
		for ($i1 = $this->n-2; $i1 >= 1; $i1--)
			$this->r[$i1] = $u[$i1] - $g[$i1] * $this->r[$i1+1];
					// ステップ４
		for ($i1 = 0; $i1 < $this->n; $i1++) {
			$this->q[$i1] = ($this->y[$i1+1] - $this->y[$i1]) / $h[$i1] - $h[$i1] * ($this->r[$i1+1] + 2.0 * $this->r[$i1]) / 3.0;
			$this->s[$i1] = ($this->r[$i1+1] - $this->r[$i1]) / (3.0 * $h[$i1]);
		}
	}

	/******************************/
	/* 補間値の計算               */
	/*      x1 : 補間値を求める値 */
	/*      return : 補間値       */
	/******************************/
	function value($x1)
	{
		$i = -1;
					// 区間の決定
		for ($i1 = 1; $i1 < $this->n && $i < 0; $i1++) {
			if ($x1 < $this->x[$i1])
				$i = $i1 - 1;
		}
		if ($i < 0)
			$i = $this->n - 1;
					// 計算
		$xx = $x1 - $this->x[$i];
		$y1 = $this->y[$i] + $xx * ($this->q[$i] + $xx * ($this->r[$i]  + $this->s[$i] * $xx));

		return $y1;
	}
}

					// 設定
	$n = 5;
	$x = array($n+1);
	$y = array($n+1);

	$x[0] = 0.0;
	$x[1] = 1.0;
	$x[2] = 2.0;
	$x[3] = 3.0;
	$x[4] = 4.0;
	$x[5] = 5.0;
	$y[0] = 0.0;
	$y[1] = 1.0;
	$y[2] = 4.0;
	$y[3] = 9.0;
	$y[4] = 16.0;
	$y[5] = 25.0;

	$spline = new Spline($n, $x, $y);
					// 実行と出力
	$x1 = 0.0;
	$sp = 0.1;
	while ($x1 <= 5.01) {
		$y1 = $spline->value($x1);
		printf("%f %f\n", $x1, $y1);
		$x1 += $sp;
	}

?>