<?php /**************************************/ /* シンプレックス法による最小値の計算 */ /* coded by Y.Suganuma */ /**************************************/ $sw = 0; // データの入力 fscanf(STDIN, "%*s %d %*s %d %*s %d %*s %lf", $fun, $n, $max, $k); fscanf(STDIN, "%*s %lf %*s %lf %*s %lf", $eps, $r1, $r2); $x = array($n); $str = trim(fgets(STDIN)); strtok($str, " "); for ($i1 = 0; $i1 < $n; $i1++) $x[$i1] = floatval(strtok(" ")); // 実行 switch ($fun) { case 1: $sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx1"); break; case 2: $sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx2"); break; case 3: $sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx3"); break; } // 結果の出力 if ($sw < 0) printf(" 収束しませんでした!"); else { printf(" 結果="); for ($i1 = 0; $i1 < $n; $i1++) printf("%f ", $x[$i1]); printf(" 最小値=%f 回数=%d\n", $y, $sw); } /******************************************/ /* シンプレックス法 */ /* max : 最大繰り返し回数 */ /* n : 次元 */ /* k : 初期値設定係数 */ /* r1 : 縮小比率 */ /* r2 : 拡大比率 */ /* y : 最小値 */ /* x : x(初期値と答え) */ /* snx : 関数値を計算する関数名 */ /* return : >=0 : 正常終了(収束回数) */ /* =-1 : 収束せず */ /******************************************/ function simplex($max, $n, $k, $eps, $r1, $r2, &$y, &$x, $snx) { // 初期値の設定 $fh = -1; $fg = -1; $fl = -1; $count = 0; $sw = -1; $yy = array($n+1); $xg = array($n); $xr = array($n); $xn = array($n); $xx = array($n+1); for ($i1 = 0; $i1 < $n+1; $i1++) $xx[$i1] = array($n); for ($i1 = 0; $i1 < $n+1; $i1++) { for ($i2 = 0; $i2 < $n; $i2++) $xx[$i1][$i2] = $x[$i2]; if ($i1 > 0) $xx[$i1][$i1-1] += $k; $yy[$i1] = $snx($xx[$i1]); } // 最大値,最小値の計算 for ($i1 = 0; $i1 < $n+1; $i1++) { if ($fh < 0 || $yy[$i1] > $yy[$fh]) $fh = $i1; if ($fl < 0 || $yy[$i1] < $yy[$fl]) $fl = $i1; } for ($i1 = 0; $i1 < $n+1; $i1++) { if ($i1 != $fh && ($fg < 0 || $yy[$i1] > $yy[$fg])) $fg = $i1; } // 最小値の計算 while ($count < $max && $sw < 0) { $count++; // 重心の計算 for ($i1 = 0; $i1 < $n; $i1++) $xg[$i1] = 0.0; for ($i1 = 0; $i1 < $n+1; $i1++) { if ($i1 != $fh) { for ($i2 = 0; $i2 < $n; $i2++) $xg[$i2] += $xx[$i1][$i2]; } } for ($i1 = 0; $i1 < $n; $i1++) $xg[$i1] /= $n; // 最大点の置き換え for ($i1 = 0; $i1 < $n; $i1++) $xr[$i1] = 2.0 * $xg[$i1] - $xx[$fh][$i1]; $yr = $snx($xr); if ($yr >= $yy[$fh]) { // 縮小 for ($i1 = 0; $i1 < $n; $i1++) $xr[$i1] = (1.0 - $r1) * $xx[$fh][$i1] + $r1 * $xr[$i1]; $yr = $snx($xr); } else if ($yr < ($yy[$fl]+($r2-1.0)*$yy[$fh])/$r2) { // 拡大 for ($i1 = 0; $i1 < $n; $i1++) $xn[$i1] = $r2 * $xr[$i1] - ($r2 - 1.0) * $xx[$fh][$i1]; $yn = $snx($xn); if ($yn <= $yr) { for ($i1 = 0; $i1 < $n; $i1++) $xr[$i1] = $xn[$i1]; $yr = $yn; } } for ($i1 = 0; $i1 < $n; $i1++) $xx[$fh][$i1] = $xr[$i1]; $yy[$fh] = $yr; // シンプレックス全体の縮小 if ($yy[$fh] >= $yy[$fg]) { for ($i1 = 0; $i1 < $n+1; $i1++) { if ($i1 != $fl) { for ($i2 = 0; $i2 < $n; $i2++) $xx[$i1][$i2] = 0.5 * ($xx[$i1][$i2] + $xx[$fl][$i2]); $yy[$i1] = $snx($xx[$i1]); } } } // 最大値,最小値の計算 $fh = -1; $fg = -1; $fl = -1; for ($i1 = 0; $i1 < $n+1; $i1++) { if ($fh < 0 || $yy[$i1] > $yy[$fh]) $fh = $i1; if ($fl < 0 || $yy[$i1] < $yy[$fl]) $fl = $i1; } for ($i1 = 0; $i1 < $n+1; $i1++) { if ($i1 != $fh && ($fg < 0 || $yy[$i1] > $yy[$fg])) $fg = $i1; } // 収束判定 $e = 0.0; for ($i1 = 0; $i1 < $n+1; $i1++) { if ($i1 != $fl) { $yr = $yy[$i1] - $yy[$fl]; $e += $yr * $yr; } } if ($e < $eps) $sw = 0; } if ($sw == 0) { $sw = $count; $y = $yy[$fl]; for ($i1 = 0; $i1 < $n; $i1++) $x[$i1] = $xx[$fl][$i1]; } return $sw; } /*******************/ /* 関数値の計算 */ /* y = f(x) */ /* return : y */ /*******************/ // 関数1 function snx1($x) { $x1 = $x[0] - 1.0; $y1 = $x[1] - 2.0; $y = $x1 * $x1 + $y1 * $y1; return $y; } // 関数2 function snx2($x) { $x1 = $x[1] - $x[0] * $x[0]; $y1 = 1.0 - $x[0]; $y = 100.0 * $x1 * $x1 + $y1 * $y1; return $y; } // 関数3 function snx3($x) { $x1 = 1.5 - $x[0] * (1.0 - $x[1]); $y1 = 2.25 - $x[0] * (1.0 - $x[1] * $x[1]); $z1 = 2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]); $y = $x1 * $x1 + $y1 * $y1 + $z1 * $z1; return $y; } ?>