| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/********************************/
/* 共役勾配法による最小値の計算 */
/* coded by Y.Suganuma */
/********************************/
#include <stdio.h>
void dsnx1(double *, double *);
double snx1(double, double *, double *);
void dsnx2(double *, double *);
double snx2(double, double *, double *);
void dsnx3(double *, double *);
double snx3(double, double *, double *);
int conjugate(int, int, int, double, double, double *, double *, double *,
double (*)(double, double *, double *),
void (*)(double *, double *));
int main()
{
double eps, step, *x, *dx, y;
int fun, i1, max, n, opt_1, sw = 0;
// データの入力
scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
scanf("%*s %lf %*s %lf", &eps, &step);
x = new double [n];
dx = new double [n];
scanf("%*s");
for (i1 = 0; i1 < n; i1++)
scanf("%lf", &x[i1]);
// 実行
switch (fun) {
case 1:
sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, snx1, dsnx1);
break;
case 2:
sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, snx2, dsnx2);
break;
case 3:
sw = conjugate(opt_1, max, n, eps, step, &y, x, dx, snx3, dsnx3);
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);
}
return 0;
}
/********************************************************/
/* 共役勾配法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* snx : 関数値を計算する関数名 */
/* dsnx : 関数の微分を計算する関数名(符号を変える) */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
#include <math.h>
double gold(double, double, double, double *, int *, int, double *, double *,
double (*)(double, double *, double *));
int conjugate(int opt_1, int max, int n, double eps, double step, double *y,
double *x, double *dx, double (*snx)(double, double *, double *),
void (*dsnx)(double *, double *))
{
double b = 0.0, f1, f2, *g, k, s1 = 1.0, s2, sp, y1, y2;
int count = 0, i1, sw = 0, sw1;
g = new double [n];
y1 = snx(0.0, x, dx);
while (count < max && sw == 0) {
// 傾きの計算
dsnx(x, g);
// 傾きのチェック
s2 = 0.0;
for (i1 = 0; i1 < n; i1++)
s2 += g[i1] * g[i1];
// 収束していない
if (sqrt(s2) > eps) {
// 方向の計算
count++;
if (count%n == 1)
b = 0.0;
else
b = s2 / s1;
for (i1 = 0; i1 < n; i1++)
dx[i1] = g[i1] + b * dx[i1];
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (i1 = 0; i1 < n; i1++)
x[i1] += step * dx[i1];
// 新しい関数値
y2 = snx(0.0, x, dx);
// 関数値の変化が大きい
if (fabs(y2-y1) > eps) {
y1 = y2;
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
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 (fabs(y1-y2) > eps) {
y1 = y2;
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
*y = y2;
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
sw = count;
*y = y1;
}
}
if (sw == 0)
sw = -1;
delete [] g;
return sw;
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* snx : 関数値を計算する関数の名前 */
/* return : 結果(w＋y＊dwのy) */
/****************************************************************/
#include <math.h>
double gold(double a, double b, double eps, double *val, int *ind, int max,
double *w, double *dw, double (*snx)(double, double *, double *))
{
double f1, f2, fa, fb, tau = (sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
int 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 (fabs(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 (fabs(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;
}
/*********************************************/
/* 与えられた点xから，dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する */
/*********************************************/
// 関数１
double snx1(double k, double *x, double *dx)
{
double x1, y1, y, w[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;
}
// 関数２
double snx2(double k, double *x, double *dx)
{
double x1, y1, y, w[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;
}
// 関数３
double snx3(double k, double *x, double *dx)
{
double x1, y1, z1, y, w[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;
}
/********************/
/* 微係数を計算する */
/********************/
// 関数１
void dsnx1(double *x, double *dx)
{
dx[0] = -2.0 * (x[0] - 1.0);
dx[1] = -2.0 * (x[1] - 2.0);
}
// 関数２
void dsnx2(double *x, double *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]);
}
// 関数３
void dsnx3(double *x, double *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]));
}
/********************************/
/* 共役勾配法による最小値の計算 */
/* coded by Y.Suganuma */
/********************************/
import java.io.*;
import java.util.*;
public class Test {
public static void main(String args[]) throws IOException
{
double eps, step, x[], dx[], y[];
int fun, i1, max, n, opt_1, sw = 0;
StringTokenizer str;
String line;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
// データの入力
// 1 行目
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
line = str.nextToken();
fun = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
n = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
max = Integer.parseInt(line);
line = str.nextToken();
line = str.nextToken();
opt_1 = Integer.parseInt(line);
// 2 行目
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
line = str.nextToken();
eps = Double.parseDouble(line);
line = str.nextToken();
line = str.nextToken();
step = Double.parseDouble(line);
// 3 行目
x = new double [n];
dx = new double [n];
y = new double [1];
line = in.readLine();
str = new StringTokenizer(line, " ");
line = str.nextToken();
for (i1 = 0; i1 < n; i1++) {
line = str.nextToken();
x[i1] = Double.parseDouble(line);
}
// 実行
Kansu kn = new Kansu(fun);
sw = kn.conjugate(opt_1, max, n, eps, step, y, x, dx);
// 結果の出力
if (sw < 0) {
System.out.print(" 収束しませんでした！");
switch (sw) {
case -1:
System.out.println("（収束回数）");
break;
case -2:
System.out.print("（１次元最適化の区間）");
break;
case -3:
System.out.print("（黄金分割法）");
break;
}
}
else {
System.out.print(" 結果＝");
for (i1 = 0; i1 < n; i1++)
System.out.print(x[i1] + " ");
System.out.println(" 最小値＝" + y[0] + " 回数＝" + sw);
}
}
}
/************************/
/* 関数値，微係数の計算 */
/************************/
class Kansu extends Conjugate
{
private int sw;
// コンストラクタ
Kansu (int s) {sw = s;}
// 与えられた点xから，dx方向へk*dxだけ進んだ
// 点における関数値を計算する
double snx(double k, double x[], double dx[])
{
double x1, y1, z1, y = 0.0, w[] = new double [2];
switch (sw) {
case 1:
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;
break;
case 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;
break;
case 3:
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;
break;
}
return y;
}
// 関数の微分
void dsnx(double x[], double dx[])
{
switch (sw) {
case 1:
dx[0] = -2.0 * (x[0] - 1.0);
dx[1] = -2.0 * (x[1] - 2.0);
break;
case 2:
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]);
break;
case 3:
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]));
break;
}
}
}
abstract class Conjugate
{
/********************************************************/
/* 共役勾配法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
abstract double snx(double k, double x[], double dx[]);
abstract void dsnx(double x[], double dx[]);
int conjugate(int opt_1, int max, int n, double eps, double step, double y[], double x[], double dx[])
{
double b = 0.0, f1, f2, g[], k, s1 = 1.0, s2, sp,
y1[] = new double [1], y2[] = new double [1];
int count = 0, i1, sw = 0, sw1[] = new int [1];
g = new double [n];
y1[0] = snx(0.0, x, dx);
while (count < max && sw == 0) {
// 傾きの計算
dsnx(x, g);
// 傾きのチェック
s2 = 0.0;
for (i1 = 0; i1 < n; i1++)
s2 += g[i1] * g[i1];
// 収束していない
if (Math.sqrt(s2) > eps) {
// 方向の計算
count++;
if (count%n == 1)
b = 0.0;
else
b = s2 / s1;
for (i1 = 0; i1 < n; i1++)
dx[i1] = g[i1] + b * dx[i1];
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (i1 = 0; i1 < n; i1++)
x[i1] += step * dx[i1];
// 新しい関数値
y2[0] = snx(0.0, x, dx);
// 関数値の変化が大きい
if (Math.abs(y2[0]-y1[0]) > eps) {
y1[0] = y2[0];
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y[0] = y2[0];
}
}
// １次元最適化を行う
else {
// 区間を決める
sw1[0] = 0;
f1 = y1[0];
sp = step;
f2 = snx(sp, x, dx);
if (f2 > f1)
sp = -step;
for (i1 = 0; i1 < max && sw1[0] == 0; i1++) {
f2 = snx(sp, x, dx);
if (f2 > f1)
sw1[0] = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1[0] == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x, dx);
// 黄金分割法が失敗
if (sw1[0] < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (i1 = 0; i1 < n; i1++)
x[i1] += k * dx[i1];
// 関数値の変化が大きい
if (Math.abs(y1[0]-y2[0]) > eps) {
y1[0] = y2[0];
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y[0] = y2[0];
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
sw = count;
y[0] = y1[0];
}
}
if (sw == 0)
sw = -1;
return sw;
}
/****************************************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* return : 結果(w＋y＊dwのy) */
/****************************************************************/
double gold(double a, double b, double eps, double val[], int ind[], int max, double w[], double dw[])
{
double f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
int count = 0;
// 初期設定
ind[0] = -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] < 0) {
count += 1;
if (f2 > f1) {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x1;
val[0] = f1;
}
else {
b = x2;
x2 = x1;
x1 = a + (1.0 - tau) * (b - a);
f2 = f1;
f1 = snx(x1, w, dw);
}
}
else {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x2;
val[0] = 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] == 0) {
ind[0] = count;
fa = snx(a, w, dw);
fb = snx(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val[0] = fa;
}
}
return x;
}
}
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>最適化（共役勾配法）</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
str1 = "";
str2 = new Array();
function main()
{
// データの設定
let n = parseInt(document.getElementById("n").value);
let max = parseInt(document.getElementById("max").value);
let step = parseFloat(document.getElementById("step").value);
let x = document.getElementById("x0").value.split(/ {1,}/)
for (let i1 = 0; i1 < n; i1++)
x[i1] = parseFloat(x[i1]);
let opt_1 = 1;
if (document.getElementById("line").options[1].selected)
opt_1 = 0;
str1 = document.getElementById("func").value + ";";
str2 = document.getElementById("dfunc").value.split("\n")
for (let i1 = 0; i1 < n; i1++)
str2[i1] = str2[i1] + ";";
let eps = 1.0e-10;
let y = new Array();
let dx = new Array();
for (let i1 = 0; i1 < n; i1++)
dx[i1] = 0.0;
// 実行
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, snx, dsnx);
// 結果
let res;
if (sw < 0) {
res = " 収束しませんでした！";
switch (sw) {
case -1:
res += "（収束回数）";
break;
case -2:
res += "（１次元最適化の区間）";
break;
case -3:
res += "（黄金分割法）";
break;
}
document.getElementById("res").value = res;
}
else {
let c = 100000;
res = " x ＝";
for (let i1 = 0; i1 < n; i1++) {
let s1 = Math.round(x[i1] * c) / c;
res = res + " " + s1;
}
res += "\n";
let s2 = Math.round(y[0] * c) / c;
res = res + " 最小値 ＝ " + s2 + " 回数 ＝ " + sw;
document.getElementById("res").value = res;
}
}
/****************/
/* 関数値の計算 */
/****************/
function snx(k, x, dx)
{
let y = eval(str1);
return y;
}
/************************/
/* 関数値の微分値の計算 */
/************************/
function dsnx(n, x, dx)
{
for (let i1 = 0; i1 < n; i1++)
dx[i1] = eval(str2[i1]);
}
/******************************************************/
/* 共役傾斜法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* fn : 関数値を計算する関数 */
/* dfn : 関数の微分値を計算する関数 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/******************************************************/
function conjugate(opt_1, max, n, eps, step, y, x, dx, fn, dfn)
{
let b = 0.0, f1, f2, k, s1 = 1.0, s2, sp, count = 0, sw = 0;
let y1 = new Array();
let y2 = new Array();
let sw1 = new Array(9);
let g = new Array();
y1[0] = fn(0.0, x, dx);
while (count < max && sw == 0) {
// 傾きの計算
dfn(n, x, g);
// 傾きのチェック
s2 = 0.0;
for (let i1 = 0; i1 < n; i1++)
s2 += g[i1] * g[i1];
// 収束していない
if (Math.sqrt(s2) > eps) {
// 方向の計算
count++;
if (count%n == 1)
b = 0.0;
else
b = s2 / s1;
for (let i1 = 0; i1 < n; i1++)
dx[i1] = g[i1] + b * dx[i1];
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (let i1 = 0; i1 < n; i1++)
x[i1] += step * dx[i1];
// 新しい関数値
y2[0] = fn(0.0, x, dx);
// 関数値の変化が大きい
if (Math.abs(y2[0]-y1[0]) > eps) {
y1[0] = y2[0];
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y[0] = y2[0];
}
}
// １次元最適化を行う
else {
// 区間を決める
sw1[0] = 0;
f1 = y1[0];
sp = step;
f2 = fn(sp, x, dx);
if (f2 > f1)
sp = -step;
for (let i1 = 0; i1 < max && sw1[0] == 0; i1++) {
f2 = fn(sp, x, dx);
if (f2 > f1)
sw1[0] = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1[0] == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x, dx, fn);
// 黄金分割法が失敗
if (sw1[0] < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (let i1 = 0; i1 < n; i1++)
x[i1] += k * dx[i1];
// 関数値の変化が大きい
if (Math.abs(y1[0]-y2[0]) > eps) {
y1[0] = y2[0];
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y[0] = y2[0];
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
sw = count;
y[0] = y1[0];
}
}
if (sw == 0)
sw = -1;
return sw;
}
/******************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* fn : 関数値を計算する関数 */
/* return : 結果(w＋y＊dwのy) */
/******************************************/
function gold(a, b, eps, val, ind, max, w, dw, fn)
{
let f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
let count = 0;
// 初期設定
ind[0] = -1;
x1 = b - tau * (b - a);
x2 = a + tau * (b - a);
f1 = fn(x1, w, dw);
f2 = fn(x2, w, dw);
// 計算
while (count < max && ind[0] < 0) {
count += 1;
if (f2 > f1) {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x1;
val[0] = f1;
}
else {
b = x2;
x2 = x1;
x1 = a + (1.0 - tau) * (b - a);
f2 = f1;
f1 = fn(x1, w, dw);
}
}
else {
if (Math.abs(b-a) < eps) {
ind[0] = 0;
x = x2;
val[0] = f2;
f1 = f2;
}
else {
a = x1;
x1 = x2;
x2 = b - (1.0 - tau) * (b - a);
f1 = f2;
f2 = fn(x2, w, dw);
}
}
}
// 収束した場合の処理
if (ind[0] == 0) {
ind[0] = count;
fa = fn(a, w, dw);
fb = fn(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val[0] = fa;
}
}
return x;
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>最適化（共役勾配法）</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，以下に示すような関数の最小値を求める場合に対する値（一次元最適化を行わない場合は，刻み幅を 0.003 にしてみて下さい）が設定されています（他の問題を実行する場合は，それらを適切に修正してください）．
<P STYLE="text-align:center"><IMG SRC="conjugate.gif"></P>
<DT>n 変数関数 f(<B>x</B>) に対する共役勾配法においては，現在の点 <B>x</B> における傾き，
<P STYLE="text-align:center"><IMG SRC="conjugate1.gif"></P>
<DT>の結果から共役勾配の方向 <B>p</B> を求め、その方向にαだけ進んだ点の関数値 f( <B>x</B> + α<B>p</B>) を計算する，という処理を繰り返します．目的関数の箇所には，この例に対する f( <B>x</B> + α<B>p</B>) の計算式，また，目的関数の微分の箇所には，<B>g</B> の計算式が与えられています．ただし，プログラム内では，変数 x は，この例の x と y からなる配列，変数 dx は <B>g</B> に対応し，関数 f を x 及び y で微分したものからなる配列，また，k はαに対応する変数として表現されています．なお，目的関数及び目的関数の微分は，JavaScript の仕様に適合した形式で記述してあることに注意してください．
</DL>
<DIV STYLE="text-align:center">
変数の数(n)：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="2">
最大繰り返し回数：<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="5000">
刻み幅：<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.001"><BR>
初期値(n個)：<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="20" VALUE="0.0 0.0">
１次元最適化：<SELECT STYLE="font-size: 100%" SIZE="1" ID="line">
<OPTION VALUE="yes"> 行う </OPTION><BR>
<OPTION VALUE="no"> 行わない </OPTION><BR>
</SELECT><BR><BR>
目的関数： f(x+αp) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="70" VALUE="100.0 * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) + (1.0 - (x[0] + k * dx[0])) * (1.0 - (x[0] + k * dx[0]))"><BR><BR>
目的関数の微分： f'(x) = <TEXTAREA ID="dfunc" STYLE="font-size: 110%" COLS="70" ROWS="10">
400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
-200.0 * (x[1] - x[0] * x[0])
</TEXTAREA><BR><BR>
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
結果：<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="40" ROWS="2"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/********************************/
/* 共役勾配法による最小値の計算 */
/* 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);
$str = trim(fgets(STDIN));
strtok($str, " ");
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] = floatval(strtok(" "));
// 実行
switch ($fun) {
case 1:
$sw = conjugate($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx1", "dsnx1");
break;
case 2:
$sw = conjugate($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx2", "dsnx2");
break;
case 3:
$sw = conjugate($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx3", "dsnx3");
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);
}
/********************************************************/
/* 共役勾配法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* snx : 関数値を計算する関数名 */
/* dsnx : 関数の微分を計算する関数名(符号を変える) */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
function conjugate($opt_1, $max, $n, $eps, $step, &$y, &$x, $dx, $snx, $dsnx)
{
$b = 0.0;
$s1 = 1.0;
$count = 0;
$sw = 0;
$g = array($n);
$y1 = $snx(0.0, $x, $dx);
while ($count < $max && $sw == 0) {
// 傾きの計算
$dsnx($x, $g);
// 傾きのチェック
$s2 = 0.0;
for ($i1 = 0; $i1 < $n; $i1++)
$s2 += $g[$i1] * $g[$i1];
// 収束していない
if (sqrt($s2) > $eps) {
// 方向の計算
$count++;
if ($count%$n == 1)
$b = 0.0;
else
$b = $s2 / $s1;
for ($i1 = 0; $i1 < $n; $i1++)
$dx[$i1] = $g[$i1] + $b * $dx[$i1];
// １次元最適化を行わない
if ($opt_1 == 0) {
// 新しい点
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] += $step * $dx[$i1];
// 新しい関数値
$y2 = $snx(0.0, $x, $dx);
// 関数値の変化が大きい
if (abs($y2-$y1) > $eps) {
$y1 = $y2;
$s1 = $s2;
}
// 収束（関数値の変化＜ｅｐｓ）
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;
$s1 = $s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
$sw = $count;
$y = $y2;
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
$sw = $count;
$y = $y1;
}
}
if ($sw == 0)
$sw = -1;
return $sw;
}
/*********************************************/
/* 与えられた点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;
}
?>
#*******************************/
# 共役勾配法による最小値の計算 */
# coded by Y.Suganuma */
#*******************************/
#********************************************/
# 与えられた点xから，dx方向へk*dxだけ進んだ */
# 点における関数値及び微係数を計算する */
#********************************************/
# 関数１
snx1 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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
x1 * x1 + y1 * y1
else
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
end
}
# 関数２
snx2 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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]
100.0 * x1 * x1 + y1 * y1
else
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])
end
}
# 関数３
snx3 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(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])
x1 * x1 + y1 * y1 + z1 * z1
else
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]))
end
}
#***************************************************************/
# 黄金分割法(与えられた方向での最小値) */
# a,b : 初期区間 a < b */
# eps : 許容誤差 */
# val : 関数値 */
# ind : 計算状況 */
# >= 0 : 正常終了(収束回数) */
# = -1 : 収束せず */
# max : 最大試行回数 */
# w : 位置 */
# dw : 傾きの成分 */
# snx : 関数値を計算する関数の名前 */
# return : 結果(w＋y＊dwのy) */
#***************************************************************/
def gold(a, b, eps, val, ind, max, w, dw, &snx)
# 初期設定
tau = (Math.sqrt(5.0) - 1.0) / 2.0
x = 0.0
count = 0
ind[0] = -1
x1 = b - tau * (b - a)
x2 = a + tau * (b - a)
f1 = snx.call(0, x1, w, dw)
f2 = snx.call(0, x2, w, dw)
# 計算
while count < max && ind[0] < 0
count += 1
if f2 > f1
if (b-a).abs() < eps
ind[0] = 0
x = x1
val[0] = f1
else
b = x2
x2 = x1
x1 = a + (1.0 - tau) * (b - a)
f2 = f1
f1 = snx.call(0, x1, w, dw)
end
else
if (b-a).abs() < eps
ind[0] = 0
x = x2
val[0] = f2
f1 = f2
else
a = x1
x1 = x2
x2 = b - (1.0 - tau) * (b - a)
f1 = f2
f2 = snx.call(0, x2, w, dw)
end
end
end
# 収束した場合の処理
if ind[0] == 0
ind[0] = count
fa = snx.call(0, a, w, dw)
fb = snx.call(0, b, w, dw)
if fb < fa
a = b
fa = fb
end
if fa < f1
x = a
val[0] = fa
end
end
return x
end
#*******************************************************/
# 共役勾配法 */
# opt_1 : =0 : １次元最適化を行わない */
# =1 : １次元最適化を行う */
# max : 最大繰り返し回数 */
# n : 次元 */
# eps : 収束判定条件 */
# step : きざみ幅 */
# y : 最小値 */
# x : x(初期値と答え) */
# dx : 関数の微分値 */
# snx : 関数値とその微分を計算する関数名 */
# （微分は，その符号を変える) */
# return : >=0 : 正常終了(収束回数) */
# =-1 : 収束せず */
# =-2 : １次元最適化の区間が求まらない */
# =-3 : 黄金分割法が失敗 */
#*******************************************************/
def conjugate(opt_1, max, n, eps, step, y, x, dx, &snx)
b = 0.0
s1 = 1.0
count = 0
sw = 0
y2 = Array.new(1)
sw1 = Array.new(1)
g = Array.new(n)
y1 = snx.call(0, 0.0, x, dx)
while count < max && sw == 0
# 傾きの計算
snx.call(1, 0.0, x, g)
# 傾きのチェック
s2 = 0.0
for i1 in 0 ... n
s2 += g[i1] * g[i1]
end
# 収束していない
if Math.sqrt(s2) > eps
# 方向の計算
count += 1
if count%n == 1
b = 0.0
else
b = s2 / s1
end
for i1 in 0 ... n
dx[i1] = g[i1] + b * dx[i1]
end
# １次元最適化を行わない
if opt_1 == 0
# 新しい点
for i1 in 0 ... n
x[i1] += step * dx[i1]
end
# 新しい関数値
y2[0] = snx.call(0, 0.0, x, dx)
# 関数値の変化が大きい
if (y2[0]-y1).abs() > eps
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else
sw = count
y[0] = y2[0]
end
# １次元最適化を行う
else
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx.call(0, sp, x, dx)
if f2 > f1
sp = -step
end
for i1 in 0 ... max
f2 = snx.call(0, sp, x, dx)
if f2 > f1
sw1[0] = 1
break
else
sp *= 2.0
f1 = f2
end
end
# 区間が求まらない
if sw1[0] == 0
sw = -2
# 区間が求まった
else
# 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x, dx, &snx)
# 黄金分割法が失敗
if sw1[0] < 0
sw = -3
# 黄金分割法が成功
else
# 新しい点
for i1 in 0 ... n
x[i1] += k * dx[i1]
end
# 関数値の変化が大きい
if (y1-y2[0]).abs() > eps
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else
sw = count
y[0] = y2[0]
end
end
end
end
# 収束（傾き＜ｅｐｓ）
else
sw = count
y[0] = y1
end
end
if sw == 0
sw = -1
end
return sw
end
# データの入力
str = gets()
a = str.split(" ")
fun = Integer(a[1])
n = Integer(a[3])
max = Integer(a[5])
opt_1 = Integer(a[7])
str = gets()
a = str.split(" ")
eps = Float(a[1])
step = Float(a[3])
x = Array.new(n)
dx = Array.new(n)
y = Array.new(1)
sw = 0
str = gets()
a = str.split(" ")
for i1 in 0 ... n
x[i1] = Float(a[i1+1])
dx[i1] = 0.0
end
# 実行
case fun
when 1
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &snx1)
when 2
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &snx2)
when 3
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &snx3)
end
# 結果の出力
if sw < 0
printf(" 収束しませんでした！")
case sw
when -1
printf("（収束回数）\n")
when -2
printf("（１次元最適化の区間）\n")
when -3
printf("（黄金分割法）\n")
end
else
printf(" 結果＝")
for i1 in 0 ... n
printf("%f ", x[i1])
end
printf(" 最小値＝%f 回数＝%d\n", y[0], sw)
end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
############################################
# 黄金分割法(与えられた方向での最小値)
# a,b : 初期区間 a < b
# eps : 許容誤差
# val : 関数値
# ind : 計算状況
# >= 0 : 正常終了(収束回数)
# = -1 : 収束せず
# max : 最大試行回数
# w : 位置
# dw : 傾きの成分
# snx : 関数値を計算する関数の名前
# return : 結果(w＋y＊dwのy)
# coded by Y.Suganuma
############################################
def gold(a, b, eps, val, ind, max, w, dw, snx) :
# 初期設定
tau = (sqrt(5.0) - 1.0) / 2.0
x = 0.0
count = 0
ind[0] = -1
x1 = b - tau * (b - a)
x2 = a + tau * (b - a)
f1 = snx(x1, w, dw)
f2 = snx(x2, w, dw)
# 計算
while count < max and ind[0] < 0 :
count += 1
if f2 > f1 :
if abs(b-a) < eps :
ind[0] = 0
x = x1
val[0] = 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] = 0
x = x2
val[0] = 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] == 0 :
ind[0] = count
fa = snx(a, w, dw)
fb = snx(b, w, dw)
if fb < fa :
a = b
fa = fb
if fa < f1 :
x = a
val[0] = fa
return x
############################################
# 共役勾配法
# opt_1 : =0 : １次元最適化を行わない
# =1 : １次元最適化を行う
# max : 最大繰り返し回数
# n : 次元
# eps : 収束判定条件
# step : きざみ幅
# y : 最小値
# x : x(初期値と答え)
# dx : 関数の微分値
# snx : 関数値を計算する関数名
# dsnx : 関数の微分を計算する関数名(符号を変える)
# return : >=0 : 正常終了(収束回数)
# =-1 : 収束せず
# =-2 : １次元最適化の区間が求まらない
# =-3 : 黄金分割法が失敗
############################################
def conjugate(opt_1, max, n, eps, step, y, x, dx, snx, dsnx) :
b = 0.0
s1 = 1.0
count = 0
sw = 0
g = np.empty(n, np.float)
y1 = snx(0.0, x, dx)
y2 = np.empty(1, np.float)
sw1 = np.empty(1, np.int)
while count < max and sw == 0 :
# 傾きの計算
dsnx(x, g)
# 傾きのチェック
s2 = 0.0
for i1 in range(0, n) :
s2 += g[i1] * g[i1]
# 収束していない
if sqrt(s2) > eps :
# 方向の計算
count += 1
if count%n == 1 :
b = 0.0
else :
b = s2 / s1
for i1 in range(0, n) :
dx[i1] = g[i1] + b * dx[i1]
# １次元最適化を行わない
if opt_1 == 0 :
# 新しい点
for i1 in range(0, n) :
x[i1] += step * dx[i1]
# 新しい関数値
y2[0] = snx(0.0, x, dx)
# 関数値の変化が大きい
if abs(y2[0]-y1) > eps :
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else :
sw = count
y[0] = y2[0]
# １次元最適化を行う
else :
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx(sp, x, dx)
if f2 > f1 :
sp = -step
for i1 in range(0, max) :
f2 = snx(sp, x, dx)
if f2 > f1 :
sw1[0] = 1
break
else :
sp *= 2.0
f1 = f2
# 区間が求まらない
if sw1[0] == 0 :
sw = -2
# 区間が求まった
else :
# 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x, dx, snx)
# 黄金分割法が失敗
if sw1[0] < 0 :
sw = -3
# 黄金分割法が成功
else :
# 新しい点
for i1 in range(0, n) :
x[i1] += k * dx[i1]
# 関数値の変化が大きい
if abs(y1-y2[0]) > eps :
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else :
sw = count
y[0] = y2[0]
# 収束（傾き＜ｅｐｓ）
else :
sw = count
y[0] = y1
if sw == 0 :
sw = -1
return sw
############################################
# 共役勾配法による最小値の計算
# coded by Y.Suganuma
############################################
############################################
# 与えられた点xから，dx方向へk*dxだけ進んだ
# 点における関数値を計算する
############################################
# 関数１
def snx1(k, x, dx) :
w = np.empty(2, np.float)
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
# 関数２
def snx2(k, x, dx) :
w = np.empty(2, np.float)
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
# 関数３
def snx3(k, x, dx) :
w = np.empty(2, np.float)
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
############################################
# 微係数を計算する
############################################
# 関数１
def dsnx1(x, dx) :
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
# 関数２
def 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])
# 関数３
def 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]))
# データの入力
sw = 0
line = sys.stdin.readline()
s = line.split()
fun = int(s[1])
n = int(s[3])
max = int(s[5])
opt_1 = int(s[7])
line = sys.stdin.readline()
s = line.split()
eps = float(s[1])
step = float(s[3])
x = np.empty(n, np.float)
dx = np.empty(n, np.float)
y = np.empty(1, np.float)
line = sys.stdin.readline()
s = line.split()
for i1 in range(0, n) :
x[i1] = float(s[i1+1])
# 実行
if fun == 1 :
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, snx1, dsnx1)
elif fun == 2 :
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, snx2, dsnx2)
else :
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, snx3, dsnx3)
# 結果の出力
if sw < 0 :
print(" 収束しませんでした！", end="")
if sw == -1 :
print("（収束回数）")
elif sw == -2 :
print("（１次元最適化の区間）")
else :
print("（黄金分割法）")
else :
print(" 結果＝", end="")
for i1 in range(0, n) :
print(str(x[i1]) + " ", end="")
print(" 最小値＝" + str(y[0]) + " 回数＝" + str(sw))
/********************************/
/* 共役勾配法による最小値の計算 */
/* coded by Y.Suganuma */
/********************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
public Test1()
{
// データの入力
// 1 行目
string[] str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
int fun = int.Parse(str[1]);
int n = int.Parse(str[3]);
int max = int.Parse(str[5]);
int opt_1 = int.Parse(str[7]);
// 2 行目
str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
double eps = double.Parse(str[1]);
double step = double.Parse(str[3]);
// 3 行目
double[] x = new double [n];
double[] dx = new double [n];
double y = 0.0;
str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
for (int i1 = 0; i1 < n; i1++)
x[i1] = double.Parse(str[i1+1]);
// 実行
int sw;
if (fun == 1)
sw = conjugate(opt_1, max, n, eps, step, ref y, x, dx, snx1, dsnx1);
else if (fun == 2)
sw = conjugate(opt_1, max, n, eps, step, ref y, x, dx, snx2, dsnx2);
else
sw = conjugate(opt_1, max, n, eps, step, ref y, x, dx, snx3, dsnx3);
// 結果の出力
if (sw < 0) {
Console.Write(" 収束しませんでした！");
switch (sw) {
case -1:
Console.WriteLine("（収束回数）");
break;
case -2:
Console.WriteLine("（１次元最適化の区間）");
break;
case -3:
Console.WriteLine("（黄金分割法）");
break;
}
}
else {
Console.Write(" 結果＝");
for (int i1 = 0; i1 < n; i1++)
Console.Write(x[i1] + " ");
Console.WriteLine(" 最小値＝" + y + " 回数＝" + sw);
}
}
// 関数１の値の計算
double snx1(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = w[0] - 1.0;
double y1 = w[1] - 2.0;
double y = x1 * x1 + y1 * y1;
return y;
}
// 関数２の値の計算
double snx2(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = w[1] - w[0] * w[0];
double y1 = 1.0 - w[0];
double y = 100.0 * x1 * x1 + y1 * y1;
return y;
}
// 関数３の値の計算
double snx3(double k, double[] x, double[] dx)
{
double[] w = new double [2];
w[0] = x[0] + k * dx[0];
w[1] = x[1] + k * dx[1];
double x1 = 1.5 - w[0] * (1.0 - w[1]);
double y1 = 2.25 - w[0] * (1.0 - w[1] * w[1]);
double z1 = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
double y = x1 * x1 + y1 * y1 + z1 * z1;
return y;
}
// 関数１の微分
int dsnx1(double[] x, double[] dx)
{
dx[0] = -2.0 * (x[0] - 1.0);
dx[1] = -2.0 * (x[1] - 2.0);
return 0;
}
// 関数２の微分
int dsnx2(double[] x, double[] 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]);
return 0;
}
// 関数３の微分
int dsnx3(double[] x, double[] 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]));
return 0;
}
/********************************************************/
/* 共役勾配法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* fn : 関数値を計算する関数 */
/* dfn : 関数の微分値を計算する関数 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
int conjugate(int opt_1, int max, int n, double eps, double step, ref double y,
double[] x, double[] dx, Func<double, double[], double[], double> fn,
Func<double[], double[], int> dfn)
{
double[] g = new double [n];
double s1 = 1.0;
double s2;
double y1 = fn(0.0, x, dx);
double y2 = 0.0;
int count = 0;
int sw = 0;
while (count < max && sw == 0) {
// 傾きの計算
dfn(x, g);
// 傾きのチェック
s2 = 0.0;
for (int i1 = 0; i1 < n; i1++)
s2 += g[i1] * g[i1];
// 収束していない
if (Math.Sqrt(s2) > eps) {
// 方向の計算
count++;
double b;
if (count%n == 1)
b = 0.0;
else
b = s2 / s1;
for (int i1 = 0; i1 < n; i1++)
dx[i1] = g[i1] + b * dx[i1];
// １次元最適化を行わない
if (opt_1 == 0) {
// 新しい点
for (int i1 = 0; i1 < n; i1++)
x[i1] += step * dx[i1];
// 新しい関数値
y2 = fn(0.0, x, dx);
// 関数値の変化が大きい
if (Math.Abs(y2-y1) > eps) {
y1 = y2;
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y = y2;
}
}
// １次元最適化を行う
else {
// 区間を決める
int sw1 = 0;
double f1 = y1;
double sp = step;
double f2 = fn(sp, x, dx);
if (f2 > f1)
sp = -step;
for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
f2 = fn(sp, x, dx);
if (f2 > f1)
sw1 = 1;
else {
sp *= 2.0;
f1 = f2;
}
}
// 区間が求まらない
if (sw1 == 0)
sw = -2;
// 区間が求まった
else {
// 黄金分割法
double k = gold(0.0, sp, eps, ref y2, ref sw1, max, x, dx, fn);
// 黄金分割法が失敗
if (sw1 < 0)
sw = -3;
// 黄金分割法が成功
else {
// 新しい点
for (int i1 = 0; i1 < n; i1++)
x[i1] += k * dx[i1];
// 関数値の変化が大きい
if (Math.Abs(y1-y2) > eps) {
y1 = y2;
s1 = s2;
}
// 収束（関数値の変化＜ｅｐｓ）
else {
sw = count;
y = y2;
}
}
}
}
}
// 収束（傾き＜ｅｐｓ）
else {
sw = count;
y = y1;
}
}
if (sw == 0)
sw = -1;
return sw;
}
/******************************************/
/* 黄金分割法(与えられた方向での最小値) */
/* a,b : 初期区間 a < b */
/* eps : 許容誤差 */
/* val : 関数値 */
/* ind : 計算状況 */
/* >= 0 : 正常終了(収束回数) */
/* = -1 : 収束せず */
/* max : 最大試行回数 */
/* w : 位置 */
/* dw : 傾きの成分 */
/* fn : 関数値を計算する関数 */
/* return : 結果(w＋y＊dwのy) */
/******************************************/
double gold(double a, double b, double eps, ref double val, ref int ind,
int max, double[] w, double[] dw,
Func<double, double[], double[], double> fn)
{
// 初期設定
double tau = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0;
double x1 = b - tau * (b - a);
double x2 = a + tau * (b - a);
double f1 = fn(x1, w, dw);
double f2 = fn(x2, w, dw);
int count = 0;
ind = -1;
// 計算
while (count < max && ind < 0) {
count += 1;
if (f2 > f1) {
if (Math.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 = fn(x1, w, dw);
}
}
else {
if (Math.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 = fn(x2, w, dw);
}
}
}
// 収束した場合の処理
if (ind == 0) {
ind = count;
double fa = fn(a, w, dw);
double fb = fn(b, w, dw);
if (fb < fa) {
a = b;
fa = fb;
}
if (fa < f1) {
x = a;
val = fa;
}
}
return x;
}
}
''''''''''''''''''''''''''''''''
' 共役勾配法による最小値の計算 '
' coded by Y.Suganuma '
''''''''''''''''''''''''''''''''
Imports System.Text.RegularExpressions
Module Test
Sub Main()
' データの入力
' 1 行目
Dim MS As Regex = New Regex("\s+")
Dim str() As String = MS.Split(Console.ReadLine().Trim())
Dim fun As Integer = Integer.Parse(str(1))
Dim n As Integer = Integer.Parse(str(3))
Dim max As Integer = Integer.Parse(str(5))
Dim opt_1 As Integer = Integer.Parse(str(7))
' 2 行目
str = MS.Split(Console.ReadLine().Trim())
Dim eps As Double = Double.Parse(str(1))
Dim step1 As Double = Double.Parse(str(3))
' 3 行目
Dim x(n) As Double
Dim dx(n) As Double
Dim y As Double = 0.0
str = MS.Split(Console.ReadLine().Trim())
For i1 As Integer = 0 To n-1
x(i1) = Double.Parse(str(i1+1))
Next
' 関数１の値の計算（ラムダ式）
Dim snx1 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = w(0) - 1.0
Dim y1 As Double = w(1) - 2.0
Dim y0 As Double = x1 * x1 + y1 * y1
Return y0
End Function
' 関数２の値の計算（ラムダ式）
Dim snx2 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = w(1) - w(0) * w(0)
Dim y1 As Double = 1.0 - w(0)
Dim y0 As Double = 100.0 * x1 * x1 + y1 * y1
Return y0
End Function
' 関数３の値の計算（ラムダ式）
Dim snx3 = Function(k0, x0, dx0) As Double
Dim w(2) As Double
w(0) = x0(0) + k0 * dx0(0)
w(1) = x0(1) + k0 * dx0(1)
Dim x1 As Double = 1.5 - w(0) * (1.0 - w(1))
Dim y1 As Double = 2.25 - w(0) * (1.0 - w(1) * w(1))
Dim z1 As Double = 2.625 - w(0) * (1.0 - w(1) * w(1) * w(1))
Dim y0 As Double = x1 * x1 + y1 * y1 + z1 * z1
Return y0
End Function
' 関数１の微分（ラムダ式）
Dim dsnx1 = Function(x0, dx0) As Integer
dx0(0) = -2.0 * (x0(0) - 1.0)
dx0(1) = -2.0 * (x0(1) - 2.0)
Return 0
End Function
' 関数２の微分（ラムダ式）
Dim dsnx2 = Function(x0, dx0) As Integer
dx0(0) = 400.0 * x0(0) * (x0(1) - x0(0) * x0(0)) + 2.0 * (1.0 - x0(0))
dx0(1) = -200.0 * (x0(1) - x0(0) * x0(0))
Return 0
End Function
' 関数３の微分（ラムダ式）
Dim dsnx3 = Function(x0, dx0) As Integer
dx0(0) = 2.0 * (1.0 - x0(1)) * (1.5 - x0(0) * (1.0 - x0(1))) +
2.0 * (1.0 - x0(1) * x0(1)) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) +
2.0 * (1.0 - x0(1) * x0(1) * x0(1)) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
dx0(1) = -2.0 * x0(0) * (1.5 - x0(0) * (1.0 - x0(1))) -
4.0 * x0(0) * x0(1) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) -
6.0 * x0(0) * x0(1) * x0(1) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
Return 0
End Function
' 実行
Dim sw As Integer
If fun = 1
sw = conjugate(opt_1, max, n, eps, step1, y, x, dx, snx1, dsnx1)
ElseIf fun = 2
sw = conjugate(opt_1, max, n, eps, step1, y, x, dx, snx2, dsnx2)
Else
sw = conjugate(opt_1, max, n, eps, step1, y, x, dx, snx3, dsnx3)
End If
' 結果の出力
If sw < 0
Console.Write(" 収束しませんでした！")
If sw = -1
Console.WriteLine("（収束回数）")
ElseIf sw = -2
Console.WriteLine("（１次元最適化の区間）")
ElseIf sw = -3
Console.WriteLine("（黄金分割法）")
End If
Else
Console.Write(" 結果＝")
For i1 As Integer = 0 To n-1
Console.Write(x(i1) & " ")
Next
Console.WriteLine(" 最小値＝" & y & " 回数＝" & sw)
End If
End Sub
''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 共役勾配法 '
' opt_1 : =0 : １次元最適化を行わない '
' =1 : １次元最適化を行う '
' max : 最大繰り返し回数 '
' n : 次元 '
' eps : 収束判定条件 '
' step1 : きざみ幅 '
' y : 最小値 '
' x : x(初期値と答え) '
' dx : 関数の微分値 '
' fn : 関数値を計算する関数 '
' dfn : 関数の微分値を計算する関数 '
' return : >=0 : 正常終了(収束回数) '
' =-1 : 収束せず '
' =-2 : １次元最適化の区間が求まらない '
' =-3 : 黄金分割法が失敗 '
''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function conjugate(opt_1 As Integer, max As Integer, n As Integer, eps As Double,
step1 As Double, ByRef y As Double, x() As Double, dx() As Double,
fn As Func(Of Double, Double(), Double(), Double),
dfn As Func(Of Double(), Double(), Integer))
Dim g(n) As Double
Dim s1 As Double = 1.0
Dim s2 As Double
Dim y1 As Double = fn(0.0, x, dx)
Dim y2 As Double = 0.0
Dim count As Integer = 0
Dim sw As Integer = 0
Do While count < max and sw = 0
' 傾きの計算
dfn(x, g)
' 傾きのチェック
s2 = 0.0
For i1 As Integer = 0 To n-1
s2 += g(i1) * g(i1)
Next
' 収束していない
If Math.Sqrt(s2) > eps
' 方向の計算
count += 1
Dim b As Double
If (count Mod n) = 1
b = 0.0
Else
b = s2 / s1
End If
For i1 As Integer = 0 To n-1
dx(i1) = g(i1) + b * dx(i1)
Next
' １次元最適化を行わない
If opt_1 = 0
' 新しい点
For i1 As Integer = 0 To n-1
x(i1) += step1 * dx(i1)
Next
' 新しい関数値
y2 = fn(0.0, x, dx)
' 関数値の変化が大きい
If Math.Abs(y2-y1) > eps
y1 = y2
s1 = s2
' 収束（関数値の変化＜ｅｐｓ）
Else
sw = count
y = y2
End If
' １次元最適化を行う
Else
' 区間を決める
Dim sw1 As Integer = 0
Dim f1 As Double = y1
Dim sp As Double = step1
Dim f2 As Double = fn(sp, x, dx)
If f2 > f1
sp = -step1
End If
Dim ii As Integer = 0
Do While ii < max and sw1 = 0
f2 = fn(sp, x, dx)
If f2 > f1
sw1 = 1
Else
sp *= 2.0
f1 = f2
End If
ii += 1
Loop
' 区間が求まらない
If sw1 = 0
sw = -2
' 区間が求まった
Else
' 黄金分割法
Dim k As Double = gold(0.0, sp, eps, y2, sw1, max, x, dx, fn)
' 黄金分割法が失敗
If sw1 < 0
sw = -3
' 黄金分割法が成功
Else
' 新しい点
For i1 As Integer = 0 To n-1
x(i1) += k * dx(i1)
Next
' 関数値の変化が大きい
If Math.Abs(y1-y2) > eps
y1 = y2
s1 = s2
' 収束（関数値の変化＜ｅｐｓ）
Else
sw = count
y = y2
End If
End If
End If
End If
' 収束（傾き＜ｅｐｓ）
Else
sw = count
y = y1
End If
Loop
If sw = 0
sw = -1
End If
Return sw
End Function
''''''''''''''''''''''''''''''''''''''''''
' 黄金分割法(与えられた方向での最小値) '
' a,b : 初期区間 a < b '
' eps : 許容誤差 '
' val : 関数値 '
' ind : 計算状況 '
' >= 0 : 正常終了(収束回数) '
' = -1 : 収束せず '
' max : 最大試行回数 '
' w : 位置 '
' dw : 傾きの成分 '
' fn : 関数値を計算する関数 '
' return : 結果(w＋y＊dwのy) '
''''''''''''''''''''''''''''''''''''''''''
Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
ByRef ind As Integer, max As Integer, w() As Double, dw() As Double,
fn As Func(Of Double, Double(), Double(), Double))
' 初期設定
Dim tau As Double = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0
Dim x1 As Double = b - tau * (b - a)
Dim x2 As Double = a + tau * (b - a)
Dim f1 As Double = fn(x1, w, dw)
Dim f2 As Double = fn(x2, w, dw)
Dim count As Integer = 0
ind = -1
' 計算
Do While count < max and ind < 0
count += 1
If f2 > f1
If Math.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 = fn(x1, w, dw)
End If
Else
If Math.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 = fn(x2, w, dw)
End If
End If
Loop
' 収束した場合の処理
If ind = 0
ind = count
Dim fa As Double = fn(a, w, dw)
Dim fb As Double = fn(b, w, dw)
If fb < fa
a = b
fa = fb
End If
If fa < f1
x = a
val = fa
End If
End If
Return x
End Function
End Module
| 情報学部 | 菅沼ホーム | 目次 | 索引 |