最適化（多項式近似法）

<?php

/*********************************************/
/* 多項式近似によるy=x^4+3x^3+2x^2+1の最小値 */
/*      coded by Y.Suganuma                  */
/*********************************************/

	$x   = -3.0;
	$d   = 0.1;
	$eps = 1.0e-7;
	$max = 100;

	$x = approx($x, $d, $eps, $val, $ind, $max, "snx");

	printf("x %f val %f ind %d\n", $x, $val, $ind);

/****************/
/* 関数値の計算 */
/****************/
function snx($x)
{
	return $x * $x * $x * $x + 3.0 * $x * $x * $x + 2.0 * $x * $x + 1.0;
}

/******************************************/
/* 多項式近似(関数の最小値)               */
/*      x0  : 初期値                      */
/*      d0  : 初期ステップ                */
/*      eps : 許容誤差                    */
/*      val : 関数値                      */
/*      ind : 計算状況                    */
/*              >= 0 : 正常終了(収束回数) */
/*              = -1 : 収束せず           */
/*      max : 最大試行回数                */
/*      fun : 関数値を計算する関数の名前  */
/*      return : 結果                     */
/******************************************/
function approx($x0, $d0, $eps, &$val, &$ind, $max, $fun)
{
	$xx    = 0.0;
	$k     = 0;
	$count = 0;
	$d     = $d0;
	$x     = array(4);
	$x[1]  = $x0;
	$f     = array(4);
	$f[1]  = $fun($x0);
	$ind   = -1;

	while ($count < $max && $ind < 0) {
		$x[3] = $x[1] + $d;
		$f[3] = $fun($x[3]);
		while ($k < $max && $f[3] <= $f[1]) {
			$k++;
			$d *= 2.0;
			$x[0] = $x[1];
			$f[0] = $f[1];
			$x[1] = $x[3];
			$f[1] = $f[3];
			$x[3] = $x[1] + $d;
			$f[3] = $fun($x[3]);
		}
					// 初期値が不適当
		if ($k >= $max)
			$count = $max;
		else {
					// ３点の選択
			$sw = 0;
			if ($k > 0) {
				$x[2] = $x[3] - 0.5 * $d;
				$f[2] = $fun($x[2]);
				$min  = -1;
				for ($i1 = 0; $i1 < 4; $i1++) {
					if ($min < 0 || $f[$i1] < $f[$min])
						$min = $i1;
				}
				if ($min >= 2) {
					for ($i1 = 0; $i1 < 3; $i1++) {
						$x[$i1] = $x[$i1+1];
						$f[$i1] = $f[$i1+1];
					}
				}
				$sw = 1;
			}
			else {
				$x[0] = $x[1] - $d0;
				$f[0] = $fun($x[0]);
				if ($f[0] > $f[1]) {
					$x[2] = $x[3];
					$f[2] = $f[3];
					$sw   = 1;
				}
				else {
					$x[1] = $x[0];
					$f[1] = $f[0];
					$d0   = -$d0;
					$d    = 2.0 * $d0;
					$k    = 1;
				}
			}
					// 収束？
			if ($sw > 0) {
				$count++;
				$dl  = 0.5 * $d * ($f[2] - $f[0]) / ($f[0] - 2.0 * $f[1] + $f[2]);
				$xx  = $x[1] - $dl;
				$val = $fun($xx);
				if (abs($dl) < $eps)
					$ind = $count;
				else {
					$k  = 0;
					$d0 = 0.5 * $d;
					$d  = $d0;
					if ($val < $f[1]) {
						$x[1] = $xx;
						$f[1] = $val;
					}
				}
			}
		}
	}

	return $xx;
}

?>