最適化（多項式近似法）

/*********************************************/
/* 多項式近似による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 eps, x, d, val[] = new double [1];
		int max, ind[] = new int [1];

		x   = -2.0;
		d   = 0.1;
		eps = 1.0e-10;
		max = 100;

		Kansu kn = new Kansu();
		x = kn.approx(x, d, eps, val, ind, max);

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

/****************/
/* 関数値の計算 */
/****************/
class Kansu extends Approximation
{
	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 Approximation {

	/************************************************/
	/* 多項式近似(関数の最小値)                     */
	/*      x0  : 初期値                            */
	/*      d0  : 初期ステップ                      */
	/*      eps : 許容誤差                          */
	/*      val : 間数値                            */
	/*      ind : 計算状況                          */
	/*              >= 0 : 正常終了(収束回数)       */
	/*              = -1 : 収束せず                 */
	/*      max : 最大試行回数                      */
	/*      return : 結果                           */
	/************************************************/

	abstract double snx(double x);

	double approx(double x0, double d0, double eps, double val[], int ind[], int max)
	{
		double f[] = new double [4];
		double x[] = new double [4];
		double xx = 0.0, d, dl;
		int i1, k = 0, count = 0, min, sw = 0;

		d      = d0;
		x[1]   = x0;
		f[1]   = snx(x0);
		ind[0] = -1;

		while (count < max && ind[0] < 0) {
			x[3] = x[1] + d;
			f[3] = snx(x[3]);
			while (k < max && f[3] <= f[1]) {
				k++;
				d *= 2.0;
				x[0] = x[1];
				f[0] = f[1];
				x[1] = x[3];
				f[1] = f[3];
				x[3] = x[1] + d;
				f[3] = snx(x[3]);
			}
						// 初期値が不適当
			if (k >= max)
				count = max;
			else {
						// ３点の選択
				sw = 0;
				if (k > 0) {
					x[2] = x[3] - 0.5 * d;
					f[2] = snx(x[2]);
					min  = -1;
					for (i1 = 0; i1 < 4; i1++) {
						if (min < 0 || f[i1] < f[min])
							min = i1;
					}
					if (min >= 2) {
						for (i1 = 0; i1 < 3; i1++) {
							x[i1] = x[i1+1];
							f[i1] = f[i1+1];
						}
					}
					sw = 1;
				}
				else {
					x[0] = x[1] - d0;
					f[0] = snx(x[0]);
					if (f[0] > f[1]) {
						x[2] = x[3];
						f[2] = f[3];
						sw = 1;
					}
					else {
						x[1] = x[0];
						f[1] = f[0];
						d0   = -d0;
						d    = 2.0 * d0;
						k    = 1;
					}
				}
						// 収束？
				if (sw > 0) {
					count++;
					dl = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2]);
					xx     = x[1] - dl;
					val[0] = snx(xx);
					if (Math.abs(dl) < eps)
						ind[0] = count;
					else {
						k  = 0;
						d0 = 0.5 * d;
						d  = d0;
						if (val[0] < f[1]) {
							x[1] = xx;
							f[1] = val[0];
						}
					}
				}
			}
		}

		return xx;
	}
}