<?php /*********************************************/ /* 黄金分割法によるy=x^4+3x^3+2x^2+1の最小値 */ /* coded by Y.Suganuma */ /*********************************************/ $a = -2.0; $b = -1.0; $eps = 1.0e-10; $max = 100; $x = gold($a, $b, $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; } /******************************************/ /* 黄金分割法(関数の最小値) */ /* a,b : 初期区間 a < b */ /* eps : 許容誤差 */ /* val : 間数値 */ /* ind : 計算状況 */ /* >= 0 : 正常終了(収束回数) */ /* = -1 : 収束せず */ /* max : 最大試行回数 */ /* fun : 関数値を計算する関数の名前 */ /* return : 結果 */ /******************************************/ function gold($a, $b, $eps, &$val, &$ind, $max, $fun) { $x = 0.0; $tau = (sqrt(5.0) - 1.0) / 2.0; $count = 0; $ind = -1; $x1 = $b - $tau * ($b - $a); $x2 = $a + $tau * ($b - $a); $f1 = $fun($x1); $f2 = $fun($x2); while ($count < $max && $ind < 0) { $count += 1; if ($f2 > $f1) { if (abs($b-$a) < $eps && abs($b-$a) < $eps*abs($b)) { $ind = 0; $x = $x1; $val = $f1; } else { $b = $x2; $x2 = $x1; $x1 = $a + (1.0 - $tau) * ($b - $a); $f2 = $f1; $f1 = $fun($x1); } } else { if (abs($b-$a) < $eps && abs($b-$a) < $eps*abs($b)) { $ind = 0; $x = $x2; $val = $f2; $f1 = $f2; } else { $a = $x1; $x1 = $x2; $x2 = $b - (1.0 - $tau) * ($b - $a); $f1 = $f2; $f2 = $fun($x2); } } } if ($ind == 0) { $ind = $count; $fa = $fun($a); $fb = $fun($b); if ($fb < $fa) { $a = $b; $fa = $fb; } if ($fa < $f1) { $x = $a; $val = $fa; } } return $x; } ?>