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

最適化( 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 *);
    int hesse1(double *, double **, double);
    void dsnx2(double *, double *);
    double snx2(double, double *, double *);
    int hesse2(double *, double **, double);
    void dsnx3(double *, double *);
    double snx3(double, double *, double *);
    int hesse3(double *, double **, double);
    int Newton(int, int, int, double, double, double *, double *, double *, double **,
              double (*)(double, double *, double *),
              void (*)(double *, double *),
              int (*)(double *, double **, double));
    
    int main()
    {
    	double eps, **H, step, *x, *dx, y;
    	int fun, i1, max, n, opt_1, sw = 0;
    					// データの入力
    	scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
    	scanf("%*s %lf %*s %lf", &eps, &step);
    	x  = new double [n];
    	dx = new double [n];
    	H  = new double * [n];
    	scanf("%*s");
    	for (i1 = 0; i1 < n; i1++) {
    		scanf("%lf", &x[i1]);
    		H[i1] = new double [2*n];
    	}
    					// 実行
    	switch (fun) {
    		case 1:
    			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx1, dsnx1, hesse1);
    			break;
    		case 2:
    			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx2, dsnx2, hesse2);
    			break;
    		case 3:
    			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx3, dsnx3, hesse3);
    			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;
    }
    
    /********************************************************/
    /* Newton法                                             */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x : x(初期値と答え)                             */
    /*      dx : 関数の微分値                               */
    /*      H : Hesse行列の逆行列                           */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      hesse : Hesse行列の逆行列を計算する関数名       */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    #include <math.h>
    
    double gold(double, double, double, double *, int *, int, double *, double *,
                double (*)(double, double *, double *));
    
    int Newton(int opt_1, int max, int n, double eps, double step, double *y,
              double *x, double *dx, double **H, double (*snx)(double, double *, double *), 
              void (*dsnx)(double *, double *), int (*hesse)(double *, double **, double))
    {
    	double f1, f2, k, sp, y1, y2;
    	double *wk = new double [n];
    	int count = 0, i1, i2, sw = 0, sw1;
    
    	y1 = snx(0.0, x, dx);
    
    	while (count < max && sw == 0) {
    					// 傾きの計算
    		dsnx(x, wk);
    					// Hesse行列の逆行列の計算
    		sw1 = hesse(x, H, eps);
    					// 収束していない
    		if (sw1 == 0) {
    						// 方向の計算
    			count++;
    			for (i1 = 0; i1 < n; i1++) {
    				dx[i1] = 0.0;
    				for (i2 = 0; i2 < n; i2++)
    					dx[i1] += H[i1][n+i2] * wk[i2];
    			}
    						// 1次元最適化を行わない
    			if (opt_1 == 0) {
    							// 新しい点
    				for (i1 = 0; i1 < n; i1++)
    					x[i1] += dx[i1];
    							// 新しい関数値
    				y2 = snx(0.0, x, dx);
    								// 関数値の変化が大きい
    				if (fabs(y2-y1) > eps)
    					y1 = y2;
    								// 収束(関数値の変化<eps)
    				else {
    					sw = count;
    					*y = y2;
    				}
    			}
    						// 1次元最適化を行う
    			else {
    							// 区間を決める
    				sw1 = 0;
    				f1  = y1;
    				sp  = step;
    				f2  = snx(sp, x, dx);
    				if (f2 > f1)
    					sp = -step;
    				for (i1 = 0; i1 < max && sw1 == 0; i1++) {
    					f2 = snx(sp, x, dx);
    					if (f2 > f1)
    						sw1 = 1;
    					else {
    						sp *= 2.0;
    						f1  = f2;
    					}
    				}
    							// 区間が求まらない
    				if (sw1 == 0)
    					sw = -2;
    							// 区間が求まった
    				else {
    								// 黄金分割法
    					k = gold(0.0, sp, eps, &y2, &sw1, max, x, dx, fn);
    								// 黄金分割法が失敗
    					if (sw1 < 0)
    						sw = -3;
    								// 黄金分割法が成功
    					else {
    									// 新しい点
    						for (i1 = 0; i1 < n; i1++)
    							x[i1] += k * dx[i1];
    										// 関数値の変化が大きい
    						if (fabs(y1-y2) > eps)
    							y1 = y2;
    										// 収束(関数値の変化<eps)
    						else {
    							sw = count;
    							*y = y2;
    						}
    					}
    				}
    			}
    		}
    					// 収束(傾き<eps)
    		else {
    			sw = count;
    			*y = y1;
    		}
    	}
    
    	if (sw == 0)
    		sw = -1;
    
    	delete [] wk;
    
    	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]));
    }
    
    /*******************************/
    /* Hesse行列の逆行列を計算する */
    /*******************************/
    int Gauss(double **, int, int, double);
    					// 関数1
    int hesse1(double *x, double **H, double eps)
    {
    	int sw;
    
    	H[0][0] = 2.0;
    	H[0][1] = 0.0;
    	H[1][0] = 0.0;
    	H[1][1] = 2.0;
    	H[0][2] = 1.0;
    	H[0][3] = 0.0;
    	H[1][2] = 0.0;
    	H[1][3] = 1.0;
    
    	sw = Gauss(H, 2, 2, eps);
    
    	return sw;
    }
    					// 関数2
    int hesse2(double *x, double **H, double eps)
    {
    	int sw;
    
    	H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
    	H[0][1] = -400.0 * x[0];
    	H[1][0] = H[0][1];
    	H[1][1] = 200.0;
    	H[0][2] = 1.0;
    	H[0][3] = 0.0;
    	H[1][2] = 0.0;
    	H[1][3] = 1.0;
    
    	sw = Gauss(H, 2, 2, eps);
    
    	return sw;
    }
    					// 関数3
    int hesse3(double *x, double **H, double eps)
    {
    	double x1, x2, x3;
    	int sw;
    
    	x1      = 1.0 - x[1];
    	x2      = 1.0 - x[1] * x[1];
    	x3      = 1.0 - x[1] * x[1] * x[1];
    	H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
    	H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                  4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                  6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
    	H[1][0] = H[0][1];
    	H[1][1] = 2.0 * x[0] * x[0] +
                  4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                  12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                  18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
    	H[0][2] = 1.0;
    	H[0][3] = 0.0;
    	H[1][2] = 0.0;
    	H[1][3] = 1.0;
    
    	sw = Gauss(H, 2, 2, eps);
    
    	return sw;
    }
    
    /*******************************************************/
    /* 線形連立方程式を解く(逆行列を求める)              */
    /*      w : 方程式の左辺及び右辺                       */
    /*      n : 方程式の数                                 */
    /*      m : 方程式の右辺の列の数                       */
    /*      eps : 正則性を判定する規準                     */
    /*      return : =0 : 正常                             */
    /*               =1 : 逆行列が存在しない               */
    /*******************************************************/
    #include <math.h>
    
    int Gauss(double **w, int n, int m, double eps)
    {
    	double y1, y2;
    	int ind = 0, nm, m1, m2, i1, i2, i3;
    
    	nm = n + m;
    
    	for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    		y1 = .0;
    		m1 = i1 + 1;
    		m2 = 0;
    
    		for (i2 = i1; i2 < n; i2++) {
    			y2 = fabs(w[i2][i1]);
    			if (y1 < y2) {
    				y1 = y2;
    				m2 = i2;
    			}
    		}
    
    		if (y1 < eps)
    			ind = 1;
    
    		else {
    
    			for (i2 = i1; i2 < nm; i2++) {
    				y1        = w[i1][i2];
    				w[i1][i2] = w[m2][i2];
    				w[m2][i2] = y1;
    			}
    
    			y1 = 1.0 / w[i1][i1];
    
    			for (i2 = m1; i2 < nm; i2++)
    				w[i1][i2] *= y1;
    
    			for (i2 = 0; i2 < n; i2++) {
    				if (i2 != i1) {
    					for (i3 = m1; i3 < nm; i3++)
    						w[i2][i3] -= w[i2][i1] * w[i1][i3];
    				}
    			}
    		}
    	}
    
    	return(ind);
    }
    			

  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, n, opt_1, sw = 0;
    		StringTokenizer str;
    		String line;
    		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
    					// データの入力
    							// 1 行目
    		line = in.readLine();
    		str  = new StringTokenizer(line, " ");
    		line = str.nextToken();
    		line = str.nextToken();
    		fun = Integer.parseInt(line);
    		line = str.nextToken();
    		line = str.nextToken();
    		n = Integer.parseInt(line);
    		line = str.nextToken();
    		line = str.nextToken();
    		max = Integer.parseInt(line);
    		line = str.nextToken();
    		line = str.nextToken();
    		opt_1 = Integer.parseInt(line);
    							// 2 行目
    		line = in.readLine();
    		str  = new StringTokenizer(line, " ");
    		line = str.nextToken();
    		line = str.nextToken();
    		eps = Double.parseDouble(line);
    		line = str.nextToken();
    		line = str.nextToken();
    		step = Double.parseDouble(line);
    							// 3 行目
    		x  = new double [n];
    		dx = new double [n];
    		y  = new double [1];
    		H  = new double [n][2*n];
    		line = in.readLine();
    		str  = new StringTokenizer(line, " ");
    		line = str.nextToken();
    		for (i1 = 0; i1 < n; i1++) {
    			line = str.nextToken();
    			x[i1] = Double.parseDouble(line);
    		}
    					// 実行
    		Kansu kn = new Kansu(fun);
    		sw = kn.Newton(opt_1, max, n, eps, step, y, x, dx, H);
    					// 結果の出力
    		if (sw < 0) {
    			System.out.print("   収束しませんでした!");
    			switch (sw) {
    				case -1:
    					System.out.println("(収束回数)");
    					break;
    				case -2:
    					System.out.print("(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 Newton
    {
    	private int sw;
    					// コンストラクタ
    	Kansu (int s) {sw = s;}
    					// 与えられた点xから,dx方向へk*dxだけ進んだ
    					// 点における関数値を計算する
    	double snx(double k, double x[], double dx[])
    	{
    		double x1, y1, z1, y = 0.0, w[] = new double [2];
    
    		switch (sw) {
    			case 1:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = w[0] - 1.0;
    				y1   = w[1] - 2.0;
    				y    = x1 * x1 + y1 * y1;
    				break;
    			case 2:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = w[1] - w[0] * w[0];
    				y1   = 1.0 - w[0];
    				y    = 100.0 * x1 * x1 + y1 * y1;
    				break;
    			case 3:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = 1.5 - w[0] * (1.0 - w[1]);
    				y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
    				z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
    				y    = x1 * x1 + y1 * y1 + z1 * z1;
    				break;
    		}
    
    		return y;
    	}
    					// 関数の微分
    	void dsnx(double x[], double dx[])
    	{
    		switch (sw) {
    			case 1:
    				dx[0] = -2.0 * (x[0] - 1.0);
    				dx[1] = -2.0 * (x[1] - 2.0);
    				break;
    			case 2:
    				dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
    				dx[1] = -200.0 * (x[1] - x[0] * x[0]);
    				break;
    			case 3:
    				dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                            2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                            2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    				dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                            4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                            6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    				break;
    		}
    	}
    					// Hesse行列の逆行列
    	int hesse(double x[], double H[][], double eps)
    	{
    		double x1, x2, x3;
    		int sw1 = 0;
    
    		switch (sw) {
    			case 1:
    				H[0][0] = 2.0;
    				H[0][1] = 0.0;
    				H[1][0] = 0.0;
    				H[1][1] = 2.0;
    				H[0][2] = 1.0;
    				H[0][3] = 0.0;
    				H[1][2] = 0.0;
    				H[1][3] = 1.0;
    				sw1 = gauss(H, 2, 2, eps);
    				break;
    			case 2:
    				H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
    				H[0][1] = -400.0 * x[0];
    				H[1][0] = H[0][1];
    				H[1][1] = 200.0;
    				H[0][2] = 1.0;
    				H[0][3] = 0.0;
    				H[1][2] = 0.0;
    				H[1][3] = 1.0;
    				sw1 = gauss(H, 2, 2, eps);
    				break;
    			case 3:
    				x1      = 1.0 - x[1];
    				x2      = 1.0 - x[1] * x[1];
    				x3      = 1.0 - x[1] * x[1] * x[1];
    				H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
    				H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                              4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                              6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
    				H[1][0] = H[0][1];
    				H[1][1] = 2.0 * x[0] * x[0] +
                              4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                              12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                              18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
    				H[0][2] = 1.0;
    				H[0][3] = 0.0;
    				H[1][2] = 0.0;
    				H[1][3] = 1.0;
    				sw1 = gauss(H, 2, 2, eps);
    				break;
    		}
    
    		return sw1;
    	}
    }
    
    abstract class Newton {
    
    	/********************************************************/
    	/* Newton法                                             */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x : 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[]);
    	abstract int hesse(double x[], double H[][], double eps);
    
    	int Newton(int opt_1, int max, int n, double eps, double step, double y[], double x[], double dx[], double H[][])
    	{
    		double f1, f2, k, sp, y1[] = new double [1], y2[] = new double [1];
    		double wk[] = new double [n];
    		int count = 0, i1, i2, sw = 0, sw1[] = new int [1];
    
    		y1[0] = snx(0.0, x, dx);
    
    		while (count < max && sw == 0) {
    					// 傾きの計算
    			dsnx(x, wk);
    					// Hesse行列の逆行列の計算
    			sw1[0] = hesse(x, H, eps);
    					// 収束していない
    			if (sw1[0] == 0) {
    						// 方向の計算
    				count++;
    				for (i1 = 0; i1 < n; i1++) {
    					dx[i1] = 0.0;
    					for (i2 = 0; i2 < n; i2++)
    						dx[i1] += H[i1][n+i2] * wk[i2];
    				}
    						// 1次元最適化を行わない
    				if (opt_1 == 0) {
    							// 新しい点
    					for (i1 = 0; i1 < n; i1++)
    						x[i1] += dx[i1];
    							// 新しい関数値
    					y2[0] = snx(0.0, x, dx);
    								// 関数値の変化が大きい
    					if (Math.abs(y2[0]-y1[0]) > eps)
    						y1[0] = y2[0];
    								// 収束(関数値の変化<eps)
    					else {
    						sw   = count;
    						y[0] = y2[0];
    					}
    				}
    						// 1次元最適化を行う
    				else {
    							// 区間を決める
    					sw1[0] = 0;
    					f1     = y1[0];
    					sp     = step;
    					f2     = snx(sp, x, dx);
    					if (f2 > f1)
    						sp = -step;
    					for (i1 = 0; i1 < max && sw1[0] == 0; i1++) {
    						f2 = snx(sp, x, dx);
    						if (f2 > f1)
    							sw1[0] = 1;
    						else {
    							sp *= 2.0;
    							f1  = f2;
    						}
    					}
    							// 区間が求まらない
    					if (sw1[0] == 0)
    						sw = -2;
    							// 区間が求まった
    					else {
    								// 黄金分割法
    						k = gold(0.0, sp, eps, y2, sw1, max, x, dx);
    								// 黄金分割法が失敗
    						if (sw1[0] < 0)
    							sw = -3;
    								// 黄金分割法が成功
    						else {
    									// 新しい点
    							for (i1 = 0; i1 < n; i1++)
    								x[i1] += k * dx[i1];
    										// 関数値の変化が大きい
    							if (Math.abs(y1[0]-y2[0]) > eps)
    								y1[0] = y2[0];
    										// 収束(関数値の変化<eps)
    							else {
    								sw   = count;
    								y[0] = y2[0];
    							}
    						}
    					}
    				}
    			}
    					// 収束(傾き<eps)
    			else {
    				sw   = count;
    				y[0] = y1[0];
    			}
    		}
    
    		if (sw == 0)
    			sw = -1;
    
    		return sw;
    	}
    
    	/****************************************************************/
    	/* 黄金分割法(与えられた方向での最小値)                         */
    	/*      a,b : 初期区間 a < b                                    */
    	/*      eps : 許容誤差                                          */
    	/*      val : 関数値                                            */
    	/*      ind : 計算状況                                          */
    	/*              >= 0 : 正常終了(収束回数)                       */
    	/*              = -1 : 収束せず                                 */
    	/*      max : 最大試行回数                                      */
    	/*      w : 位置                                                */
    	/*      dw : 傾きの成分                                         */
    	/*      return : 結果(w+y*dwのy)                              */
    	/****************************************************************/
    	double gold(double a, double b, double eps, double val[], int ind[], int max, double w[], double dw[])
    	{
    		double f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
    		int count = 0;
    					// 初期設定
    		ind[0] = -1;
    		x1     = b - tau * (b - a);
    		x2     = a + tau * (b - a);
    		f1     = snx(x1, w, dw);
    		f2     = snx(x2, w, dw);
    					// 計算
    		while (count < max && ind[0] < 0) {
    			count += 1;
    			if (f2 > f1) {
    				if (Math.abs(b-a) < eps) {
    					ind[0] = 0;
    					x      = x1;
    					val[0] = f1;
    				}
    				else {
    					b  = x2;
    					x2 = x1;
    					x1 = a + (1.0 - tau) * (b - a);
    					f2 = f1;
    					f1 = snx(x1, w, dw);
    				}
    			}
    			else {
    				if (Math.abs(b-a) < eps) {
    					ind[0] = 0;
    					x      = x2;
    					val[0] = f2;
    					f1     = f2;
    				}
    				else {
    					a  = x1;
    					x1 = x2;
    					x2 = b - (1.0 - tau) * (b - a);
    					f1 = f2;
    					f2 = snx(x2, w, dw);
    				}
    			}
    		}
    					// 収束した場合の処理
    		if (ind[0] == 0) {
    			ind[0] = count;
    			fa     = snx(a, w, dw);
    			fb     = snx(b, w, dw);
    			if (fb < fa) {
    				a  = b;
    				fa = fb;
    			}
    			if (fa < f1) {
    				x      = a;
    				val[0] = fa;
    			}
    		}
    
    		return x;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 正則性を判定する規準                     */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	int gauss(double w[][], int n, int m, double eps)
    	{
    		double y1, y2;
    		int ind = 0, nm, m1, m2, i1, i2, i3;
    
    		nm = n + m;
    
    		for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    			y1 = .0;
    			m1 = i1 + 1;
    			m2 = 0;
    
    			for (i2 = i1; i2 < n; i2++) {
    				y2 = Math.abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    
    			if (y1 < eps)
    				ind = 1;
    
    			else {
    
    				for (i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    
    				y1 = 1.0 / w[i1][i1];
    
    				for (i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return(ind);
    	}
    }
    			

  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();
    		str3 = 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;
    			str1 = document.getElementById("func").value + ";";
    			str2 = document.getElementById("dfunc").value.split(/\n{1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				str2[i1] = str2[i1] + ";";
    			str3 = document.getElementById("hesse").value.split(/\n{1,}/)
    			for (let i1 = 0; i1 < n*n; i1++)
    				str3[i1] = str3[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 = newton(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]);
    		}
    
    		/******************************************************/
    		/* Newton法                                           */
    		/*      opt_1 : =0 : 1次元最適化を行わない           */
    		/*              =1 : 1次元最適化を行う               */
    		/*      max : 最大繰り返し回数                        */
    		/*      n : 次元                                      */
    		/*      eps : 収束判定条件                            */
    		/*      step : きざみ幅                               */
    		/*      y : 最小値                                    */
    		/*      x : x(初期値と答え)                           */
    		/*      dx : 関数の微分値                             */
    		/*      H : Hesse行列の逆行列                         */
    		/*      fn : 関数値を計算する関数                     */
    		/*      dfn : 関数の微分値を計算する関数              */
    		/*      return : >=0 : 正常終了(収束回数)             */
    		/*               =-1 : 収束せず                       */
    		/*               =-2 : 1次元最適化の区間が求まらない */
    		/*               =-3 : 黄金分割法が失敗               */
    		/******************************************************/
    		function newton(opt_1, max, n, eps, step, y, x, dx, H, fn, dfn)
    		{
    			let f1, f2, k, sp, count = 0, sw = 0;
    			let y1 = new Array();
    			let y2 = new Array();
    			let wk = new Array();
    			let sw1 = new Array();
    
    			y1[0] = fn(0.0, x, dx);
    
    			while (count < max && sw == 0) {
    					// 傾きの計算
    				dfn(n, x, wk);
    					// Hesse行列の逆行列の計算
    				sw1[0] = hesse(n, x, H, eps);
    					// 収束していない
    				if (sw1[0] == 0) {
    						// 方向の計算
    					count++;
    					for (let i1 = 0; i1 < n; i1++) {
    						dx[i1] = 0.0;
    						for (let i2 = 0; i2 < n; i2++)
    							dx[i1] += H[i1][n+i2] * wk[i2];
    					}
    						// 1次元最適化を行わない
    					if (opt_1 == 0) {
    							// 新しい点
    						for (let i1 = 0; i1 < n; i1++)
    							x[i1] += dx[i1];
    							// 新しい関数値
    						y2[0] = fn(0.0, x, dx);
    								// 関数値の変化が大きい
    						if (Math.abs(y2[0]-y1[0]) > eps)
    							y1[0] = y2[0];
    								// 収束(関数値の変化<eps)
    						else {
    							sw   = count;
    							y[0] = y2[0];
    						}
    					}
    						// 1次元最適化を行う
    					else {
    							// 区間を決める
    						sw1[0] = 0;
    						f1     = y1[0];
    						sp     = step;
    						f2     = fn(sp, x, dx);
    						if (f2 > f1)
    							sp = -step;
    						for (let i1 = 0; i1 < max && sw1[0] == 0; i1++) {
    							f2 = fn(sp, x, dx);
    							if (f2 > f1)
    								sw1[0] = 1;
    							else {
    								sp *= 2.0;
    								f1  = f2;
    							}
    						}
    							// 区間が求まらない
    						if (sw1[0] == 0)
    							sw = -2;
    							// 区間が求まった
    						else {
    								// 黄金分割法
    							k = gold(0.0, sp, eps, y2, sw1, max, x, dx, snx);
    								// 黄金分割法が失敗
    							if (sw1[0] < 0)
    								sw = -3;
    								// 黄金分割法が成功
    							else {
    									// 新しい点
    								for (let i1 = 0; i1 < n; i1++)
    									x[i1] += k * dx[i1];
    										// 関数値の変化が大きい
    								if (Math.abs(y1[0]-y2[0]) > eps)
    									y1[0] = y2[0];
    										// 収束(関数値の変化<eps)
    								else {
    									sw   = count;
    									y[0] = y2[0];
    								}
    							}
    						}
    					}
    				}
    					// 収束(傾き<eps)
    				else {
    					sw   = count;
    					y[0] = y1[0];
    				}
    			}
    
    			if (sw == 0)
    				sw = -1;
    
    			return sw;
    		}
    
    		/******************************************/
    		/* 黄金分割法(与えられた方向での最小値)   */
    		/*      a,b : 初期区間 a < b              */
    		/*      eps : 許容誤差                    */
    		/*      val : 関数値                      */
    		/*      ind : 計算状況                    */
    		/*              >= 0 : 正常終了(収束回数) */
    		/*              = -1 : 収束せず           */
    		/*      max : 最大試行回数                */
    		/*      w : 位置                          */
    		/*      dw : 傾きの成分                   */
    		/*      fn : 関数値を計算する関数         */
    		/*      return : 結果(w+y*dwのy)        */
    		/******************************************/
    		function gold(a, b, eps, val, ind, max, w, dw, fn)
    		{
    			let f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
    			let count = 0;
    					// 初期設定
    			ind[0] = -1;
    			x1     = b - tau * (b - a);
    			x2     = a + tau * (b - a);
    			f1     = fn(x1, w, dw);
    			f2     = fn(x2, w, dw);
    					// 計算
    			while (count < max && ind[0] < 0) {
    				count += 1;
    				if (f2 > f1) {
    					if (Math.abs(b-a) < eps) {
    						ind[0] = 0;
    						x      = x1;
    						val[0] = f1;
    					}
    					else {
    						b  = x2;
    						x2 = x1;
    						x1 = a + (1.0 - tau) * (b - a);
    						f2 = f1;
    						f1 = fn(x1, w, dw);
    					}
    				}
    				else {
    					if (Math.abs(b-a) < eps) {
    						ind[0] = 0;
    						x      = x2;
    						val[0] = f2;
    						f1     = f2;
    					}
    					else {
    						a  = x1;
    						x1 = x2;
    						x2 = b - (1.0 - tau) * (b - a);
    						f1 = f2;
    						f2 = fn(x2, w, dw);
    					}
    				}
    			}
    					// 収束した場合の処理
    			if (ind[0] == 0) {
    				ind[0] = count;
    				fa     = fn(a, w, dw);
    				fb     = fn(b, w, dw);
    				if (fb < fa) {
    					a  = b;
    					fa = fb;
    				}
    				if (fa < f1) {
    					x      = a;
    					val[0] = fa;
    				}
    			}
    
    			return x;
    		}
    
    		/*********************/
    		/* Hesse行列の逆行列 */
    		/*********************/
    		function hesse(n, x, H, eps)
    		{
    			let sw1 = 0;
    
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++)
    					H[i1][i2] = eval(str3[i1*n+i2]);
    			}
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++) {
    					if (i1 == i2)
    						H[i1][n+i2] = 1.0;
    					else
    						H[i1][n+i2] = 0.0;
    				}
    			}
    
    			sw1 = Gauss(H, n, n, eps);
    
    			return sw1;
    		}
    
    		/******************************************/
    		/* 線形連立方程式を解く(逆行列を求める) */
    		/*      w : 方程式の左辺及び右辺          */
    		/*      n : 方程式の数                    */
    		/*      m : 方程式の右辺の列の数          */
    		/*      eps : 正則性を判定する規準        */
    		/*      return : =0 : 正常                */
    		/*               =1 : 逆行列が存在しない  */
    		/******************************************/
    		function Gauss(w, n, m, eps)
    		{
    			let y1, y2, ind = 0, nm, m1, m2;
    
    			nm = n + m;
    
    			for (let i1 = 0; i1 < n && ind == 0; i1++) {
    
    				y1 = .0;
    				m1 = i1 + 1;
    				m2 = 0;
    
    				for (let i2 = i1; i2 < n; i2++) {
    					y2 = Math.abs(w[i2][i1]);
    					if (y1 < y2) {
    						y1 = y2;
    						m2 = i2;
    					}
    				}
    
    				if (y1 < eps)
    					ind = 1;
    
    				else {
    
    					for (let i2 = i1; i2 < nm; i2++) {
    						y1        = w[i1][i2];
    						w[i1][i2] = w[m2][i2];
    						w[m2][i2] = y1;
    					}
    
    					y1 = 1.0 / w[i1][i1];
    
    					for (let i2 = m1; i2 < nm; i2++)
    						w[i1][i2] *= y1;
    
    					for (let i2 = 0; i2 < n; i2++) {
    						if (i2 != i1) {
    							for (let i3 = m1; i3 < nm; i3++)
    								w[i2][i3] -= w[i2][i1] * w[i1][i3];
    						}
    					}
    				}
    			}
    
    			return(ind);
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>最適化(Newton 法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,以下に示すような関数の最小値を求める場合に対する値が設定されています(他の問題を実行する場合は,それらを適切に修正してください).
    		<P STYLE="text-align:center"><IMG SRC="newton.gif"></P>
    		<DT>n 変数関数 f(<B>x</B>) に対する Newton 法においては,現在の点 <B>x</B> における傾き,
    		<P STYLE="text-align:center"><IMG SRC="newton1.gif"></P>
    		<DT>と Hesse 行列の計算結果から方向 <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 はαに対応する変数として表現されています.また,Hesse 行列は,1 行 1 列,1 行 2 列,・・・,n 行 n 列の順で入力して下さい. なお,目的関数,目的関数の微分,及び,Hesse 行列は,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="200"> 
    		刻み幅:<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.1"><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><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>
    		Hesse行列: H = <TEXTAREA ID="hesse" STYLE="font-size: 110%" COLS="70" ROWS="10">
    400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
    -400.0 * x[0]
    400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
    200.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 %lf %*s %lf", $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(2*$n);
    	}
    					// 実行
    	switch ($fun) {
    		case 1:
    			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx1", "dsnx1", "hesse1");
    			break;
    		case 2:
    			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx2", "dsnx2", "hesse2");
    			break;
    		case 3:
    			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx3", "dsnx3", "hesse3");
    			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);
    	}
    
    /********************************************************/
    /* Newton法                                             */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x : x(初期値と答え)                             */
    /*      dx : 関数の微分値                               */
    /*      H : Hesse行列の逆行列                           */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      hesse : Hesse行列の逆行列を計算する関数名       */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    
    function Newton($opt_1, $max, $n, $eps, $step, &$y, &$x, $dx, $H, $snx, $dsnx, $hesse)
    {
    	$wk    = array($n);
    	$count = 0;
    	$sw    = 0;
    
    	$y1 = $snx(0.0, $x, $dx);
    
    	while ($count < $max && $sw == 0) {
    					// 傾きの計算
    		$dsnx($x, $wk);
    					// Hesse行列の逆行列の計算
    		$sw1 = $hesse($x, $H, $eps);
    					// 収束していない
    		if ($sw1 == 0) {
    						// 方向の計算
    			$count++;
    			for ($i1 = 0; $i1 < $n; $i1++) {
    				$dx[$i1] = 0.0;
    				for ($i2 = 0; $i2 < $n; $i2++)
    					$dx[$i1] += $H[$i1][$n+$i2] * $wk[$i2];
    			}
    						// 1次元最適化を行わない
    			if ($opt_1 == 0) {
    							// 新しい点
    				for ($i1 = 0; $i1 < $n; $i1++)
    					$x[$i1] += $dx[$i1];
    							// 新しい関数値
    				$y2 = $snx(0.0, $x, $dx);
    								// 関数値の変化が大きい
    				if (abs($y2-$y1) > $eps)
    					$y1 = $y2;
    								// 収束(関数値の変化<eps)
    				else {
    					$sw = $count;
    					$y = $y2;
    				}
    			}
    						// 1次元最適化を行う
    			else {
    							// 区間を決める
    				$sw1 = 0;
    				$f1  = $y1;
    				$sp  = $step;
    				$f2  = $snx($sp, $x, $dx);
    				if ($f2 > $f1)
    					$sp = -$step;
    				for ($i1 = 0; $i1 < $max && $sw1 == 0; $i1++) {
    					$f2 = $snx($sp, $x, $dx);
    					if ($f2 > $f1)
    						$sw1 = 1;
    					else {
    						$sp *= 2.0;
    						$f1  = $f2;
    					}
    				}
    							// 区間が求まらない
    				if ($sw1 == 0)
    					$sw = -2;
    							// 区間が求まった
    				else {
    								// 黄金分割法
    					$k = gold(0.0, $sp, $eps, $y2, $sw1, $max, $x, $dx, $snx);
    								// 黄金分割法が失敗
    					if ($sw1 < 0)
    						$sw = -3;
    								// 黄金分割法が成功
    					else {
    									// 新しい点
    						for ($i1 = 0; $i1 < $n; $i1++)
    							$x[$i1] += $k * $dx[$i1];
    										// 関数値の変化が大きい
    						if (abs($y1-$y2) > $eps)
    							$y1 = $y2;
    										// 収束(関数値の変化<eps)
    						else {
    							$sw = $count;
    							$y = $y2;
    						}
    					}
    				}
    			}
    		}
    					// 収束(傾き<eps)
    		else {
    			$sw = $count;
    			$y = $y1;
    		}
    	}
    
    	if ($sw == 0)
    		$sw = -1;
    
    	return $sw;
    }
    
    /*******************************/
    /* Hesse行列の逆行列を計算する */
    /*******************************/
    					// 関数1
    function hesse1($x, &$H, $eps)
    {
    	$H[0][0] = 2.0;
    	$H[0][1] = 0.0;
    	$H[1][0] = 0.0;
    	$H[1][1] = 2.0;
    	$H[0][2] = 1.0;
    	$H[0][3] = 0.0;
    	$H[1][2] = 0.0;
    	$H[1][3] = 1.0;
    
    	$sw = Gauss($H, 2, 2, $eps);
    
    	return $sw;
    }
    					// 関数2
    function hesse2($x, &$H, $eps)
    {
    	$H[0][0] = 400.0 * (3.0 * $x[0] * $x[0] - $x[1]) + 2.0;
    	$H[0][1] = -400.0 * $x[0];
    	$H[1][0] = $H[0][1];
    	$H[1][1] = 200.0;
    	$H[0][2] = 1.0;
    	$H[0][3] = 0.0;
    	$H[1][2] = 0.0;
    	$H[1][3] = 1.0;
    
    	$sw = Gauss($H, 2, 2, $eps);
    
    	return $sw;
    }
    					// 関数3
    function hesse3($x, &$H, $eps)
    {
    	$x1      = 1.0 - $x[1];
    	$x2      = 1.0 - $x[1] * $x[1];
    	$x3      = 1.0 - $x[1] * $x[1] * $x[1];
    	$H[0][0] = 2.0 * $x1 * $x1 + 2.0 * $x2 * $x2 + 2.0 * $x3 * $x3;
    	$H[0][1] = 2.0 * (1.5 - $x[0] * $x1) - 2.0 * $x[0] * $x1 +
                  4.0 * $x[1] * (2.25 - $x[0] * $x2) - 4.0 * $x[0] * $x[1] * $x2 +
                  6.0 * $x[1] * $x[1] * (2.625 - $x[0] * $x3) - 6.0 * $x[0] * $x[1] * $x[1] * $x3;
    	$H[1][0] = $H[0][1];
    	$H[1][1] = 2.0 * $x[0] * $x[0] +
                  4.0 * $x[0] * (2.25 - $x[0] * $x2) + 8.0 * $x[0] * $x[0] * $x[1] * $x[1] +
                  12.0 * $x[0] * $x[1] * (2.625 - $x[0] * $x3) +
                  18.0 * $x[0] * $x[0] * $x[1] * $x[1] * $x[1] * $x[1];
    	$H[0][2] = 1.0;
    	$H[0][3] = 0.0;
    	$H[1][2] = 0.0;
    	$H[1][3] = 1.0;
    
    	$sw = Gauss($H, 2, 2, $eps);
    
    	return $sw;
    }
    
    /*******************************************************/
    /* 線形連立方程式を解く(逆行列を求める)              */
    /*      w : 方程式の左辺及び右辺                       */
    /*      n : 方程式の数                                 */
    /*      m : 方程式の右辺の列の数                       */
    /*      eps : 正則性を判定する規準                     */
    /*      return : =0 : 正常                             */
    /*               =1 : 逆行列が存在しない               */
    /*******************************************************/
    function Gauss(&$w, $n, $m, $eps)
    {
    	$ind = 0;
    	$nm  = $n + $m;
    
    	for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
    
    		$y1 = .0;
    		$m1 = $i1 + 1;
    		$m2 = 0;
    
    		for ($i2 = $i1; $i2 < $n; $i2++) {
    			$y2 = abs($w[$i2][$i1]);
    			if ($y1 < $y2) {
    				$y1 = $y2;
    				$m2 = $i2;
    			}
    		}
    
    		if ($y1 < $eps)
    			$ind = 1;
    
    		else {
    
    			for ($i2 = $i1; $i2 < $nm; $i2++) {
    				$y1          = $w[$i1][$i2];
    				$w[$i1][$i2] = $w[$m2][$i2];
    				$w[$m2][$i2] = $y1;
    			}
    
    			$y1 = 1.0 / $w[$i1][$i1];
    
    			for ($i2 = $m1; $i2 < $nm; $i2++)
    				$w[$i1][$i2] *= $y1;
    
    			for ($i2 = 0; $i2 < $n; $i2++) {
    				if ($i2 != $i1) {
    					for ($i3 = $m1; $i3 < $nm; $i3++)
    						$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
    				}
    			}
    		}
    	}
    
    	return($ind);
    }
    
    /*********************************************/
    /* 与えられた点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

      元来,関数 Gauss は,終了状態を表す変数 ind を return していました.しかし,この方法では,ブロック内で関数 Gauss を呼んだとき,ブロックに引数として渡した配列の値の変化を,ブロックを呼んだ側で知ることができません.そのため,関数 Gauss によって目的とする配列を return するように修正し,その配列の 1 行目の最後の要素に ind を記憶しています.ブロック付き関数呼び出しと C++ などにおける関数名を引数とする場合との違いが大きく出てしまいました.しかし,このままの仕様では,ブロック付き関数呼び出しにおいて,複数の種類のデータを返したいような場合は不便だと思います.
    #*****************************/
    # Newton法による最小値の計算 */
    #      coded by Y.Suganuma   */
    #*****************************/
    
    #******************************************************/
    # 線形連立方程式を解く(逆行列を求める)              */
    #      w : 方程式の左辺及び右辺                       */
    #      n : 方程式の数                                 */
    #      m : 方程式の右辺の列の数                       */
    #      eps : 正則性を判定する規準                     */
    #      return : =0 : 正常                             */
    #               =1 : 逆行列が存在しない               */
    #******************************************************/
    def Gauss(w, n, m, eps)
    
    	ind = 0
    	nm  = n + m
    
    	for i1 in 0 ... n
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    
    		for i2 in i1 ... n
    			y2 = w[i2][i1].abs()
    			if y1 < y2
    				y1 = y2
    				m2 = i2
    			end
    		end
    
    		if y1 < eps
    			ind = 1
    			break
    
    		else
    
    			for i2 in i1 ... nm
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    			end
    
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in m1 ... nm
    				w[i1][i2] *= y1
    			end
    
    			for i2 in 0 ... n
    				if i2 != i1
    					for i3 in m1 ... nm
    						w[i2][i3] -= w[i2][i1] * w[i1][i3]
    					end
    				end
    			end
    		end
    	end
    
    	w[0][nm] = ind
    	return w
    end
    
    #******************************************************/
    # 与えられた点xから,dx方向へk*dxだけ進んだ点における */
    # 関数値,微係数,及び,Hesse行列の逆行列を計算する   */
    #******************************************************/
    					# 関数1
    snx1 = Proc.new { |sw, k, x, dx, h, eps|
    	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
    	elsif sw == 1
    		dx[0] = -2.0 * (x[0] - 1.0)
    		dx[1] = -2.0 * (x[1] - 2.0)
    	else
    		h = [[2.0, 0.0, 1.0, 0.0, 0.0], [0.0, 2.0, 0.0, 1.0, 0.0]]
    		Gauss(h, 2, 2, eps)
    	end
    }
    					# 関数2
    snx2 = Proc.new { |sw, k, x, dx, h, eps|
    	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
    	elsif sw == 1
    		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])
    	else
    		h[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
    		h[0][1] = -400.0 * x[0]
    		h[1][0] = h[0][1]
    		h[1][1] = 200.0
    		h[0][2] = 1.0
    		h[0][3] = 0.0
    		h[1][2] = 0.0
    		h[1][3] = 1.0
    		Gauss(h, 2, 2, eps)
    	end
    }
    					# 関数3
    snx3 = Proc.new { |sw, k, x, dx, h, eps|
    	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
    	elsif sw == 1
    		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]))
    	else
    		x1      = 1.0 - x[1]
    		x2      = 1.0 - x[1] * x[1]
    		x3      = 1.0 - x[1] * x[1] * x[1]
    		h[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
    		h[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                      4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                      6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
    		h[1][0] = h[0][1]
    		h[1][1] = 2.0 * x[0] * x[0] +
                      4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                      12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                      18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
    		h[0][2] = 1.0
    		h[0][3] = 0.0
    		h[1][2] = 0.0
    		h[1][3] = 1.0
    		Gauss(h, 2, 2, eps)
    	end
    }
    
    #***************************************************************/
    # 黄金分割法(与えられた方向での最小値)                         */
    #      a,b : 初期区間 a < b                                    */
    #      eps : 許容誤差                                          */
    #      val : 関数値                                            */
    #      ind : 計算状況                                          */
    #              >= 0 : 正常終了(収束回数)                       */
    #              = -1 : 収束せず                                 */
    #      max : 最大試行回数                                      */
    #      w : 位置                                                */
    #      dw : 傾きの成分                                         */
    #      h : Hesse行列の逆行列                                   */
    #      snx : 関数値を計算する関数の名前                        */
    #            (微分は,その符号を変える)                       */
    #      return : 結果(w+y*dwのy)                              */
    #***************************************************************/
    def gold(a, b, eps, val, ind, max, w, dw, h, &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, h, eps)
    	f2     = snx.call(0, x2, w, dw, h, eps)
    					# 計算
    	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, h, eps)
    			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, h, eps)
    			end
    		end
    	end
    					# 収束した場合の処理
    	if ind[0] == 0
    		ind[0] = count
    		fa     = snx.call(0, a, w, dw, h, eps)
    		fb     = snx.call(0, b, w, dw, h, eps)
    		if fb < fa
    			a  = b
    			fa = fb
    		end
    		if fa < f1
    			x      = a
    			val[0] = fa
    		end
    	end
    
    	return x
    end
    
    #*******************************************************/
    # Newton法                                             */
    #      opt_1 : =0 : 1次元最適化を行わない             */
    #              =1 : 1次元最適化を行う                 */
    #      max : 最大繰り返し回数                          */
    #      n : 次元                                        */
    #      eps : 収束判定条件                              */
    #      step : きざみ幅                                 */
    #      y : 最小値                                      */
    #      x : x(初期値と答え)                             */
    #      dx : 関数の微分値                               */
    #      h : Hesse行列の逆行列                           */
    #      snx : 関数値,その微分,及び,Hesse行列の逆行列 */
    #      return : >=0 : 正常終了(収束回数)               */
    #               =-1 : 収束せず                         */
    #               =-2 : 1次元最適化の区間が求まらない   */
    #               =-3 : 黄金分割法が失敗                 */
    #*******************************************************/
    def Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx)
    
    	sw1 = Array.new(n)
    	y2  = Array.new(n)
    	wk  = Array.new(n)
    	for i1 in 0 ... n
    		wk[i1] = 0.0
    	end
    	count = 0
    	sw    = 0
    
    	y1 = snx.call(0, 0.0, x, dx, h, eps)
    
    	while count < max && sw == 0
    					# 傾きの計算
    		snx.call(1, 0.0, x, wk, h, eps)
    					# Hesse行列の逆行列の計算
    		h      = snx.call(2, 0.0, x, wk, h, eps)
    		sw1[0] = h[0][2*n]   # もっと良い方法はあるか?
    					# 収束していない
    		if sw1[0] == 0
    						# 方向の計算
    			count += 1
    			for i1 in 0 ... n
    				dx[i1] = 0.0
    				for i2 in 0 ... n
    					dx[i1] += h[i1][n+i2] * wk[i2]
    				end
    			end
    						# 1次元最適化を行わない
    			if opt_1 == 0
    							# 新しい点
    				for i1 in 0 ... n
    					x[i1] += dx[i1]
    				end
    							# 新しい関数値
    				y2[0] = snx.call(0, 0.0, x, dx, h, eps)
    								# 関数値の変化が大きい
    				if (y2[0]-y1).abs() > eps
    					y1 = y2[0]
    								# 収束(関数値の変化<eps)
    				else
    					sw   = count
    					y[0] = y2[0]
    				end
    						# 1次元最適化を行う
    			else
    							# 区間を決める
    				sw1[0] = 0
    				f1     = y1
    				sp     = step
    				f2     = snx.call(0, sp, x, dx, h, eps)
    				if f2 > f1
    					sp = -step
    				end
    				for i1 in 0 ... max
    					f2 = snx.call(0, sp, x, dx, h, eps)
    					if f2 > f1
    						sw1[0] = 1
    						break
    					else
    						sp *= 2.0
    						f1  = f2
    					end
    				end
    							# 区間が求まらない
    				if sw1[0] == 0
    					sw = -2
    							# 区間が求まった
    				else
    								# 黄金分割法
    					k = gold(0.0, sp, eps, y2, sw1, max, x, dx, h, &snx)
    								# 黄金分割法が失敗
    					if sw1[0] < 0
    						sw = -3
    								# 黄金分割法が成功
    					else
    									# 新しい点
    						for i1 in 0 ... n
    							x[i1] += k * dx[i1]
    						end
    										# 関数値の変化が大きい
    						if (y1-y2[0]).abs() > eps
    							y1 = y2[0]
    										# 収束(関数値の変化<eps)
    						else
    							sw   = count
    							y[0] = y2[0]
    						end
    					end
    				end
    			end
    					# 収束(傾き<eps)
    		else
    			sw = count
    			y[0] = y1
    		end
    	end
    
    	if sw == 0
    		sw = -1
    	end
    
    	return sw
    end
    				# データの入力
    str   = gets()
    a     = str.split(" ")
    fun   = Integer(a[1])
    n     = Integer(a[3])
    max   = Integer(a[5])
    opt_1 = Integer(a[7])
    str   = gets()
    a     = str.split(" ")
    eps   = Float(a[1])
    step  = Float(a[3])
    x     = Array.new(n)
    dx    = Array.new(n)
    y     = Array.new(1)
    sw    = 0
    str   = gets()
    a     = str.split(" ")
    h     = Array.new(n)
    for i1 in 0 ... n
    	x[i1]  = Float(a[i1+1])
    	dx[i1] = 0.0
    	h[i1]  = Array.new(2*n+1)
    end
    				# 実行
    case fun
    	when 1
    		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx1)
    	when 2
    		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx2)
    	when 3
    		sw = Newton(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
    
    ############################################
    # 線形連立方程式を解く(逆行列を求める)
    #      w : 方程式の左辺及び右辺
    #      n : 方程式の数
    #      m : 方程式の右辺の列の数
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    #      coded by Y.Suganuma
    ############################################
    
    def gauss(w, n, m, eps) :
    
    	nm  = n + m
    	ind = 0
    
    	for i1 in range(0, n) :
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    					# ピボット要素の選択
    		for i2 in range(i1, n) :
    			y2 = abs(w[i2][i1])
    			if y1 < y2 :
    				y1 = y2
    				m2 = i2
    					# 逆行列が存在しない
    		if y1 < eps :
    			ind = 1
    			break
    					# 逆行列が存在する
    		else :   # 行の入れ替え
    			for i2 in range(i1, nm) :
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    						# 掃き出し操作
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in range(m1, nm) :
    				w[i1][i2] *= y1
    
    			for i2 in range(0, n) :
    				if i2 != i1 :
    					for i3 in range(m1, nm) :
    						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
    
    	return ind
    
    ############################################
    # Newton法
    #      opt_1 : =0 : 1次元最適化を行わない
    #              =1 : 1次元最適化を行う
    #      max : 最大繰り返し回数
    #      n : 次元
    #      eps : 収束判定条件
    #      step : きざみ幅
    #      y : 最小値
    #      x : x(初期値と答え)
    #      dx : 関数の微分値
    #      H : Hesse行列の逆行列
    #      snx : 関数値を計算する関数名
    #      dsnx : 関数の微分を計算する関数名(符号を変える)
    #      hesse : Hesse行列の逆行列を計算する関数名
    #      return : >=0 : 正常終了(収束回数)
    #               =-1 : 収束せず
    #               =-2 : 1次元最適化の区間が求まらない
    #               =-3 : 黄金分割法が失敗
    #      coded by Y.Suganuma
    ############################################
    
    def Newton(opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx, hesse) :
    
    	wk    = np.empty(n, np.float)
    	count = 0
    	sw    = 0
    	y2    = np.empty(1, np.float)
    	sw1   = np.empty(1, np.int)
    
    	y1    = snx(0.0, x, dx)
    
    	while count < max and sw == 0 : 
    			# 傾きの計算
    		dsnx(x, wk)
    			# Hesse行列の逆行列の計算
    		sw1[0] = hesse(x, H, eps)
    			# 収束していない
    		if sw1[0] == 0 :
    				# 方向の計算
    			count += 1
    			for i1 in range(0, n) :
    				dx[i1] = 0.0
    				for i2 in range(0, n) :
    					dx[i1] += H[i1][n+i2] * wk[i2]
    				# 1次元最適化を行わない
    			if opt_1 == 0 :
    					# 新しい点
    				for i1 in range(0, n) :
    					x[i1] += dx[i1]
    					# 新しい関数値
    				y2[0] = snx(0.0, x, dx)
    						# 関数値の変化が大きい
    				if abs(y2[0]-y1) > eps :
    					y1 = y2[0]
    						# 収束(関数値の変化<eps)
    				else :
    					sw   = count
    					y[0] = y2[0]
    				# 1次元最適化を行う
    			else :
    					# 区間を決める
    				sw1[0] = 0
    				f1     = y1
    				sp     = step
    				f2     = snx(sp, x, dx)
    				if f2 > f1 :
    					sp = -step
    				for i1 in range(0, max) :
    					f2 = snx(sp, x, dx)
    					if f2 > f1 :
    						sw1[0] = 1
    						break
    					else :
    						sp *= 2.0
    						f1  = f2
    					# 区間が求まらない
    				if sw1[0] == 0 :
    					sw = -2
    					# 区間が求まった
    				else :
    						# 黄金分割法
    					k = gold(0.0, sp, eps, y2, sw1, max, x, dx, snx)
    						# 黄金分割法が失敗
    					if sw1[0] < 0 :
    						sw = -3
    						# 黄金分割法が成功
    					else :
    							# 新しい点
    						for i1 in range(0, n) :
    							x[i1] += k * dx[i1]
    								# 関数値の変化が大きい
    						if abs(y1-y2[0]) > eps :
    							y1 = y2[0]
    								# 収束(関数値の変化<eps)
    						else :
    							sw   = count
    							y[0] = y2[0]
    			# 収束(傾き<eps)
    		else :
    			sw   = count
    			y[0] = y1
    
    	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]))
    
    ############################################
    # Hesse行列の逆行列を計算する */
    ############################################
    
    			# 関数1
    def hesse1(x, H, eps) :
    
    	H[0][0] = 2.0
    	H[0][1] = 0.0
    	H[1][0] = 0.0
    	H[1][1] = 2.0
    	H[0][2] = 1.0
    	H[0][3] = 0.0
    	H[1][2] = 0.0
    	H[1][3] = 1.0
    	sw = gauss(H, 2, 2, eps)
    	return sw
    
    			# 関数2
    def hesse2(x, H, eps) :
    
    	H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
    	H[0][1] = -400.0 * x[0]
    	H[1][0] = H[0][1]
    	H[1][1] = 200.0
    	H[0][2] = 1.0
    	H[0][3] = 0.0
    	H[1][2] = 0.0
    	H[1][3] = 1.0
    	sw = gauss(H, 2, 2, eps)
    	return sw
    
    			# 関数3
    def hesse3(x, H, eps) :
    
    	x1      = 1.0 - x[1]
    	x2      = 1.0 - x[1] * x[1]
    	x3      = 1.0 - x[1] * x[1] * x[1]
    	H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
    	H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 + 4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 + 6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
    	H[1][0] = H[0][1]
    	H[1][1] = 2.0 * x[0] * x[0] + 4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] + 12.0 * x[0] * x[1] * (2.625 - x[0] * x3) + 18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
    	H[0][2] = 1.0
    	H[0][3] = 0.0
    	H[1][2] = 0.0
    	H[1][3] = 1.0
    	sw = gauss(H, 2, 2, eps)
    	return sw
    
    			# データの入力
    sw    = 0
    line  = sys.stdin.readline()
    s     = line.split()
    fun   = int(s[1])
    n     = int(s[3])
    max   = int(s[5])
    opt_1 = int(s[7])
    
    line  = sys.stdin.readline()
    s     = line.split()
    eps   = float(s[1])
    step  = float(s[3])
    
    x     = np.empty(n, np.float)
    dx    = np.empty(n, np.float)
    H     = np.empty((n, 2*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 = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1, hesse1)
    elif fun == 2 :
    	sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2, hesse2)
    else :
    	sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3, hesse3)
    			# 結果の出力
    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);
    		double eps  = double.Parse(str[1]);
    		double step = double.Parse(str[3]);
    							// 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 [2*n];
    		}
    					// 実行
    		int sw;
    		if (fun == 1)
    			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx1, dsnx1, hesse1);
    		else if (fun == 2)
    			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx2, dsnx2, hesse2);
    		else
    			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx3, dsnx3, hesse3);
    					// 結果の出力
    		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;
    	}
    					// Hesse行列の逆行列1
    	int hesse1(double[] x, double[][] H, double eps)
    	{
    		H[0][0] = 2.0;
    		H[0][1] = 0.0;
    		H[1][0] = 0.0;
    		H[1][1] = 2.0;
    		H[0][2] = 1.0;
    		H[0][3] = 0.0;
    		H[1][2] = 0.0;
    		H[1][3] = 1.0;
    		int sw1 = gauss(H, 2, 2, eps);
    		return sw1;
    	}
    					// Hesse行列の逆行列2
    	int hesse2(double[] x, double[][] H, double eps)
    	{
    		H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
    		H[0][1] = -400.0 * x[0];
    		H[1][0] = H[0][1];
    		H[1][1] = 200.0;
    		H[0][2] = 1.0;
    		H[0][3] = 0.0;
    		H[1][2] = 0.0;
    		H[1][3] = 1.0;
    		int sw1 = gauss(H, 2, 2, eps);
    		return sw1;
    	}
    					// Hesse行列の逆行列3
    	int hesse3(double[] x, double[][] H, double eps)
    	{
    		double x1 = 1.0 - x[1];
    		double x2 = 1.0 - x[1] * x[1];
    		double x3 = 1.0 - x[1] * x[1] * x[1];
    		H[0][0]   = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
    		H[0][1]   = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
    	                4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
    	                6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
    		H[1][0] = H[0][1];
    		H[1][1] = 2.0 * x[0] * x[0] +
    	              4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
    	              12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
    	              18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
    		H[0][2] = 1.0;
    		H[0][3] = 0.0;
    		H[1][2] = 0.0;
    		H[1][3] = 1.0;
    		int sw1 = gauss(H, 2, 2, eps);
    		return sw1;
    	}
    
    	/********************************************************/
    	/* Newton法                                             */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x : x(初期値と答え)                             */
    	/*      dx : 関数の微分値                               */
    	/*      H : Hesse行列の逆行列                           */
    	/*      fn : 関数値を計算する関数                       */
    	/*      dfn : 関数の微分値を計算する関数                */
    	/*      hesse : Hesse行列の逆行列を計算する関数         */
    	/*      return : >=0 : 正常終了(収束回数)               */
    	/*               =-1 : 収束せず                         */
    	/*               =-2 : 1次元最適化の区間が求まらない   */
    	/*               =-3 : 黄金分割法が失敗                 */
    	/********************************************************/
    	int Newton(int opt_1, int max, int n, double eps, double step, ref double y,
    	           double[] x, double[] dx, double[][] H,
    	           Func<double, double[], double[], double> fn,
    	           Func<double[], double[], int> dfn,
    	           Func<double[], double[][], double, int> hesse)
    	{
    		double[] wk = new double [n];
    		double y1   = fn(0.0, x, dx);
    		double y2   = 0.0;
    		int count   = 0;
    		int sw      = 0;
    
    		while (count < max && sw == 0) {
    					// 傾きの計算
    			dfn(x, wk);
    					// Hesse行列の逆行列の計算
    			int sw1 = hesse(x, H, eps);
    					// 収束していない
    			if (sw1 == 0) {
    						// 方向の計算
    				count++;
    				for (int i1 = 0; i1 < n; i1++) {
    					dx[i1] = 0.0;
    					for (int i2 = 0; i2 < n; i2++)
    						dx[i1] += H[i1][n+i2] * wk[i2];
    				}
    						// 1次元最適化を行わない
    				if (opt_1 == 0) {
    							// 新しい点
    					for (int i1 = 0; i1 < n; i1++)
    						x[i1] += dx[i1];
    							// 新しい関数値
    					y2 = fn(0.0, x, dx);
    								// 関数値の変化が大きい
    					if (Math.Abs(y2-y1) > eps)
    						y1 = y2;
    								// 収束(関数値の変化<eps)
    					else {
    						sw = count;
    						y  = y2;
    					}
    				}
    						// 1次元最適化を行う
    				else {
    							// 区間を決める
    					sw1       = 0;
    					double f1 = y1;
    					double sp = step;
    					double f2 = fn(sp, x, dx);
    					if (f2 > f1)
    						sp = -step;
    					for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
    						f2 = fn(sp, x, dx);
    						if (f2 > f1)
    							sw1 = 1;
    						else {
    							sp *= 2.0;
    							f1  = f2;
    						}
    					}
    							// 区間が求まらない
    					if (sw1 == 0)
    						sw = -2;
    							// 区間が求まった
    					else {
    								// 黄金分割法
    						double k = gold(0.0, sp, eps, ref y2, ref sw1, max, x, dx, fn);
    								// 黄金分割法が失敗
    						if (sw1 < 0)
    							sw = -3;
    								// 黄金分割法が成功
    						else {
    									// 新しい点
    							for (int i1 = 0; i1 < n; i1++)
    								x[i1] += k * dx[i1];
    										// 関数値の変化が大きい
    							if (Math.Abs(y1-y2) > eps)
    								y1 = y2;
    										// 収束(関数値の変化<eps)
    							else {
    								sw = count;
    								y  = y2;
    							}
    						}
    					}
    				}
    			}
    					// 収束(傾き<eps)
    			else {
    				sw = count;
    				y  = y1;
    			}
    		}
    
    		if (sw == 0)
    			sw = -1;
    
    		return sw;
    	}
    
    	/******************************************/
    	/* 黄金分割法(与えられた方向での最小値)   */
    	/*      a,b : 初期区間 a < b              */
    	/*      eps : 許容誤差                    */
    	/*      val : 関数値                      */
    	/*      ind : 計算状況                    */
    	/*              >= 0 : 正常終了(収束回数) */
    	/*              = -1 : 収束せず           */
    	/*      max : 最大試行回数                */
    	/*      w : 位置                          */
    	/*      dw : 傾きの成分                   */
    	/*      fn : 関数値を計算する関数         */
    	/*      return : 結果(w+y*dwのy)        */
    	/******************************************/
    	double gold(double a, double b, double eps, ref double val, ref int ind,
    	            int max, double[] w, double[] dw,
    	            Func<double, double[], double[], double> fn)
    	{
    					// 初期設定
    		double tau = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0;
    		double x1  = b - tau * (b - a);
    		double x2  = a + tau * (b - a);
    		double f1  = fn(x1, w, dw);
    		double f2  = fn(x2, w, dw);
    		int count  = 0;
    		ind        = -1;
    					// 計算
    		while (count < max && ind < 0) {
    			count += 1;
    			if (f2 > f1) {
    				if (Math.Abs(b-a) < eps) {
    					ind = 0;
    					x   = x1;
    					val = f1;
    				}
    				else {
    					b  = x2;
    					x2 = x1;
    					x1 = a + (1.0 - tau) * (b - a);
    					f2 = f1;
    					f1 = fn(x1, w, dw);
    				}
    			}
    			else {
    				if (Math.Abs(b-a) < eps) {
    					ind = 0;
    					x   = x2;
    					val = f2;
    					f1  = f2;
    				}
    				else {
    					a  = x1;
    					x1 = x2;
    					x2 = b - (1.0 - tau) * (b - a);
    					f1 = f2;
    					f2 = fn(x2, w, dw);
    				}
    			}
    		}
    					// 収束した場合の処理
    		if (ind == 0) {
    			ind       = count;
    			double fa = fn(a, w, dw);
    			double fb = fn(b, w, dw);
    			if (fb < fa) {
    				a  = b;
    				fa = fb;
    			}
    			if (fa < f1) {
    				x   = a;
    				val = fa;
    			}
    		}
    
    		return x;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 逆行列の存在を判定する規準               */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	int gauss(double[][] w, int n, int m, double eps)
    	{
    		int ind = 0;
    		int nm  = n + m;
    
    		for (int i1 = 0; i1 < n && ind == 0; i1++) {
    
    			double y1 = .0;
    			int m1    = i1 + 1;
    			int m2    = 0;
    						// ピボット要素の選択
    			for (int i2 = i1; i2 < n; i2++) {
    				double y2 = Math.Abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (int i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (int i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (int i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (int i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    			

  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 eps As Double   = Double.Parse(str(1))
    		Dim step1 As Double = Double.Parse(str(3))
    							' 3 行目
    		Dim H(n,2*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
    					' Hesse行列の逆行列1(ラムダ式)
    		Dim hesse1 = Function(x0, H0, eps0)
    			H0(0,0) = 2.0
    			H0(0,1) = 0.0
    			H0(1,0) = 0.0
    			H0(1,1) = 2.0
    			H0(0,2) = 1.0
    			H0(0,3) = 0.0
    			H0(1,2) = 0.0
    			H0(1,3) = 1.0
    			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
    			Return sw1
    		End Function
    					' Hesse行列の逆行列2(ラムダ式)
    		Dim hesse2 = Function(x0, H0, eps0)
    			H0(0,0) = 400.0 * (3.0 * x0(0) * x0(0) - x0(1)) + 2.0
    			H0(0,1) = -400.0 * x0(0)
    			H0(1,0) = H0(0,1)
    			H0(1,1) = 200.0
    			H0(0,2) = 1.0
    			H0(0,3) = 0.0
    			H0(1,2) = 0.0
    			H0(1,3) = 1.0
    			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
    			Return sw1
    		End Function
    					' Hesse行列の逆行列3(ラムダ式)
    		Dim hesse3 = Function(x0, H0, eps0)
    			Dim x1 As Double = 1.0 - x0(1)
    			Dim x2 As Double = 1.0 - x0(1) * x0(1)
    			Dim x3 As Double = 1.0 - x0(1) * x0(1) * x0(1)
    			H0(0,0) = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
    			H0(0,1) = 2.0 * (1.5 - x0(0) * x1) - 2.0 * x0(0) * x1 +
    			            4.0 * x0(1) * (2.25 - x0(0) * x2) - 4.0 * x0(0) * x0(1) * x2 +
    			            6.0 * x0(1) * x0(1) * (2.625 - x0(0) * x3) - 6.0 * x0(0) * x0(1) * x0(1) * x3
    			H0(1,0) = H0(0,1)
    			H0(1,1) = 2.0 * x0(0) * x0(0) +
    			          4.0 * x0(0) * (2.25 - x0(0) * x2) + 8.0 * x0(0) * x0(0) * x0(1) * x0(1) +
    			          12.0 * x0(0) * x0(1) * (2.625 - x0(0) * x3) +
    			          18.0 * x0(0) * x0(0) * x0(1) * x0(1) * x0(1) * x0(1)
    			H0(0,2) = 1.0
    			H0(0,3) = 0.0
    			H0(1,2) = 0.0
    			H0(1,3) = 1.0
    			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
    			Return sw1
    		End Function
    					' 実行
    		Dim sw As Integer
    		If fun = 1
    			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx1, dsnx1, hesse1)
    		ElseIf fun = 2
    			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx2, dsnx2, hesse2)
    		Else
    			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx3, dsnx3, hesse3)
    		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
    
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' Newton法                                             '
    	'      opt_1 : =0 : 1次元最適化を行わない             '
    	'              =1 : 1次元最適化を行う                 '
    	'      max : 最大繰り返し回数                          '
    	'      n : 次元                                        '
    	'      eps : 収束判定条件                              '
    	'      step1 : きざみ幅                                '
    	'      y : 最小値                                      '
    	'      x : x(初期値と答え)                             '
    	'      dx : 関数の微分値                               '
    	'      H : Hesse行列の逆行列                           '
    	'      fn : 関数値を計算する関数                       '
    	'      dfn : 関数の微分値を計算する関数                '
    	'      hesse : Hesse行列の逆行列を計算する関数         '
    	'      return : >=0 : 正常終了(収束回数)               '
    	'               =-1 : 収束せず                         '
    	'               =-2 : 1次元最適化の区間が求まらない   '
    	'               =-3 : 黄金分割法が失敗                 '
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function Newton(opt_1 As Integer, max As Integer, n As Integer, eps As Double,
    	                step1 As Double, ByRef y As Double, x() As Double, dx() As Double,
    	                H(,) As Double, fn As Func(Of Double, Double(), Double(), Double),
    	                dfn As Func(Of Double(), Double(), Integer),
    	                hesse As Func(Of Double(), Double(,), Double, Integer))
    
    		Dim wk(n) As Double
    		Dim y1 As Double     = fn(0.0, x, dx)
    		Dim y2 As Double     = 0.0
    		Dim count As Integer = 0
    		Dim sw As Integer    = 0
    
    		Do while count < max and sw = 0
    					' 傾きの計算
    			dfn(x, wk)
    					' Hesse行列の逆行列の計算
    			Dim sw1 As Integer = hesse(x, H, eps)
    					' 収束していない
    			If sw1 = 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,n+i2) * wk(i2)
    					Next
    				Next
    						' 1次元最適化を行わない
    				If opt_1 = 0
    							' 新しい点
    					For i1 As Integer = 0 To n-1
    						x(i1) += dx(i1)
    					Next
    							' 新しい関数値
    					y2 = fn(0.0, x, dx)
    								' 関数値の変化が大きい
    					If Math.Abs(y2-y1) > eps
    						y1 = y2
    								' 収束(関数値の変化<eps)
    					Else
    						sw = count
    						y  = y2
    					End If
    						' 1次元最適化を行う
    				Else
    							' 区間を決める
    					sw1 = 0
    					Dim f1 As Double = y1
    					Dim sp As Double = step1
    					Dim f2 As Double = fn(sp, x, dx)
    					If f2 > f1
    						sp = -step1
    					End If
    					Dim ii As Integer = 0
    					Do While ii < max and sw1 = 0
    						f2 = fn(sp, x, dx)
    						If f2 > f1
    							sw1 = 1
    						Else
    							sp *= 2.0
    							f1  = f2
    						End If
    						ii += 1
    					Loop
    							' 区間が求まらない
    					If sw1 = 0
    						sw = -2
    							' 区間が求まった
    					Else
    								' 黄金分割法
    						Dim k As Double = gold(0.0, sp, eps, y2, sw1, max, x, dx, fn)
    								' 黄金分割法が失敗
    						If sw1 < 0
    							sw = -3
    								' 黄金分割法が成功
    						Else
    									' 新しい点
    							For i1 As Integer = 0 To n-1
    								x(i1) += k * dx(i1)
    							Next
    										' 関数値の変化が大きい
    							If Math.Abs(y1-y2) > eps
    								y1 = y2
    										' 収束(関数値の変化<eps)
    							Else
    								sw = count
    								y  = y2
    							End If
    						End If
    					End If
    				End If
    					' 収束(傾き<eps)
    			Else
    				sw = count
    				y  = y1
    			End If
    		Loop
    
    		If sw = 0
    			sw = -1
    		End If
    
    		Return sw
    
    	End Function
    
    	''''''''''''''''''''''''''''''''''''''''''
    	' 黄金分割法(与えられた方向での最小値)   '
    	'      a,b : 初期区間 a < b              '
    	'      eps : 許容誤差                    '
    	'      val : 関数値                      '
    	'      ind : 計算状況                    '
    	'              >= 0 : 正常終了(収束回数) '
    	'              = -1 : 収束せず           '
    	'      max : 最大試行回数                '
    	'      w : 位置                          '
    	'      dw : 傾きの成分                   '
    	'      fn : 関数値を計算する関数         '
    	'      return : 結果(w+y*dwのy)        '
    	''''''''''''''''''''''''''''''''''''''''''
    	Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
    	              ByRef ind As Integer, max As Integer, w() As Double, dw() As Double,
    	              fn As Func(Of Double, Double(), Double(), Double))
    
    					' 初期設定
    		Dim tau As Double    = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0
    		Dim x1 As Double     = b - tau * (b - a)
    		Dim x2 As Double     = a + tau * (b - a)
    		Dim f1 As Double     = fn(x1, w, dw)
    		Dim f2 As Double     = fn(x2, w, dw)
    		Dim count As Integer = 0
    		ind = -1
    					' 計算
    		Do While count < max and ind < 0
    			count += 1
    			If f2 > f1
    				If Math.Abs(b-a) < eps
    					ind = 0
    					x   = x1
    					val = f1
    				Else
    					b  = x2
    					x2 = x1
    					x1 = a + (1.0 - tau) * (b - a)
    					f2 = f1
    					f1 = fn(x1, w, dw)
    				End If
    			Else
    				If Math.Abs(b-a) < eps
    					ind = 0
    					x   = x2
    					val = f2
    					f1  = f2
    				Else
    					a  = x1
    					x1 = x2
    					x2 = b - (1.0 - tau) * (b - a)
    					f1 = f2
    					f2 = fn(x2, w, dw)
    				End If
    			End If
    		Loop
    					' 収束した場合の処理
    		If ind = 0
    			ind = count
    			Dim fa As Double = fn(a, w, dw)
    			Dim fb As Double = fn(b, w, dw)
    			If fb < fa
    				a  = b
    				fa = fb
    			End If
    			If fa < f1
    				x   = a
    				val = fa
    			End If
    		End If
    
    		Return x
    
    	End Function
    
    	''''''''''''''''''''''''''''''''''''''''''
    	' 線形連立方程式を解く(逆行列を求める) '
    	'      w : 方程式の左辺及び右辺          '
    	'      n : 方程式の数                    '
    	'      m : 方程式の右辺の列の数          '
    	'      eps : 逆行列の存在を判定する規準  '
    	'      return : =0 : 正常                '
    	'               =1 : 逆行列が存在しない  '
    	''''''''''''''''''''''''''''''''''''''''''
    	Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer
    
    		Dim ind As Integer = 0
    		Dim nm As Integer  = n + m
    
    		Dim i1 As Integer = 0
    		Do While i1 < n and ind = 0
    
    			Dim y1 As Double  = 0.0
    			Dim m1 As Integer = i1 + 1
    			Dim m2 As Integer = 0
    						' ピボット要素の選択
    			For i2 As Integer = i1 To n-1
    				Dim y2 As Double = Math.Abs(w(i2,i1))
    				If y1 < y2
    					y1 = y2
    					m2 = i2
    				End If
    			Next
    						' 逆行列が存在しない
    			If y1 < eps
    				ind = 1
    						' 逆行列が存在する
    			Else
    							' 行の入れ替え
    				For i2 As Integer = i1 To nm-1
    					y1       = w(i1,i2)
    					w(i1,i2) = w(m2,i2)
    					w(m2,i2) = y1
    				Next
    							' 掃き出し操作
    				y1 = 1.0 / w(i1,i1)
    
    				For i2 As Integer = m1 To nm-1
    					w(i1,i2) *= y1
    				Next
    
    				For i2 As Integer = 0 To n-1
    					If i2 <> i1
    						For i3 As Integer = m1 To nm-1
    							w(i2,i3) -= w(i2,i1) * w(i1,i3)
    						Next
    					End If
    				Next
    			End If
    			i1 += 1
    		Loop
    
    		Return ind
    
    	End Function
    
    End Module
    			

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