データ例とプログラム
-------------------------i_data-------------------------
// 関数 a，一次元最適化を使用しない
関数 1 変数の数 2 最大試行回数 100 一次元最適化 0
許容誤差 1.0e-10 刻み幅 1.0
初期値 0.0 0.0
// 関数 a，一次元最適化を使用する
関数 1 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 1.0
初期値 0.0 0.0
// 関数 b，一次元最適化を使用しない
関数 2 変数の数 2 最大試行回数 100 一次元最適化 0
許容誤差 1.0e-10 刻み幅 1.0
初期値 0.0 0.0
// 関数 b，一次元最適化を使用する
関数 2 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 0.1
初期値 0.0 0.0
// 関数 c，一次元最適化を使用しない
関数 3 変数の数 2 最大試行回数 100 一次元最適化 0
許容誤差 1.0e-10 刻み幅 1.0
初期値 1.0 0.0
// 関数 c，一次元最適化を使用する
関数 3 変数の数 2 最大試行回数 100 一次元最適化 1
許容誤差 1.0e-10 刻み幅 1.0
初期値 1.0 0.0
-------------------------プログラム-------------------------
/******************************/
/* Newton法による最小値の計算 */
/* coded by Y.Suganuma */
/******************************/
import java.io.*;
import java.util.*;
public class Test {
public static void main(String args[]) throws IOException
{
double eps, H[][], 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];
H = new double [n][2*n];
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.Newton(opt_1, max, n, eps, step, y, x, dx, H);
// 結果の出力
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 Newton
{
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;
}
}
// Hesse行列の逆行列
int hesse(double x[], double H[][], double eps)
{
double x1, x2, x3;
int sw1 = 0;
switch (sw) {
case 1:
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;
sw1 = gauss(H, 2, 2, eps);
break;
case 2:
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;
sw1 = gauss(H, 2, 2, eps);
break;
case 3:
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;
sw1 = gauss(H, 2, 2, eps);
break;
}
return sw1;
}
}
abstract class Newton {
/********************************************************/
/* Newton法 */
/* opt_1 : =0 : １次元最適化を行わない */
/* =1 : １次元最適化を行う */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* eps : 収束判定条件 */
/* step : きざみ幅 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* dx : 関数の微分値 */
/* H : Hesse行列の逆行列 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/* =-2 : １次元最適化の区間が求まらない */
/* =-3 : 黄金分割法が失敗 */
/********************************************************/
abstract double snx(double k, double x[], double dx[]);
abstract void dsnx(double x[], double dx[]);
abstract int hesse(double x[], double H[][], double eps);
int Newton(int opt_1, int max, int n, double eps, double step, double y[], double x[], double dx[], double H[][])
{
double f1, f2, k, sp, y1[] = new double [1], y2[] = new double [1];
double wk[] = new double [n];
int count = 0, i1, i2, sw = 0, sw1[] = new int [1];
y1[0] = snx(0.0, x, dx);
while (count < max && sw == 0) {
// 傾きの計算
dsnx(x, wk);
// Hesse行列の逆行列の計算
sw1[0] = hesse(x, H, eps);
// 収束していない
if (sw1[0] == 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[0] = snx(0.0, x, dx);
// 関数値の変化が大きい
if (Math.abs(y2[0]-y1[0]) > eps)
y1[0] = y2[0];
// 収束（関数値の変化＜ｅｐｓ）
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];
// 収束（関数値の変化＜ｅｐｓ）
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;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
int gauss(double w[][], int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
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 = Math.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);
}
}