最適化(黄金分割法)

/*********************************************/
/* 黄金分割法によるy=x^4+3x^3+2x^2+1の最小値 */
/*        coded by Y.Suganuma                */
/*********************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double a, b, eps, x, val[] = new double [1];
		int max, ind[] = new int [1];

		a   = -2.0;
		b   = -1.0;
		eps = 1.0e-10;
		max = 100;

		Kansu kn = new Kansu();
		x = kn.gold(a, b, eps, val, ind, max);

		System.out.println("x " + x + " val " + val[0] + " ind " + ind[0]);
	}
}

/****************/
/* 関数値の計算 */
/****************/
class Kansu extends Gold
{
	double snx(double x)
	{
		double y = x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0;
		return y;
	}
}

abstract class Gold 
{
	/************************************************/
	/* 黄金分割法(関数の最小値)                     */
	/*      a,b : 初期区間 a < b                    */
	/*      eps : 許容誤差                          */
	/*      val : 間数値                            */
	/*      ind : 計算状況                          */
	/*              >= 0 : 正常終了(収束回数)       */
	/*              = -1 : 収束せず                 */
	/*      max : 最大試行回数                      */
	/*      return : 結果                           */
	/************************************************/

	abstract double snx(double x);

	double gold(double a, double b, double eps, double val[], int ind[], int max)
	{
		double f1, f2, fa, fb, tau, x = 0.0, x1, x2;
		int count;

		tau    = (Math.sqrt(5.0) - 1.0) / 2.0;
		count  = 0;
		ind[0] = -1;
		x1     = b - tau * (b - a);
		x2     = a + tau * (b - a);
		f1     = snx(x1);
		f2     = snx(x2);

		while (count < max && ind[0] < 0) {
			count += 1;
			if (f2 > f1) {
				if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
					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);
				}
			}
			else {
				if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
					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);
				}
			}
		}

		if (ind[0] == 0) {
			ind[0] = count;
			fa     = snx(a);
			fb     = snx(b);
			if (fb < fa) {
				a  = b;
				fa = fb;
			}
			if (fa < f1) {
				x      = a;
				val[0] = fa;
			}
		}

		return x;
	}
}