静岡理工科大学 菅沼ホーム 目次 索引

最適化(準 Newton 法)

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,準 Newton 法を使用して,非線形関数の最小値を求めるためのものです( プログラムの使用方法).

  1. C++

    /**********************************/
    /* 準Newton法によるの最小値の計算 */
    /*      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 DFP_BFGS(int, int, int, int, double, double, double *, double *, double *, double **,
              double (*)(double, double *, double *),
              void (*)(double *, double *));
    
    int main()
    {
    	double eps, **H, step, *x, *dx, y;
    	int fun, i1, max, method, n, opt_1, sw = 0;
    					// データの入力
    	scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
    	scanf("%*s %d %*s %lf %*s %lf", &method, &eps, &step);
    	x  = new double [n];
    	dx = new double [n];
    	H  = new double * [n];
    	scanf("%*s");
    	for (i1 = 0; i1 < n; i1++) {
    		scanf("%lf", &x[i1]);
    		H[i1] = new double [n];
    	}
    					// 実行
    	switch (fun) {
    		case 1:
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx1, dsnx1);
    			break;
    		case 2:
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx2, dsnx2);
    			break;
    		case 3:
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, &y, x, dx, H, snx3, dsnx3);
    			break;
    	}
    					// 結果の出力
    	if (sw < 0) {
    		printf("   収束しませんでした!");
    		switch (sw) {
    			case -1:
    				printf("(収束回数)\n");
    				break;
    			case -2:
    				printf("(1次元最適化の区間)\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;
    }
    
    /********************************************************/
    /* DFP法 or BFGS法                                      */
    /*      method : =0 : DFP法                             */
    /*               =1 : BFGS法                            */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x1 : x(初期値と答え)                            */
    /*      dx : 関数の微分値                               */
    /*      H : Hesse行列の逆行列                           */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    #include <math.h>
    
    double gold(double, double, double, double *, int *, int, double *, double *,
                double (*)(double, double *, double *));
    
    int DFP_BFGS(int method, int opt_1, int max, int n, double eps, double step, double *y,
              double *x1, double *dx, double **H, double (*snx)(double, double *, double *), 
              void (*dsnx)(double *, double *))
    {
    	double f1, f2, *g, *g1, *g2, **H1, **H2, **I, k, *s, sm, sm1, sm2,
               sp, *w, *x2, y1, y2;
    	int count = 0, i1, i2, i3, sw = 0, sw1;
    
    	x2 = new double [n];
    	g  = new double [n];
    	g1 = new double [n];
    	g2 = new double [n];
    	s  = new double [n];
    	w  = new double [n];
    
    	y1 = snx(0.0, x1, dx);
    	dsnx(x1, g1);
    
    	H1 = new double * [n];
    	H2 = new double * [n];
    	I  = new double * [n];
    	for (i1 = 0; i1 < n; i1++) {
    		H1[i1] = new double [n];
    		H2[i1] = new double [n];
    		I[i1]  = new double [n];
    		for (i2 = 0; i2 < n; i2++) {
    			H[i1][i2] = 0.0;
    			I[i1][i2] = 0.0;
    		}
    		H[i1][i1] = 1.0;
    		I[i1][i1] = 1.0;
    	}
    
    	while (count < max && sw == 0) {
    		count++;
    					// 方向の計算
    		for (i1 = 0; i1 < n; i1++) {
    			dx[i1] = 0.0;
    			for (i2 = 0; i2 < n; i2++)
    				dx[i1] -= H[i1][i2] * g1[i2];
    		}
    					// 1次元最適化を行わない
    		if (opt_1 == 0) {
    						// 新しい点
    			for (i1 = 0; i1 < n; i1++)
    				x2[i1] = x1[i1] + step * dx[i1];
    						// 新しい関数値
    			y2 = snx(0.0, x2, dx);
    		}
    					// 1次元最適化を行う
    		else {
    						// 区間を決める
    			sw1 = 0;
    			f1  = y1;
    			sp  = step;
    			f2  = snx(sp, x1, dx);
    			if (f2 > f1)
    				sp = -step;
    			for (i1 = 0; i1 < max && sw1 == 0; i1++) {
    				f2 = snx(sp, x1, 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, x1, dx, snx);
    							// 黄金分割法が失敗
    				if (sw1 < 0)
    					sw = -3;
    							// 黄金分割法が成功
    				else {
    								// 新しい点
    					for (i1 = 0; i1 < n; i1++)
    						x2[i1] = x1[i1] + k * dx[i1];
    				}
    			}
    		}
    
    		if (sw == 0) {
    					// 収束(関数値の変化<eps)
    			if (fabs(y2-y1) < eps) {
    				sw = count;
    				*y = y2;
    				for (i1 = 0; i1 < n; i1++)
    					x1[i1] = x2[i1];
    			}
    					// 関数値の変化が大きい
    			else {
    						// 傾きの計算
    				dsnx(x2, g2);
    				sm = 0.0;
    				for (i1 = 0; i1 < n; i1++)
    					sm += g2[i1] * g2[i1];
    				sm = sqrt(sm);
    						// 収束(傾き<eps)
    				if (sm < eps) {
    					sw = count;
    					*y = y2;
    					for (i1 = 0; i1 < n; i1++)
    						x1[i1] = x2[i1];
    				}
    						// 収束していない
    				else {
    					for (i1 = 0; i1 < n; i1++) {
    						g[i1] = g2[i1] - g1[i1];
    						s[i1] = x2[i1] - x1[i1];
    					}
    					sm1 = 0.0;
    					for (i1 = 0; i1 < n; i1++)
    						sm1 += s[i1] * g[i1];
    							// DFP法
    					if (method == 0) {
    						for (i1 = 0; i1 < n; i1++) {
    							w[i1] = 0.0;
    							for (i2 = 0; i2 < n; i2++)
    								w[i1] += g[i2] * H[i2][i1];
    						}
    						sm2 = 0.0;
    						for (i1 = 0; i1 < n; i1++)
    							sm2 += w[i1] * g[i1];
    						for (i1 = 0; i1 < n; i1++) {
    							w[i1] = 0.0;
    							for (i2 = 0; i2 < n; i2++)
    								w[i1] += H[i1][i2] * g[i2];
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++)
    								H1[i1][i2] = w[i1] * g[i2];
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++) {
    								H2[i1][i2] = 0.0;
    								for (i3 = 0; i3 < n; i3++)
    									H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    							}
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++)
    								H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
    						}
    					}
    							// BFGS法
    					else {
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++)
    								H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++) {
    								H2[i1][i2] = 0.0;
    								for (i3 = 0; i3 < n; i3++)
    									H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    							}
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++) {
    								H[i1][i2] = 0.0;
    								for (i3 = 0; i3 < n; i3++)
    									H[i1][i2] += H2[i1][i3] * H1[i3][i2];
    							}
    						}
    						for (i1 = 0; i1 < n; i1++) {
    							for (i2 = 0; i2 < n; i2++)
    								H[i1][i2] += s[i1] * s[i2] / sm1;
    						}
    					}
    					y1 = y2;
    					for (i1 = 0; i1 < n; i1++) {
    						g1[i1] = g2[i1];
    						x1[i1] = x2[i1];
    					}
    				}
    			}
    		}
    	}
    
    	if (sw == 0)
    		sw = -1;
    
    	delete [] x2;
    	delete [] g;
    	delete [] g1;
    	delete [] g2;
    	delete [] s;
    	delete [] w;
    
    	for (i1 = 0; i1 < n; i1++) {
    		delete [] H1[i1];
    		delete [] H2[i1];
    		delete [] I[i1];
    	}
    	delete [] H1;
    	delete [] H2;
    	delete [] I;
    
    	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だけ進んだ */
    /* 点における関数値を計算する                */
    /*********************************************/
    					// 関数1
    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;
    }
    					// 関数2
    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;
    }
    					// 関数3
    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;
    }
    
    /********************/
    /* 微係数を計算する */
    /********************/
    					// 関数1
    void dsnx1(double *x, double *dx)
    {
    	dx[0] = 2.0 * (x[0] - 1.0);
    	dx[1] = 2.0 * (x[1] - 2.0);
    }
    					// 関数2
    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]);
    }
    					// 関数3
    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]));
    }
    			

  2. Java

    /**********************************/
    /* 準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, method, 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();
    		method = Integer.parseInt(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][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.DFP_BFGS(method, 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("(1次元最適化の区間)");
    					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 P_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;
    		}
    	}
    }
    
    abstract class P_Newton {
    
    	/********************************************************/
    	/* DFP法 or BFGS法                                      */
    	/*      method : =0 : DFP法                             */
    	/*               =1 : BFGS法                            */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x1 : x(初期値と答え)                            */
    	/*      dx : 関数の微分値                               */
    	/*      H : Hesse行列の逆行列                           */
    	/*      return : >=0 : 正常終了(収束回数)               */
    	/*               =-1 : 収束せず                         */
    	/*               =-2 : 1次元最適化の区間が求まらない   */
    	/*               =-3 : 黄金分割法が失敗                 */
    	/********************************************************/
    
    	abstract double snx(double k, double x[], double dx[]);
    	abstract void dsnx(double x[], double dx[]);
    
    	int DFP_BFGS(int method, int opt_1, int max, int n, double eps, double step, double y[], double x1[], double dx[], double H[][])
    	{
    		double f1, f2, g[], g1[], g2[], H1[][], H2[][], I[][], k, s[], sm, sm1, sm2,
                   sp, w[], x2[], y1[] = new double [1], y2[] = new double [1];
    		int count = 0, i1, i2, i3, sw = 0, sw1[] = new int [1];
    
    		x2 = new double [n];
    		g  = new double [n];
    		g1 = new double [n];
    		g2 = new double [n];
    		s  = new double [n];
    		w  = new double [n];
    
    		y1[0] = snx(0.0, x1, dx);
    		dsnx(x1, g1);
    
    		H1 = new double [n][n];
    		H2 = new double [n][n];
    		I  = new double [n][n];
    		for (i1 = 0; i1 < n; i1++) {
    			for (i2 = 0; i2 < n; i2++) {
    				H[i1][i2] = 0.0;
    				I[i1][i2] = 0.0;
    			}
    			H[i1][i1] = 1.0;
    			I[i1][i1] = 1.0;
    		}
    
    		while (count < max && sw == 0) {
    			count++;
    					// 方向の計算
    			for (i1 = 0; i1 < n; i1++) {
    				dx[i1] = 0.0;
    				for (i2 = 0; i2 < n; i2++)
    					dx[i1] -= H[i1][i2] * g1[i2];
    			}
    					// 1次元最適化を行わない
    			if (opt_1 == 0) {
    						// 新しい点
    				for (i1 = 0; i1 < n; i1++)
    					x2[i1] = x1[i1] + step * dx[i1];
    						// 新しい関数値
    				y2[0] = snx(0.0, x2, dx);
    			}
    					// 1次元最適化を行う
    			else {
    						// 区間を決める
    				sw1[0] = 0;
    				f1     = y1[0];
    				sp     = step;
    				f2     = snx(sp, x1, dx);
    				if (f2 > f1)
    					sp = -step;
    				for (i1 = 0; i1 < max && sw1[0] == 0; i1++) {
    					f2 = snx(sp, x1, 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, x1, dx);
    							// 黄金分割法が失敗
    					if (sw1[0] < 0)
    						sw = -3;
    							// 黄金分割法が成功
    					else {
    								// 新しい点
    						for (i1 = 0; i1 < n; i1++)
    							x2[i1] = x1[i1] + k * dx[i1];
    					}
    				}
    			}
    
    			if (sw == 0) {
    					// 収束(関数値の変化<eps)
    				if (Math.abs(y2[0]-y1[0]) < eps) {
    					sw   = count;
    					y[0] = y2[0];
    					for (i1 = 0; i1 < n; i1++)
    						x1[i1] = x2[i1];
    				}
    					// 関数値の変化が大きい
    				else {
    						// 傾きの計算
    					dsnx(x2, g2);
    					sm = 0.0;
    					for (i1 = 0; i1 < n; i1++)
    						sm += g2[i1] * g2[i1];
    					sm = Math.sqrt(sm);
    						// 収束(傾き<eps)
    					if (sm < eps) {
    						sw   = count;
    						y[0] = y2[0];
    						for (i1 = 0; i1 < n; i1++)
    							x1[i1] = x2[i1];
    					}
    						// 収束していない
    					else {
    						for (i1 = 0; i1 < n; i1++) {
    							g[i1] = g2[i1] - g1[i1];
    							s[i1] = x2[i1] - x1[i1];
    						}
    						sm1 = 0.0;
    						for (i1 = 0; i1 < n; i1++)
    							sm1 += s[i1] * g[i1];
    							// DFP法
    						if (method == 0) {
    							for (i1 = 0; i1 < n; i1++) {
    								w[i1] = 0.0;
    								for (i2 = 0; i2 < n; i2++)
    									w[i1] += g[i2] * H[i2][i1];
    							}
    							sm2 = 0.0;
    							for (i1 = 0; i1 < n; i1++)
    								sm2 += w[i1] * g[i1];
    							for (i1 = 0; i1 < n; i1++) {
    								w[i1] = 0.0;
    								for (i2 = 0; i2 < n; i2++)
    									w[i1] += H[i1][i2] * g[i2];
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++)
    									H1[i1][i2] = w[i1] * g[i2];
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++) {
    									H2[i1][i2] = 0.0;
    									for (i3 = 0; i3 < n; i3++)
    										H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    								}
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++)
    									H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
    							}
    						}
    							// BFGS法
    						else {
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++)
    									H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++) {
    									H2[i1][i2] = 0.0;
    									for (i3 = 0; i3 < n; i3++)
    										H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    								}
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++) {
    									H[i1][i2] = 0.0;
    									for (i3 = 0; i3 < n; i3++)
    										H[i1][i2] += H2[i1][i3] * H1[i3][i2];
    								}
    							}
    							for (i1 = 0; i1 < n; i1++) {
    								for (i2 = 0; i2 < n; i2++)
    									H[i1][i2] += s[i1] * s[i2] / sm1;
    							}
    						}
    						y1[0] = y2[0];
    						for (i1 = 0; i1 < n; i1++) {
    							g1[i1] = g2[i1];
    							x1[i1] = x2[i1];
    						}
    					}
    				}
    			}
    		}
    
    		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;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,JavaScript の仕様に適合した形で最小値を求めたい式を入力することによって,任意の関数の最小値を画面上で求めることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>最適化(準 Newton 法)</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 = 0;
    			if (document.getElementById("line").options[1].selected)
    				opt_1 = 1;
    			let method = 0;
    			if (document.getElementById("method").options[1].selected)
    				method = 1;
    			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();
    			let H   = new Array();
    			for (let i1 = 0; i1 < n; i1++) {
    				dx[i1] = 0.0;
    				H[i1]  = new Array();
    			}
    					// 実行
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx);
    					// 結果
    			let res;
    			if (sw < 0) {
    				res = " 収束しませんでした!";
    				switch (sw) {
    					case -1:
    						res += "(収束回数)";
    						break;
    					case -2:
    						res += "(1次元最適化の区間)";
    						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]);
    		}
    
    		/******************************************************/
    		/* DFP法 or BFGS法                                    */
    		/*      method : =0 : DFP法                           */
    		/*               =1 : BFGS法                          */
    		/*      opt_1 : =0 : 1次元最適化を行わない           */
    		/*              =1 : 1次元最適化を行う               */
    		/*      max : 最大繰り返し回数                        */
    		/*      n : 次元                                      */
    		/*      eps : 収束判定条件                            */
    		/*      step : きざみ幅                               */
    		/*      y : 最小値                                    */
    		/*      x1 : x(初期値と答え)                          */
    		/*      dx : 関数の微分値                             */
    		/*      H : Hesse行列の逆行列                         */
    		/*      fn : 関数値を計算する関数                     */
    		/*      dfn : 関数の微分値を計算する関数              */
    		/*      return : >=0 : 正常終了(収束回数)             */
    		/*               =-1 : 収束せず                       */
    		/*               =-2 : 1次元最適化の区間が求まらない */
    		/*               =-3 : 黄金分割法が失敗               */
    		/******************************************************/
    		function DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, H, fn, dfn)
    		{
    			let f1, f2, k, sm, sm1, sm2, sp, count = 0, sw = 0;
    			let g = new Array();
    			let s = new Array();
    			let w = new Array();
    			let x2 = new Array();
    			let y1 = new Array();
    			let y2 = new Array();
    			let g1 = new Array();
    			let g2 = new Array();
    			let sw1 = new Array();
    			let H1 = new Array();
    			let H2 = new Array();
    			let I  = new Array();
    			for (let i1 = 0; i1 < n; i1++) {
    				H1[i1] = new Array();
    				H2[i1] = new Array();
    				I[i1]  = new Array();
    			}
    
    			y1[0] = fn(0.0, x1, dx);
    			dfn(n, x1, g1);
    
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++) {
    					H[i1][i2] = 0.0;
    					I[i1][i2] = 0.0;
    				}
    				H[i1][i1] = 1.0;
    				I[i1][i1] = 1.0;
    			}
    
    			while (count < max && sw == 0) {
    				count++;
    					// 方向の計算
    				for (let i1 = 0; i1 < n; i1++) {
    					dx[i1] = 0.0;
    					for (let i2 = 0; i2 < n; i2++)
    						dx[i1] -= H[i1][i2] * g1[i2];
    				}
    					// 1次元最適化を行わない
    				if (opt_1 == 0) {
    						// 新しい点
    					for (let i1 = 0; i1 < n; i1++)
    						x2[i1] = x1[i1] + step * dx[i1];
    						// 新しい関数値
    					y2[0] = fn(0.0, x2, dx);
    				}
    					// 1次元最適化を行う
    				else {
    						// 区間を決める
    					sw1[0] = 0;
    					f1     = y1[0];
    					sp     = step;
    					f2     = fn(sp, x1, dx);
    					if (f2 > f1)
    						sp = -step;
    					for (let i1 = 0; i1 < max && sw1[0] == 0; i1++) {
    						f2 = fn(sp, x1, 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, x1, dx, snx);
    							// 黄金分割法が失敗
    						if (sw1[0] < 0)
    							sw = -3;
    							// 黄金分割法が成功
    						else {
    								// 新しい点
    							for (let i1 = 0; i1 < n; i1++)
    								x2[i1] = x1[i1] + k * dx[i1];
    						}
    					}
    				}
    
    				if (sw == 0) {
    					// 収束(関数値の変化<eps)
    					if (Math.abs(y2[0]-y1[0]) < eps) {
    						sw   = count;
    						y[0] = y2[0];
    						for (let i1 = 0; i1 < n; i1++)
    							x1[i1] = x2[i1];
    					}
    					// 関数値の変化が大きい
    					else {
    						// 傾きの計算
    						dfn(n, x2, g2);
    						sm = 0.0;
    						for (let i1 = 0; i1 < n; i1++)
    							sm += g2[i1] * g2[i1];
    						sm = Math.sqrt(sm);
    						// 収束(傾き<eps)
    						if (sm < eps) {
    							sw   = count;
    							y[0] = y2[0];
    							for (let i1 = 0; i1 < n; i1++)
    								x1[i1] = x2[i1];
    						}
    						// 収束していない
    						else {
    							for (let i1 = 0; i1 < n; i1++) {
    								g[i1] = g2[i1] - g1[i1];
    								s[i1] = x2[i1] - x1[i1];
    							}
    							sm1 = 0.0;
    							for (let i1 = 0; i1 < n; i1++)
    								sm1 += s[i1] * g[i1];
    							// DFP法
    							if (method == 0) {
    								for (let i1 = 0; i1 < n; i1++) {
    									w[i1] = 0.0;
    									for (let i2 = 0; i2 < n; i2++)
    										w[i1] += g[i2] * H[i2][i1];
    								}
    								sm2 = 0.0;
    								for (let i1 = 0; i1 < n; i1++)
    									sm2 += w[i1] * g[i1];
    								for (let i1 = 0; i1 < n; i1++) {
    									w[i1] = 0.0;
    									for (let i2 = 0; i2 < n; i2++)
    										w[i1] += H[i1][i2] * g[i2];
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++)
    										H1[i1][i2] = w[i1] * g[i2];
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++) {
    										H2[i1][i2] = 0.0;
    										for (let i3 = 0; i3 < n; i3++)
    											H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    									}
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++)
    										H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
    								}
    							}
    							// BFGS法
    							else {
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++)
    										H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++) {
    										H2[i1][i2] = 0.0;
    										for (let i3 = 0; i3 < n; i3++)
    											H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    									}
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++) {
    										H[i1][i2] = 0.0;
    										for (let i3 = 0; i3 < n; i3++)
    											H[i1][i2] += H2[i1][i3] * H1[i3][i2];
    									}
    								}
    								for (let i1 = 0; i1 < n; i1++) {
    									for (let i2 = 0; i2 < n; i2++)
    										H[i1][i2] += s[i1] * s[i2] / sm1;
    								}
    							}
    							y1[0] = y2[0];
    							for (let i1 = 0; i1 < n; i1++) {
    								g1[i1] = g2[i1];
    								x1[i1] = x2[i1];
    							}
    						}
    					}
    				}
    			}
    
    			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>最適化(準 Newton 法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,以下に示すような関数の最小値を求める場合に対する値( BFGS 法において,一次元最適化を行わない場合は刻み幅を 0.02,行う場合は 0.002 にしてみて下さい)が設定されています(他の問題を実行する場合は,それらを適切に修正してください).
    		<P STYLE="text-align:center"><IMG SRC="s_newton.gif"></P>
    		<DT>n 変数関数 f(<B>x</B>) に対する準 Newton 法においては,現在の点 <B>x</B> における傾き,
    		<P STYLE="text-align:center"><IMG SRC="s_newton1.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="1000"> 
    		刻み幅:<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.2"><BR>
    		初期値(n個):<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="20" VALUE="0.0 0.0"> 
    		1次元最適化:<SELECT STYLE="font-size: 100%" SIZE="1" ID="line">
    							<OPTION VALUE="no"> 行わない </OPTION><BR>
    							<OPTION VALUE="yes"> 行う </OPTION><BR>
    					  </SELECT> 
    		方法:<SELECT STYLE="font-size: 100%" SIZE="1" ID="method">
    							<OPTION VALUE="no"> DFP法 </OPTION><BR>
    							<OPTION VALUE="yes"> BFGS法 </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>
    			

  4. PHP

    <?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 %d %*s %lf %*s %lf", $method, $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($n);
    	}
    					// 実行
    	switch ($fun) {
    		case 1:
    			$sw = DFP_BFGS($method, $opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx1", "dsnx1");
    			break;
    		case 2:
    			$sw = DFP_BFGS($method, $opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx2", "dsnx2");
    			break;
    		case 3:
    			$sw = DFP_BFGS($method, $opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx3", "dsnx3");
    			break;
    	}
    					// 結果の出力
    	if ($sw < 0) {
    		printf("   収束しませんでした!");
    		switch ($sw) {
    			case -1:
    				printf("(収束回数)\n");
    				break;
    			case -2:
    				printf("(1次元最適化の区間)\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);
    	}
    
    /********************************************************/
    /* DFP法 or BFGS法                                      */
    /*      method : =0 : DFP法                             */
    /*               =1 : BFGS法                            */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x1 : x(初期値と答え)                            */
    /*      dx : 関数の微分値                               */
    /*      H : Hesse行列の逆行列                           */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    
    function DFP_BFGS($method, $opt_1, $max, $n, $eps, $step, &$y, &$x1, $dx, $H, $snx, $dsnx)
    {
    	$count = 0;
    	$sw    = 0;
    	$x2    = array($n);
    	$g     = array($n);
    	$g1    = array($n);
    	$g2    = array($n);
    	$s     = array($n);
    	$w     = array($n);
    
    	$y1 = $snx(0.0, $x1, $dx);
    	$dsnx($x1, $g1);
    
    	$H1 = array($n);
    	$H2 = array($n);
    	$I  = array($n);
    	for ($i1 = 0; $i1 < $n; $i1++) {
    		$H1[$i1] = array($n);
    		$H2[$i1] = array($n);
    		$I[$i1]  = array($n);
    		for ($i2 = 0; $i2 < $n; $i2++) {
    			$H[$i1][$i2] = 0.0;
    			$I[$i1][$i2] = 0.0;
    		}
    		$H[$i1][$i1] = 1.0;
    		$I[$i1][$i1] = 1.0;
    	}
    
    	while ($count < $max && $sw == 0) {
    		$count++;
    					// 方向の計算
    		for ($i1 = 0; $i1 < $n; $i1++) {
    			$dx[$i1] = 0.0;
    			for ($i2 = 0; $i2 < $n; $i2++)
    				$dx[$i1] -= $H[$i1][$i2] * $g1[$i2];
    		}
    					// 1次元最適化を行わない
    		if ($opt_1 == 0) {
    						// 新しい点
    			for ($i1 = 0; $i1 < $n; $i1++)
    				$x2[$i1] = $x1[$i1] + $step * $dx[$i1];
    						// 新しい関数値
    			$y2 = $snx(0.0, $x2, $dx);
    		}
    					// 1次元最適化を行う
    		else {
    						// 区間を決める
    			$sw1 = 0;
    			$f1  = $y1;
    			$sp  = $step;
    			$f2  = $snx($sp, $x1, $dx);
    			if ($f2 > $f1)
    				$sp = -$step;
    			for ($i1 = 0; $i1 < $max && $sw1 == 0; $i1++) {
    				$f2 = $snx($sp, $x1, $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, $x1, $dx, $snx);
    							// 黄金分割法が失敗
    				if ($sw1 < 0)
    					$sw = -3;
    							// 黄金分割法が成功
    				else {
    								// 新しい点
    					for ($i1 = 0; $i1 < $n; $i1++)
    						$x2[$i1] = $x1[$i1] + $k * $dx[$i1];
    				}
    			}
    		}
    
    		if ($sw == 0) {
    					// 収束(関数値の変化<eps)
    			if (abs($y2-$y1) < $eps) {
    				$sw = $count;
    				$y  = $y2;
    				for ($i1 = 0; $i1 < $n; $i1++)
    					$x1[$i1] = $x2[$i1];
    			}
    					// 関数値の変化が大きい
    			else {
    						// 傾きの計算
    				$dsnx($x2, $g2);
    				$sm = 0.0;
    				for ($i1 = 0; $i1 < $n; $i1++)
    					$sm += $g2[$i1] * $g2[$i1];
    				$sm = sqrt($sm);
    						// 収束(傾き<eps)
    				if ($sm < $eps) {
    					$sw = $count;
    					$y  = $y2;
    					for ($i1 = 0; $i1 < $n; $i1++)
    						$x1[$i1] = $x2[$i1];
    				}
    						// 収束していない
    				else {
    					for ($i1 = 0; $i1 < $n; $i1++) {
    						$g[$i1] = $g2[$i1] - $g1[$i1];
    						$s[$i1] = $x2[$i1] - $x1[$i1];
    					}
    					$sm1 = 0.0;
    					for ($i1 = 0; $i1 < $n; $i1++)
    						$sm1 += $s[$i1] * $g[$i1];
    							// DFP法
    					if ($method == 0) {
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							$w[$i1] = 0.0;
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$w[$i1] += $g[$i2] * $H[$i2][$i1];
    						}
    						$sm2 = 0.0;
    						for ($i1 = 0; $i1 < $n; $i1++)
    							$sm2 += $w[$i1] * $g[$i1];
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							$w[$i1] = 0.0;
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$w[$i1] += $H[$i1][$i2] * $g[$i2];
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$H1[$i1][$i2] = $w[$i1] * $g[$i2];
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++) {
    								$H2[$i1][$i2] = 0.0;
    								for ($i3 = 0; $i3 < $n; $i3++)
    									$H2[$i1][$i2] += $H1[$i1][$i3] * $H[$i3][$i2];
    							}
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$H[$i1][$i2] = $H[$i1][$i2]  + $s[$i1] * $s[$i2] / $sm1 - $H2[$i1][$i2] / $sm2;
    						}
    					}
    							// BFGS法
    					else {
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$H1[$i1][$i2] = $I[$i1][$i2] - $s[$i1] * $g[$i2] / $sm1;
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++) {
    								$H2[$i1][$i2] = 0.0;
    								for ($i3 = 0; $i3 < $n; $i3++)
    									$H2[$i1][$i2] += $H1[$i1][$i3] * $H[$i3][$i2];
    							}
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++) {
    								$H[$i1][$i2] = 0.0;
    								for ($i3 = 0; $i3 < $n; $i3++)
    									$H[$i1][$i2] += $H2[$i1][$i3] * $H1[$i3][$i2];
    							}
    						}
    						for ($i1 = 0; $i1 < $n; $i1++) {
    							for ($i2 = 0; $i2 < $n; $i2++)
    								$H[$i1][$i2] += $s[$i1] * $s[$i2] / $sm1;
    						}
    					}
    					$y1 = $y2;
    					for ($i1 = 0; $i1 < $n; $i1++) {
    						$g1[$i1] = $g2[$i1];
    						$x1[$i1] = $x2[$i1];
    					}
    				}
    			}
    		}
    	}
    
    	if ($sw == 0)
    		$sw = -1;
    
    	return $sw;
    }
    
    /*********************************************/
    /* 与えられた点xから,dx方向へk*dxだけ進んだ */
    /* 点における関数値を計算する                */
    /*********************************************/
    					// 関数1
    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;
    }
    					// 関数2
    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;
    }
    					// 関数3
    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;
    }
    
    /********************/
    /* 微係数を計算する */
    /********************/
    					// 関数1
    function dsnx1($x, &$dx)
    {
    	$dx[0] = -2.0 * ($x[0] - 1.0);
    	$dx[1] = -2.0 * ($x[1] - 2.0);
    }
    					// 関数2
    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]);
    }
    					// 関数3
    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;
    }
    
    ?>
    			

  5. Ruby

    #*********************************/
    # 準Newton法によるの最小値の計算 */
    #      coded by Y.Suganuma       */
    #*********************************/
    
    #********************************************/
    # 与えられた点xから,dx方向へk*dxだけ進んだ */
    # 点における関数値及び微係数を計算する      */
    #********************************************/
    					# 関数1
    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
    }
    					# 関数2
    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
    }
    					# 関数3
    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
    
    #*******************************************************/
    # DFP法 or BFGS法                                      */
    #      method : =0 : DFP法                             */
    #               =1 : BFGS法                            */
    #      opt_1 : =0 : 1次元最適化を行わない             */
    #              =1 : 1次元最適化を行う                 */
    #      max : 最大繰り返し回数                          */
    #      n : 次元                                        */
    #      eps : 収束判定条件                              */
    #      step : きざみ幅                                 */
    #      y : 最小値                                      */
    #      x1 : x(初期値と答え)                            */
    #      dx : 関数の微分値                               */
    #      H : Hesse行列の逆行列                           */
    #      snx : 関数値とその微分を計算する関数名          */
    #            (微分は,その符号を変える)               */
    #      return : >=0 : 正常終了(収束回数)               */
    #               =-1 : 収束せず                         */
    #               =-2 : 1次元最適化の区間が求まらない   */
    #               =-3 : 黄金分割法が失敗                 */
    #*******************************************************/
    def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, h, &snx)
    
    	count = 0
    	sw    = 0
    	y2    = Array.new(1)
    	sw1   = Array.new(1)
    	x2    = Array.new(n)
    	g     = Array.new(n)
    	g1    = Array.new(n)
    	g2    = Array.new(n)
    	s     = Array.new(n)
    	w     = Array.new(n)
    
    	y1 = snx.call(0, 0.0, x1, dx)
    	snx.call(1, 0.0, x1, g1)
    
    	h1 = Array.new(n)
    	h2 = Array.new(n)
    	i  = Array.new(n)
    	for i1 in 0 ... n
    		h1[i1] = Array.new(n)
    		h2[i1] = Array.new(n)
    		i[i1]  = Array.new(n)
    		for i2 in 0 ... n
    			h[i1][i2] = 0.0
    			i[i1][i2] = 0.0
    		end
    		h[i1][i1] = 1.0
    		i[i1][i1] = 1.0
    	end
    
    	while count < max && sw == 0
    		count += 1
    					# 方向の計算
    		for i1 in 0 ... n
    			dx[i1] = 0.0
    			for i2 in 0 ... n
    				dx[i1] -= h[i1][i2] * g1[i2]
    			end
    		end
    					# 1次元最適化を行わない
    		if opt_1 == 0
    						# 新しい点
    			for i1 in 0 ... n
    				x2[i1] = x1[i1] + step * dx[i1]
    			end
    						# 新しい関数値
    			y2[0] = snx.call(0, 0.0, x2, dx)
    					# 1次元最適化を行う
    		else
    						# 区間を決める
    			sw1[0] = 0
    			f1     = y1
    			sp     = step
    			f2     = snx.call(0, sp, x1, dx)
    			if f2 > f1
    				sp = -step
    			end
    			for i1 in 0 ... max
    				f2 = snx.call(0, sp, x1, 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, x1, dx, &snx)
    							# 黄金分割法が失敗
    				if sw1[0] < 0
    					sw = -3
    							# 黄金分割法が成功
    				else
    								# 新しい点
    					for i1 in 0 ... n
    						x2[i1] = x1[i1] + k * dx[i1]
    					end
    				end
    			end
    		end
    
    		if sw == 0
    					# 収束(関数値の変化<eps)
    			if (y2[0]-y1).abs() < eps
    				sw   = count
    				y[0] = y2[0]
    				for i1 in 0 ... n
    					x1[i1] = x2[i1]
    				end
    					# 関数値の変化が大きい
    			else
    						# 傾きの計算
    				snx.call(1, 0.0, x2, g2)
    				sm = 0.0
    				for i1 in 0 ... n
    					sm += g2[i1] * g2[i1]
    				end
    				sm = Math.sqrt(sm)
    						# 収束(傾き<eps)
    				if sm < eps
    					sw   = count
    					y[0] = y2[0]
    					for i1 in 0 ... n
    						x1[i1] = x2[i1]
    					end
    						# 収束していない
    				else
    					for i1 in 0 ... n
    						g[i1] = g2[i1] - g1[i1]
    						s[i1] = x2[i1] - x1[i1]
    					end
    					sm1 = 0.0
    					for i1 in 0 ... n
    						sm1 += s[i1] * g[i1]
    					end
    							# DFP法
    					if method == 0
    						for i1 in 0 ... n
    							w[i1] = 0.0
    							for i2 in 0 ... n
    								w[i1] += g[i2] * h[i2][i1]
    							end
    						end
    						sm2 = 0.0
    						for i1 in 0 ... n
    							sm2 += w[i1] * g[i1]
    						end
    						for i1 in 0 ... n
    							w[i1] = 0.0
    							for i2 in 0 ... n
    								w[i1] += h[i1][i2] * g[i2]
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h1[i1][i2] = w[i1] * g[i2]
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h2[i1][i2] = 0.0
    								for i3 in 0 ... n
    									h2[i1][i2] += h1[i1][i3] * h[i3][i2]
    								end
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h[i1][i2] = h[i1][i2]  + s[i1] * s[i2] / sm1 - h2[i1][i2] / sm2
    							end
    						end
    							# BFGS法
    					else
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h1[i1][i2] = i[i1][i2] - s[i1] * g[i2] / sm1
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h2[i1][i2] = 0.0
    								for i3 in 0 ... n
    									h2[i1][i2] += h1[i1][i3] * h[i3][i2]
    								end
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h[i1][i2] = 0.0
    								for i3 in 0 ... n
    									h[i1][i2] += h2[i1][i3] * h1[i3][i2]
    								end
    							end
    						end
    						for i1 in 0 ... n
    							for i2 in 0 ... n
    								h[i1][i2] += s[i1] * s[i2] / sm1
    							end
    						end
    					end
    					y1 = y2[0]
    					for i1 in 0 ... n
    						g1[i1] = g2[i1]
    						x1[i1] = x2[i1]
    					end
    				end
    			end
    		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(" ");
    method = Integer(a[1])
    eps    = Float(a[3])
    step   = Float(a[5])
    x      = Array.new(n)
    dx     = Array.new(n)
    h      = 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
    	h[i1]  = Array.new(n)
    end
    				# 実行
    case fun
    	when 1
    		sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx1)
    	when 2
    		sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx2)
    	when 3
    		sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx3)
    end
    				# 結果の出力
    if sw < 0
    	printf("   収束しませんでした!")
    	case sw
    		when -1
    			printf("(収束回数)\n")
    		when -2
    			printf("(1次元最適化の区間)\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
    			

  6. Python

    # -*- 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
    
    ############################################
    # DFP法 or BFGS法
    #      method : =0 : DFP法
    #               =1 : BFGS法
    #      opt_1 : =0 : 1次元最適化を行わない
    #              =1 : 1次元最適化を行う
    #      max : 最大繰り返し回数
    #      n : 次元
    #      eps : 収束判定条件
    #      step : きざみ幅
    #      y : 最小値
    #      x1 : x(初期値と答え)
    #      dx : 関数の微分値
    #      H : Hesse行列の逆行列
    #      snx : 関数値を計算する関数名
    #      dsnx : 関数の微分を計算する関数名(符号を変える)
    #      return : >=0 : 正常終了(収束回数)
    #               =-1 : 収束せず
    #               =-2 : 1次元最適化の区間が求まらない
    #               =-3 : 黄金分割法が失敗
    #      coded by Y.Suganuma
    ############################################
    
    def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, H, snx, dsnx) :
    
    	count = 0
    	sw    = 0
    	sw1   = np.empty(1, np.int)
    	y2    = np.empty(1, np.float)
    	x2    = np.empty(n, np.float)
    	g     = np.empty(n, np.float)
    	g1    = np.empty(n, np.float)
    	g2    = np.empty(n, np.float)
    	s     = np.empty(n, np.float)
    	w     = np.empty(n, np.float)
    	I     = np.identity(n, np.float)
    	H1    = np.empty((n, n), np.float)
    	H2    = np.empty((n, n), np.float)
    
    	y1    = snx(0.0, x1, dx)
    	dsnx(x1, g1)
    
    	while count < max and sw == 0 :
    		count += 1
    			# 方向の計算
    		for i1 in range(0, n) :
    			dx[i1] = 0.0
    			for i2 in range(0, n) :
    				dx[i1] -= H[i1][i2] * g1[i2]
    			# 1次元最適化を行わない
    		if opt_1 == 0 :
    				# 新しい点
    			for i1 in range(0, n) :
    				x2[i1] = x1[i1] + step * dx[i1]
    				# 新しい関数値
    			y2[0] = snx(0.0, x2, dx)
    			# 1次元最適化を行う
    		else :
    				# 区間を決める
    			sw1[0] = 0
    			f1     = y1
    			sp     = step
    			f2     = snx(sp, x1, dx)
    			if f2 > f1 :
    				sp = -step
    			for i1 in range(0, max) :
    				f2 = snx(sp, x1, 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, x1, dx, snx)
    					# 黄金分割法が失敗
    				if sw1[0] < 0 :
    					sw = -3
    					# 黄金分割法が成功
    				else :
    						# 新しい点
    					for i1 in range(0, n) :
    						x2[i1] = x1[i1] + k * dx[i1]
    		if sw == 0 :
    			# 収束(関数値の変化<eps)
    			if abs(y2[0]-y1) < eps :
    				sw   = count
    				y[0] = y2[0]
    				for i1 in range(0, n) :
    					x1[i1] = x2[i1]
    			# 関数値の変化が大きい
    			else :
    				# 傾きの計算
    				dsnx(x2, g2)
    				sm = 0.0
    				for i1 in range(0, n) :
    					sm += g2[i1] * g2[i1]
    				sm = sqrt(sm)
    				# 収束(傾き<eps)
    				if sm < eps :
    					sw   = count
    					y[0] = y2[0]
    					for i1 in range(0, n) :
    						x1[i1] = x2[i1]
    				# 収束していない
    				else :
    					for i1 in range(0, n) :
    						g[i1] = g2[i1] - g1[i1]
    						s[i1] = x2[i1] - x1[i1]
    					sm1 = 0.0
    					for i1 in range(0, n) :
    						sm1 += s[i1] * g[i1]
    							# DFP法
    					if method == 0 :
    						for i1 in range(0, n) :
    							w[i1] = 0.0
    							for i2 in range(0, n) :
    								w[i1] += g[i2] * H[i2][i1]
    						sm2 = 0.0
    						for i1 in range(0, n) :
    							sm2 += w[i1] * g[i1]
    						for i1 in range(0, n) :
    							w[i1] = 0.0
    							for i2 in range(0, n) :
    								w[i1] += H[i1][i2] * g[i2]
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H1[i1][i2] = w[i1] * g[i2]
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H2[i1][i2] = 0.0
    								for i3 in range(0, n) :
    									H2[i1][i2] += H1[i1][i3] * H[i3][i2]
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2
    					# BFGS法
    					else :
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H2[i1][i2] = 0.0
    								for i3 in range(0, n) :
    									H2[i1][i2] += H1[i1][i3] * H[i3][i2]
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H[i1][i2] = 0.0
    								for i3 in range(0, n) :
    									H[i1][i2] += H2[i1][i3] * H1[i3][i2]
    						for i1 in range(0, n) :
    							for i2 in range(0, n) :
    								H[i1][i2] += s[i1] * s[i2] / sm1
    					y1 = y2[0]
    					for i1 in range(0, n) :
    						g1[i1] = g2[i1]
    						x1[i1] = x2[i1]
    
    	if sw == 0 :
    		sw = -1
    
    	return sw
    
    ############################################
    # 準Newton法によるの最小値の計算
    #      coded by Y.Suganuma
    ############################################
    
    ############################################
    # 与えられた点xから,dx方向へk*dxだけ進んだ
    # 点における関数値を計算する
    ############################################
    			 # 関数1
    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
    			# 関数2
    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
    			# 関数3
    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
    
    ############################################
    # 微係数を計算する
    ############################################
    			# 関数1
    def dsnx1(x, dx) :
    	dx[0] = -2.0 * (x[0] - 1.0)
    	dx[1] = -2.0 * (x[1] - 2.0)
    			# 関数2
    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])
    			# 関数3
    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()
    method = int(s[1])
    eps    = float(s[3])
    step   = float(s[5])
    
    x     = np.empty(n, np.float)
    dx    = np.empty(n, np.float)
    H     = np.identity(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 = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1)
    elif fun == 2 :
    	sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2)
    else :
    	sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3)
    			# 結果の出力
    if sw < 0 :
    	print("   収束しませんでした!", end="")
    	if sw == -1 :
    		print("(収束回数)")
    	elif sw == -2 :
    		print("(1次元最適化の区間)")
    	else :
    		print("(黄金分割法)")
    else :
    	print("   結果=", end="")
    	for i1 in range(0, n) :
    		print(str(x[i1]) + " ", end="")
    	print(" 最小値=" + str(y[0]) + "  回数=" + str(sw))
    			

  7. C#

    /**********************************/
    /* 準Newton法によるの最小値の計算 */
    /*      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);
    		int method  = int.Parse(str[1]);
    		double eps  = double.Parse(str[3]);
    		double step = double.Parse(str[5]);
    							// 3 行目
    		double[][] H = new double [n][];
    		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]);
    			H[i1] = new double [n];
    		}
    					// 実行
    		int sw;
    		if (fun == 1)
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx1, dsnx1);
    		else if (fun == 2)
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx2, dsnx2);
    		else
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step, ref y, x, dx, H, snx3, dsnx3);
    					// 結果の出力
    		if (sw < 0) {
    			Console.Write("   収束しませんでした!");
    			switch (sw) {
    				case -1:
    					Console.WriteLine("(収束回数)");
    					break;
    				case -2:
    					Console.WriteLine("(1次元最適化の区間)");
    					break;
    				case -3:
    					Console.WriteLine("(黄金分割法)");
    					break;
    			}
    		}
    		else {
    			Console.Write("   結果=");
    			for (int i1 = 0; i1 < n; i1++)
    				Console.Write(x[i1] + " ");
    			Console.WriteLine(" 最小値=" + y + "  回数=" + sw);
    		}
    	}
    			// 関数1の値の計算
    	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;
    	}
    			// 関数2の値の計算
    	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;
    	}
    			// 関数3の値の計算
    	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;
    	}
    					// 関数1の微分
    	int dsnx1(double[] x, double[] dx)
    	{
    		dx[0] = -2.0 * (x[0] - 1.0);
    		dx[1] = -2.0 * (x[1] - 2.0);
    		return 0;
    	}
    					// 関数2の微分
    	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;
    	}
    					// 関数3の微分
    	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;
    	}
    
    	/********************************************************/
    	/* DFP法 or BFGS法                                      */
    	/*      method : =0 : DFP法                             */
    	/*               =1 : BFGS法                            */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x1 : x(初期値と答え)                            */
    	/*      dx : 関数の微分値                               */
    	/*      H : Hesse行列の逆行列                           */
    	/*      fn : 関数値を計算する関数                       */
    	/*      dfn : 関数の微分値を計算する関数                */
    	/*      return : >=0 : 正常終了(収束回数)               */
    	/*               =-1 : 収束せず                         */
    	/*               =-2 : 1次元最適化の区間が求まらない   */
    	/*               =-3 : 黄金分割法が失敗                 */
    	/********************************************************/
    	int DFP_BFGS(int method, int opt_1, int max, int n, double eps, double step,
    	             ref double y, double[] x1, double[] dx, double[][] H,
    	           Func<double, double[], double[], double> fn,
    	           Func<double[], double[], int> dfn)
    	{
    		double[] x2 = new double [n];
    		double[] g  = new double [n];
    		double[] g1 = new double [n];
    		double[] g2 = new double [n];
    		double[] s  = new double [n];
    		double[] w  = new double [n];
    
    		double y1 = fn(0.0, x1, dx);
    		double y2 = 0.0;
    		dfn(x1, g1);
    
    		double[][] H1 = new double [n][];
    		double[][] H2 = new double [n][];
    		double[][] I  = new double [n][];
    		for (int i1 = 0; i1 < n; i1++) {
    			H1[i1] = new double [n];
    			H2[i1] = new double [n];
    			I[i1]  = new double [n];
    			for (int i2 = 0; i2 < n; i2++) {
    				H[i1][i2] = 0.0;
    				I[i1][i2] = 0.0;
    			}
    			H[i1][i1] = 1.0;
    			I[i1][i1] = 1.0;
    		}
    
    		int count = 0;
    		int sw    = 0;
    
    		while (count < max && sw == 0) {
    			count++;
    					// 方向の計算
    			for (int i1 = 0; i1 < n; i1++) {
    				dx[i1] = 0.0;
    				for (int i2 = 0; i2 < n; i2++)
    					dx[i1] -= H[i1][i2] * g1[i2];
    			}
    					// 1次元最適化を行わない
    			if (opt_1 == 0) {
    						// 新しい点
    				for (int i1 = 0; i1 < n; i1++)
    					x2[i1] = x1[i1] + step * dx[i1];
    						// 新しい関数値
    				y2 = fn(0.0, x2, dx);
    			}
    					// 1次元最適化を行う
    			else {
    						// 区間を決める
    				int sw1   = 0;
    				double f1 = y1;
    				double sp = step;
    				double f2 = fn(sp, x1, dx);
    				if (f2 > f1)
    					sp = -step;
    				for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
    					f2 = fn(sp, x1, 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, x1, dx, fn);
    							// 黄金分割法が失敗
    					if (sw1 < 0)
    						sw = -3;
    							// 黄金分割法が成功
    					else {
    								// 新しい点
    						for (int i1 = 0; i1 < n; i1++)
    							x2[i1] = x1[i1] + k * dx[i1];
    					}
    				}
    			}
    
    			if (sw == 0) {
    					// 収束(関数値の変化<eps)
    				if (Math.Abs(y2-y1) < eps) {
    					sw = count;
    					y  = y2;
    					for (int i1 = 0; i1 < n; i1++)
    						x1[i1] = x2[i1];
    				}
    					// 関数値の変化が大きい
    				else {
    						// 傾きの計算
    					dfn(x2, g2);
    					double sm = 0.0;
    					for (int i1 = 0; i1 < n; i1++)
    						sm += g2[i1] * g2[i1];
    					sm = Math.Sqrt(sm);
    						// 収束(傾き<eps)
    					if (sm < eps) {
    						sw = count;
    						y  = y2;
    						for (int i1 = 0; i1 < n; i1++)
    							x1[i1] = x2[i1];
    					}
    						// 収束していない
    					else {
    						for (int i1 = 0; i1 < n; i1++) {
    							g[i1] = g2[i1] - g1[i1];
    							s[i1] = x2[i1] - x1[i1];
    						}
    						double sm1 = 0.0;
    						for (int i1 = 0; i1 < n; i1++)
    							sm1 += s[i1] * g[i1];
    							// DFP法
    						if (method == 0) {
    							for (int i1 = 0; i1 < n; i1++) {
    								w[i1] = 0.0;
    								for (int i2 = 0; i2 < n; i2++)
    									w[i1] += g[i2] * H[i2][i1];
    							}
    							double sm2 = 0.0;
    							for (int i1 = 0; i1 < n; i1++)
    								sm2 += w[i1] * g[i1];
    							for (int i1 = 0; i1 < n; i1++) {
    								w[i1] = 0.0;
    								for (int i2 = 0; i2 < n; i2++)
    									w[i1] += H[i1][i2] * g[i2];
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++)
    									H1[i1][i2] = w[i1] * g[i2];
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++) {
    									H2[i1][i2] = 0.0;
    									for (int i3 = 0; i3 < n; i3++)
    										H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    								}
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++)
    									H[i1][i2] = H[i1][i2]  + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2;
    							}
    						}
    							// BFGS法
    						else {
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++)
    									H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1;
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++) {
    									H2[i1][i2] = 0.0;
    									for (int i3 = 0; i3 < n; i3++)
    										H2[i1][i2] += H1[i1][i3] * H[i3][i2];
    								}
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++) {
    									H[i1][i2] = 0.0;
    									for (int i3 = 0; i3 < n; i3++)
    										H[i1][i2] += H2[i1][i3] * H1[i3][i2];
    								}
    							}
    							for (int i1 = 0; i1 < n; i1++) {
    								for (int i2 = 0; i2 < n; i2++)
    									H[i1][i2] += s[i1] * s[i2] / sm1;
    							}
    						}
    						y1 = y2;
    						for (int i1 = 0; i1 < n; i1++) {
    							g1[i1] = g2[i1];
    							x1[i1] = x2[i1];
    						}
    					}
    				}
    			}
    		}
    
    		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;
    	}
    }
    			

  8. VB

    ''''''''''''''''''''''''''''''''''
    ' 準Newton法によるの最小値の計算 '
    '      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 method As Integer = Integer.Parse(str(1))
    		Dim eps As Double     = Double.Parse(str(3))
    		Dim step1 As Double   = Double.Parse(str(5))
    							' 3 行目
    		Dim H(n,n) As Double
    		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
    					' 関数1の値の計算(ラムダ式)
    		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
    					' 関数2の値の計算(ラムダ式)
    		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
    					' 関数3の値の計算(ラムダ式)
    		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
    					' 関数1の微分(ラムダ式)
    		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
    					' 関数2の微分(ラムダ式)
    		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
    					' 関数3の微分(ラムダ式)
    		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 = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx1, dsnx1)
    		ElseIf fun = 2
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx2, dsnx2)
    		Else
    			sw = DFP_BFGS(method, opt_1, max, n, eps, step1, y, x, dx, H, snx3, dsnx3)
    		End If
    
    					' 結果の出力
    		If sw < 0
    			Console.Write("   収束しませんでした!")
    			If sw = -1
    				Console.WriteLine("(収束回数)")
    			ElseIf sw = -2
    				Console.WriteLine("(1次元最適化の区間)")
    			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
    
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' DFP法 or BFGS法                                      '
    	'      method : =0 : DFP法                             '
    	'               =1 : BFGS法                            '
    	'      opt_1 : =0 : 1次元最適化を行わない             '
    	'              =1 : 1次元最適化を行う                 '
    	'      max : 最大繰り返し回数                          '
    	'      n : 次元                                        '
    	'      eps : 収束判定条件                              '
    	'      step1 : きざみ幅                                '
    	'      y : 最小値                                      '
    	'      x1 : x(初期値と答え)                            '
    	'      dx : 関数の微分値                               '
    	'      H : Hesse行列の逆行列                           '
    	'      fn : 関数値を計算する関数                       '
    	'      dfn : 関数の微分値を計算する関数                '
    	'      return : >=0 : 正常終了(収束回数)               '
    	'               =-1 : 収束せず                         '
    	'               =-2 : 1次元最適化の区間が求まらない   '
    	'               =-3 : 黄金分割法が失敗                 '
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function DFP_BFGS(method As Integer, opt_1 As Integer, max As Integer, n As Integer,
    	                  eps As Double, step1 As Double, ByRef y As Double, x1() As Double,
    	                  dx() As Double, H(,) As Double,
    	                  fn As Func(Of Double, Double(), Double(), Double),
    	                  dfn As Func(Of Double(), Double(), Integer))
    
    		Dim x2(n) As Double
    		Dim g(n) As Double
    		Dim g1(n) As Double
    		Dim g2(n) As Double
    		Dim s(n) As Double
    		Dim w(n) As Double
    
    		Dim y1 As Double = fn(0.0, x1, dx)
    		Dim y2 As Double = 0.0
    		dfn(x1, g1)
    
    		Dim H1(n,n) As Double
    		Dim H2(n,n) As Double
    		Dim I(n,n) As Double
    		For i1 As Integer = 0 To n-1
    			For i2 As Integer = 0 To n-1
    				H(i1,i2) = 0.0
    				I(i1,i2) = 0.0
    			Next
    			H(i1,i1) = 1.0
    			I(i1,i1) = 1.0
    		Next
    
    		Dim count As Integer = 0
    		Dim sw As Integer    = 0
    
    		Do While count < max and sw = 0
    			count += 1
    					' 方向の計算
    			For i1 As Integer = 0 To n-1
    				dx(i1) = 0.0
    				For i2 As Integer = 0 To n-1
    					dx(i1) -= H(i1,i2) * g1(i2)
    				Next
    			Next
    					' 1次元最適化を行わない
    			If opt_1 = 0
    						' 新しい点
    				For i1 As Integer = 0 To n-1
    					x2(i1) = x1(i1) + step1 * dx(i1)
    				Next
    						' 新しい関数値
    				y2 = fn(0.0, x2, dx)
    					' 1次元最適化を行う
    			Else
    						' 区間を決める
    				Dim sw1 As Integer = 0
    				Dim f1 As Double   = y1
    				Dim sp As Double   = step1
    				Dim f2 As Double   = fn(sp, x1, dx)
    				If f2 > f1
    					sp = -step1
    				End If
    				Dim ii As Integer = 0
    				Do While ii < max and sw1 = 0
    					f2 = fn(sp, x1, 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, x1, dx, fn)
    							' 黄金分割法が失敗
    					If sw1 < 0
    						sw = -3
    							' 黄金分割法が成功
    					Else
    								' 新しい点
    						For i1 As Integer = 0 To n-1
    							x2(i1) = x1(i1) + k * dx(i1)
    						Next
    					End If
    				End If
    			End If
    
    			If sw = 0
    					' 収束(関数値の変化<eps)
    				If Math.Abs(y2-y1) < eps
    					sw = count
    					y  = y2
    					For i1 As Integer = 0 To n-1
    						x1(i1) = x2(i1)
    					Next
    					' 関数値の変化が大きい
    				Else
    						' 傾きの計算
    					dfn(x2, g2)
    					Dim sm As Double = 0.0
    					For i1 As Integer = 0 To n-1
    						sm += g2(i1) * g2(i1)
    					Next
    					sm = Math.Sqrt(sm)
    						' 収束(傾き<eps)
    					If sm < eps
    						sw = count
    						y  = y2
    						For i1 As Integer = 0 To n-1
    							x1(i1) = x2(i1)
    						Next
    						' 収束していない
    					Else
    						For i1 As Integer = 0 To n-1
    							g(i1) = g2(i1) - g1(i1)
    							s(i1) = x2(i1) - x1(i1)
    						Next
    						Dim sm1 As Double = 0.0
    						For i1 As Integer = 0 To n-1
    							sm1 += s(i1) * g(i1)
    						Next
    							' DFP法
    						If method = 0
    							For i1 As Integer = 0 To n-1
    								w(i1) = 0.0
    								For i2 As Integer = 0 To n-1
    									w(i1) += g(i2) * H(i2,i1)
    								Next
    							Next
    							Dim sm2 As Double = 0.0
    							For i1 As Integer = 0 To n-1
    								sm2 += w(i1) * g(i1)
    							Next
    							For i1 As Integer = 0 To n-1
    								w(i1) = 0.0
    								For i2 As Integer = 0 To n-1
    									w(i1) += H(i1,i2) * g(i2)
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H1(i1,i2) = w(i1) * g(i2)
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H2(i1,i2) = 0.0
    									For i3 As Integer = 0 To n-1
    										H2(i1,i2) += H1(i1,i3) * H(i3,i2)
    									Next
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H(i1,i2) = H(i1,i2)  + s(i1) * s(i2) / sm1 - H2(i1,i2) / sm2
    								Next
    							Next
    							' BFGS法
    						Else
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H1(i1,i2) = I(i1,i2) - s(i1) * g(i2) / sm1
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H2(i1,i2) = 0.0
    									For i3 As Integer = 0 To n-1
    										H2(i1,i2) += H1(i1,i3) * H(i3,i2)
    									Next
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H(i1,i2) = 0.0
    									For i3 As Integer = 0 To n-1
    										H(i1,i2) += H2(i1,i3) * H1(i3,i2)
    									Next
    								Next
    							Next
    							For i1 As Integer = 0 To n-1
    								For i2 As Integer = 0 To n-1
    									H(i1,i2) += s(i1) * s(i2) / sm1
    								Next
    							Next
    						End If
    						y1 = y2
    						For i1 As Integer = 0 To n-1
    							g1(i1) = g2(i1)
    							x1(i1) = x2(i1)
    						Next
    					End If
    				End If
    			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
    			

静岡理工科大学 菅沼ホーム 目次 索引