<?php
/******************************/
/* Newton法による最小値の計算 */
/* coded by Y.Suganuma */
/******************************/
$sw = 0;
// データの入力
fscanf(STDIN, "%*s %d %*s %d %*s %d %*s %d", $fun, $n, $max, $opt_1);
fscanf(STDIN, "%*s %lf %*s %lf", $eps, $step);
$x = array($n);
$dx = array(0.0, 0.0);
$H = array($n);
$str = trim(fgets(STDIN));
strtok($str, " ");
for ($i1 = 0; $i1 < $n; $i1++) {
$x[$i1] = floatval(strtok(" "));
$H[$i1] = array(2*$n);
}
// 実行
switch ($fun) {
case 1:
$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx1", "dsnx1", "hesse1");
break;
case 2:
$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx2", "dsnx2", "hesse2");
break;
case 3:
$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx3", "dsnx3", "hesse3");
break;
}
// 結果の出力
if ($sw < 0) {
printf(" 収束しませんでした！");
switch ($sw) {
case -1:
printf("（収束回数）\n");
break;
case -2:
printf("（１次元最適化の区間）\n");
break;
case -3:
printf("（黄金分割法）\n");
break;
}
}
else {
printf(" 結果＝");
for ($i1 = 0; $i1 < $n; $i1++)
printf("%f ", $x[$i1]);
printf(" 最小値＝%f 回数＝%d\n", $y, $sw);
}
/********************************************************/
/* Newton法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* snx : 関数値を計算する関数名 */
/* dsnx : 関数の微分を計算する関数名(符号を変える) */
/* hesse : Hesse行列の逆行列を計算する関数名 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
function Newton($opt_1, $max, $n, $eps, $step, &$y, &$x, $dx, $H, $snx, $dsnx, $hesse)
{
$wk = array($n);
$count = 0;
$sw = 0;
$y1 = $snx(0.0, $x, $dx);
while ($count < $max && $sw == 0) {
// 傾きの計算
$dsnx($x, $wk);
// Hesse行列の逆行列の計算
$sw1 = $hesse($x, $H, $eps);
// 収束していない
if ($sw1 == 0) {
// 方向の計算
$count++;
for ($i1 = 0; $i1 < $n; $i1++) {
$dx[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++)
$dx[$i1] += $H[$i1][$n+$i2] * $wk[$i2];
}
// １次元最適化を行わない
if ($opt_1 == 0) {
// 新しい点
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] += $dx[$i1];
// 新しい関数値
$y2 = $snx(0.0, $x, $dx);
// 関数値の変化が大きい
if (abs($y2-$y1) > $eps)
$y1 = $y2;
// 収束（関数値の変化＜ｅｐｓ）
else {
$sw = $count;
$y = $y2;
}
}
// １次元最適化を行う
else {
// 区間を決める
$sw1 = 0;
$f1 = $y1;
$sp = $step;
$f2 = $snx($sp, $x, $dx);
if ($f2 > $f1)
$sp = -$step;
for ($i1 = 0; $i1 < $max && $sw1 == 0; $i1++) {
$f2 = $snx($sp, $x, $dx);
if ($f2 > $f1)
$sw1 = 1;
else {
$sp *= 2.0;
$f1 = $f2;
}
}
// 区間が求まらない
if ($sw1 == 0)
$sw = -2;
// 区間が求まった
else {
// 黄金分割法
$k = gold(0.0, $sp, $eps, $y2, $sw1, $max, $x, $dx, $snx);
// 黄金分割法が失敗
if ($sw1 < 0)
$sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] += $k * $dx[$i1];
// 関数値の変化が大きい
if (abs($y1-$y2) > $eps)
$y1 = $y2;
// 収束（関数値の変化＜ｅｐｓ）
else {
$sw = $count;
$y = $y2;
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
$sw = $count;
$y = $y1;
}
}
if ($sw == 0)
$sw = -1;
return $sw;
}
/*******************************/
/* Hesse行列の逆行列を計算する */
/*******************************/
// 関数１
function hesse1($x, &$H, $eps)
{
$H[0][0] = 2.0;
$H[0][1] = 0.0;
$H[1][0] = 0.0;
$H[1][1] = 2.0;
$H[0][2] = 1.0;
$H[0][3] = 0.0;
$H[1][2] = 0.0;
$H[1][3] = 1.0;
$sw = Gauss($H, 2, 2, $eps);
return $sw;
}
// 関数２
function hesse2($x, &$H, $eps)
{
$H[0][0] = 400.0 * (3.0 * $x[0] * $x[0] - $x[1]) + 2.0;
$H[0][1] = -400.0 * $x[0];
$H[1][0] = $H[0][1];
$H[1][1] = 200.0;
$H[0][2] = 1.0;
$H[0][3] = 0.0;
$H[1][2] = 0.0;
$H[1][3] = 1.0;
$sw = Gauss($H, 2, 2, $eps);
return $sw;
}
// 関数３
function hesse3($x, &$H, $eps)
{
$x1 = 1.0 - $x[1];
$x2 = 1.0 - $x[1] * $x[1];
$x3 = 1.0 - $x[1] * $x[1] * $x[1];
$H[0][0] = 2.0 * $x1 * $x1 + 2.0 * $x2 * $x2 + 2.0 * $x3 * $x3;
$H[0][1] = 2.0 * (1.5 - $x[0] * $x1) - 2.0 * $x[0] * $x1 +
4.0 * $x[1] * (2.25 - $x[0] * $x2) - 4.0 * $x[0] * $x[1] * $x2 +
6.0 * $x[1] * $x[1] * (2.625 - $x[0] * $x3) - 6.0 * $x[0] * $x[1] * $x[1] * $x3;
$H[1][0] = $H[0][1];
$H[1][1] = 2.0 * $x[0] * $x[0] +
4.0 * $x[0] * (2.25 - $x[0] * $x2) + 8.0 * $x[0] * $x[0] * $x[1] * $x[1] +
12.0 * $x[0] * $x[1] * (2.625 - $x[0] * $x3) +
18.0 * $x[0] * $x[0] * $x[1] * $x[1] * $x[1] * $x[1];
$H[0][2] = 1.0;
$H[0][3] = 0.0;
$H[1][2] = 0.0;
$H[1][3] = 1.0;
$sw = Gauss($H, 2, 2, $eps);
return $sw;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
function Gauss(&$w, $n, $m, $eps)
{
$ind = 0;
$nm = $n + $m;
for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
$y1 = .0;
$m1 = $i1 + 1;
$m2 = 0;
for ($i2 = $i1; $i2 < $n; $i2++) {
$y2 = abs($w[$i2][$i1]);
if ($y1 < $y2) {
$y1 = $y2;
$m2 = $i2;
}
}
if ($y1 < $eps)
$ind = 1;
else {
for ($i2 = $i1; $i2 < $nm; $i2++) {
$y1 = $w[$i1][$i2];
$w[$i1][$i2] = $w[$m2][$i2];
$w[$m2][$i2] = $y1;
}
$y1 = 1.0 / $w[$i1][$i1];
for ($i2 = $m1; $i2 < $nm; $i2++)
$w[$i1][$i2] *= $y1;
for ($i2 = 0; $i2 < $n; $i2++) {
if ($i2 != $i1) {
for ($i3 = $m1; $i3 < $nm; $i3++)
$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
}
}
}
}
return($ind);
}
/*********************************************/
/* 与えられた点xから，dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する */
/*********************************************/
// 関数１
function snx1($k, $x, $dx)
{
$w = array(2);
$w[0] = $x[0] + $k * $dx[0];
$w[1] = $x[1] + $k * $dx[1];
$x1 = $w[0] - 1.0;
$y1 = $w[1] - 2.0;
$y = $x1 * $x1 + $y1 * $y1;
return $y;
}
// 関数２
function snx2($k, $x, $dx)
{
$w = array(2);
$w[0] = $x[0] + $k * $dx[0];
$w[1] = $x[1] + $k * $dx[1];
$x1 = $w[1] - $w[0] * $w[0];
$y1 = 1.0 - $w[0];
$y = 100.0 * $x1 * $x1 + $y1 * $y1;
return $y;
}
// 関数３
function snx3($k, $x, $dx)
{
$w = array(2);
$w[0] = $x[0] + $k * $dx[0];
$w[1] = $x[1] + $k * $dx[1];
$x1 = 1.5 - $w[0] * (1.0 - $w[1]);
$y1 = 2.25 - $w[0] * (1.0 - $w[1] * $w[1]);
$z1 = 2.625 - $w[0] * (1.0 - $w[1] * $w[1] * $w[1]);
$y = $x1 * $x1 + $y1 * $y1 + $z1 * $z1;
return $y;
}
/********************/
/* 微係数を計算する */
/********************/
// 関数１
function dsnx1($x, &$dx)
{
$dx[0] = -2.0 * ($x[0] - 1.0);
$dx[1] = -2.0 * ($x[1] - 2.0);
}
// 関数２
function dsnx2($x, &$dx)
{
$dx[0] = 400.0 * $x[0] * ($x[1] - $x[0] * $x[0]) + 2.0 * (1.0 - $x[0]);
$dx[1] = -200.0 * ($x[1] - $x[0] * $x[0]);
}
// 関数３
function dsnx3($x, &$dx)
{
$dx[0] = 2.0 * (1.0 - $x[1]) * (1.5 - $x[0] * (1.0 - $x[1])) +
2.0 * (1.0 - $x[1] * $x[1]) * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) +
2.0 * (1.0 - $x[1] * $x[1] * $x[1]) * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
$dx[1] = -2.0 * $x[0] * (1.5 - $x[0] * (1.0 - $x[1])) -
4.0 * $x[0] * $x[1] * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) -
6.0 * $x[0] * $x[1] * $x[1] * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* snx : 関数値を計算する関数の名前 */
/* return : 結果(w＋y＊dwのy) */
/****************************************************************/
function gold($a, $b, $eps, &$val, &$ind, $max, $w, $dw, $snx)
{
// 初期設定
$tau = (sqrt(5.0) - 1.0) / 2.0;
$x = 0.0;
$count = 0;
$ind = -1;
$x1 = $b - $tau * ($b - $a);
$x2 = $a + $tau * ($b - $a);
$f1 = $snx($x1, $w, $dw);
$f2 = $snx($x2, $w, $dw);
// 計算
while ($count < $max && $ind < 0) {
$count += 1;
if ($f2 > $f1) {
if (abs($b-$a) < $eps) {
$ind = 0;
$x = $x1;
$val = $f1;
}
else {
$b = $x2;
$x2 = $x1;
$x1 = $a + (1.0 - $tau) * ($b - $a);
$f2 = $f1;
$f1 = $snx($x1, $w, $dw);
}
}
else {
if (abs($b-$a) < $eps) {
$ind = 0;
$x = $x2;
$val = $f2;
$f1 = $f2;
}
else {
$a = $x1;
$x1 = $x2;
$x2 = $b - (1.0 - $tau) * ($b - $a);
$f1 = $f2;
$f2 = $snx($x2, $w, $dw);
}
}
}
// 収束した場合の処理
if ($ind == 0) {
$ind = $count;
$fa = $snx($a, $w, $dw);
$fb = $snx($b, $w, $dw);
if ($fb < $fa) {
$a = $b;
$fa = $fb;
}
if ($fa < $f1) {
$x = $a;
$val = $fa;
}
}
return $x;
}
?>