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

最適化(最急降下法)

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

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

  1. C++

    /********************************/
    /* 最急降下法による最小値の計算 */
    /*      coded by Y.Suganuma     */
    /********************************/
    #include <stdio.h>
    
    void dsnx1(double *, double *);
    double snx1(double, double *, double *);
    void dsnx2(double *, double *);
    double snx2(double, double *, double *);
    void dsnx3(double *, double *);
    double snx3(double, double *, double *);
    int steep(int, int, int, double, double, double *, double *, double *,
              double (*)(double, double *, double *),
              void (*)(double *, double *));
    
    int main()
    {
    	double eps, 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];
    	scanf("%*s");
    	for (i1 = 0; i1 < n; i1++)
    		scanf("%lf", &x[i1]);
    					// 実行
    	switch (fun) {
    		case 1:
    			sw = steep(opt_1, max, n, eps, step, &y, x, dx, snx1, dsnx1);
    			break;
    		case 2:
    			sw = steep(opt_1, max, n, eps, step, &y, x, dx, snx2, dsnx2);
    			break;
    		case 3:
    			sw = steep(opt_1, max, n, eps, step, &y, x, dx, snx3, dsnx3);
    			break;
    	}
    					// 結果の出力
    	if (sw < 0) {
    		printf("   収束しませんでした!");
    		switch (sw) {
    			case -1:
    				printf("(収束回数)\n");
    				break;
    			case -2:
    				printf("(1次元最適化の区間)\n");
    				break;
    			case -3:
    				printf("(黄金分割法)\n");
    				break;
    		}
    	}
    	else {
    		printf("   結果=");
    		for (i1 = 0; i1 < n; i1++)
    			printf("%f ", x[i1]);
    		printf(" 最小値=%f  回数=%d\n", y, sw);
    	}
    
    	return 0;
    }
    
    /********************************************************/
    /* 最急降下法(steepest descent method)                  */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x : x(初期値と答え)                             */
    /*      dx : 関数の微分値                               */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    #include <math.h>
    
    double gold(double, double, double, double *, int *, int, double *, double *,
                double (*)(double, double *, double *));
    
    int steep(int opt_1, int max, int n, double eps, double step, double *y,
              double *x, double *dx, double (*snx)(double, double *, double *), 
              void (*dsnx)(double *, double *))
    {
    	double f1, f2, k, sp, y1, y2;
    	int count = 0, i1, sw = 0, sw1;
    
    	y1 = snx(0.0, x, dx);
    
    	while (count < max && sw == 0) {
    					// 傾きの計算
    		dsnx(x, dx);
    		count++;
    					// 1次元最適化を行わない
    		if (opt_1 == 0) {
    						// 新しい点
    			for (i1 = 0; i1 < n; i1++)
    				x[i1] += step * 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, snx);
    							// 黄金分割法が失敗
    				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;
    					}
    				}
    			}
    		}
    	}
    
    	if (sw == 0)
    		sw = -1;
    
    	return sw;
    }
    
    /****************************************************************/
    /* 黄金分割法(与えられた方向での最小値)                         */
    /*      a,b : 初期区間 a < b                                    */
    /*      eps : 許容誤差                                          */
    /*      val : 関数値                                            */
    /*      ind : 計算状況                                          */
    /*              >= 0 : 正常終了(収束回数)                       */
    /*              = -1 : 収束せず                                 */
    /*      max : 最大試行回数                                      */
    /*      w : 位置                                                */
    /*      dw : 傾きの成分                                         */
    /*      snx : 関数値を計算する関数の名前                        */
    /*      return : 結果(w+y*dwのy)                              */
    /****************************************************************/
    #include <math.h>
    
    double gold(double a, double b, double eps, double *val, int *ind, int max,
                double *w, double *dw, double (*snx)(double, double *, double *))
    {
    	double f1, f2, fa, fb, tau = (sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
    	int count = 0;
    					// 初期設定
    	*ind  = -1;
    	x1    = b - tau * (b - a);
    	x2    = a + tau * (b - a);
    	f1    = snx(x1, w, dw);
    	f2    = snx(x2, w, dw);
    					// 計算
    	while (count < max && *ind < 0) {
    		count += 1;
    		if (f2 > f1) {
    			if (fabs(b-a) < eps) {
    				*ind = 0;
    				x    = x1;
    				*val = f1;
    			}
    			else {
    				b  = x2;
    				x2 = x1;
    				x1 = a + (1.0 - tau) * (b - a);
    				f2 = f1;
    				f1 = snx(x1, w, dw);
    			}
    		}
    		else {
    			if (fabs(b-a) < eps) {
    				*ind = 0;
    				x    = x2;
    				*val = f2;
    				f1   = f2;
    			}
    			else {
    				a  = x1;
    				x1 = x2;
    				x2 = b - (1.0 - tau) * (b - a);
    				f1 = f2;
    				f2 = snx(x2, w, dw);
    			}
    		}
    	}
    					// 収束した場合の処理
    	if (*ind == 0) {
    		*ind = count;
    		fa   = snx(a, w, dw);
    		fb   = snx(b, w, dw);
    		if (fb < fa) {
    			a  = b;
    			fa = fb;
    		}
    		if (fa < f1) {
    			x    = a;
    			*val = fa;
    		}
    	}
    
    	return x;
    }
    
    /*********************************************/
    /* 与えられた点xから,dx方向へk*dxだけ進んだ */
    /* 点における関数値を計算する                */
    /*********************************************/
    					// 関数1
    double snx1(double k, double *x, double *dx)
    {
    	double x1, y1, y, w[2];
    
    	w[0] = x[0] + k * dx[0];
    	w[1] = x[1] + k * dx[1];
    	x1   = w[0] - 1.0;
    	y1   = w[1] - 2.0;
    	y    = x1 * x1 + y1 * y1;
    
    	return y;
    }
    					// 関数2
    double snx2(double k, double *x, double *dx)
    {
    	double x1, y1, y, w[2];
    
    	w[0] = x[0] + k * dx[0];
    	w[1] = x[1] + k * dx[1];
    	x1   = w[1] - w[0] * w[0];
    	y1   = 1.0 - w[0];
    	y    = 100.0 * x1 * x1 + y1 * y1;
    
    	return y;
    }
    					// 関数3
    double snx3(double k, double *x, double *dx)
    {
    	double x1, y1, z1, y, w[2];
    
    	w[0] = x[0] + k * dx[0];
    	w[1] = x[1] + k * dx[1];
    	x1   = 1.5 - w[0] * (1.0 - w[1]);
    	y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
    	z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
    	y    = x1 * x1 + y1 * y1 + z1 * z1;
    
    	return y;
    }
    
    /********************/
    /* 微係数を計算する */
    /********************/
    					// 関数1
    void dsnx1(double *x, double *dx)
    {
    	dx[0] = -2.0 * (x[0] - 1.0);
    	dx[1] = -2.0 * (x[1] - 2.0);
    }
    					// 関数2
    void dsnx2(double *x, double *dx)
    {
    	dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
    	dx[1] = -200.0 * (x[1] - x[0] * x[0]);
    }
    					// 関数3
    void dsnx3(double *x, double *dx)
    {
    	dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    	dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    }
    			

  2. Java

    /********************************/
    /* 最急降下法による最小値の計算 */
    /*      coded by Y.Suganuma     */
    /********************************/
    import java.io.*;
    import java.util.*;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double eps, 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];
    		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.steep(opt_1, max, n, eps, step, y, x, dx);
    					// 結果の出力
    		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 Steep
    {
    	private int sw;
    					// コンストラクタ
    	Kansu (int s) {sw = s;}
    					// 与えられた点xから,dx方向へk*dxだけ進んだ
    					// 点における関数値を計算する
    	double snx(double k, double x[], double dx[])
    	{
    		double x1, y1, z1, y = 0.0, w[] = new double [2];
    
    		switch (sw) {
    			case 1:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = w[0] - 1.0;
    				y1   = w[1] - 2.0;
    				y    = x1 * x1 + y1 * y1;
    				break;
    			case 2:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = w[1] - w[0] * w[0];
    				y1   = 1.0 - w[0];
    				y    = 100.0 * x1 * x1 + y1 * y1;
    				break;
    			case 3:
    				w[0] = x[0] + k * dx[0];
    				w[1] = x[1] + k * dx[1];
    				x1   = 1.5 - w[0] * (1.0 - w[1]);
    				y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
    				z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
    				y    = x1 * x1 + y1 * y1 + z1 * z1;
    				break;
    		}
    
    		return y;
    	}
    					// 関数の微分
    	void dsnx(double x[], double dx[])
    	{
    		switch (sw) {
    			case 1:
    				dx[0] = -2.0 * (x[0] - 1.0);
    				dx[1] = -2.0 * (x[1] - 2.0);
    				break;
    			case 2:
    				dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
    				dx[1] = -200.0 * (x[1] - x[0] * x[0]);
    				break;
    			case 3:
    				dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                            2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                            2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    				dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                            4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                            6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    				break;
    		}
    	}
    }
    
    abstract class Steep
    {
    	/********************************************************/
    	/* 最急降下法(steepest descent method)                  */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x : x(初期値と答え)                             */
    	/*      dx : 関数の微分値                               */
    	/*      return : >=0 : 正常終了(収束回数)               */
    	/*               =-1 : 収束せず                         */
    	/*               =-2 : 1次元最適化の区間が求まらない   */
    	/*               =-3 : 黄金分割法が失敗                 */
    	/********************************************************/
    
    	abstract double snx(double k, double x[], double dx[]);
    	abstract void dsnx(double x[], double dx[]);
    
    	int steep(int opt_1, int max, int n, double eps, double step, double y[], double x[], double dx[])
    	{
    		double f1, f2, k, sp, y1[] = new double [1], y2[] = new double [1];
    		int count = 0, i1, sw = 0, sw1[] = new int [1];
    
    		y1[0] = snx(0.0, x, dx);
    
    		while (count < max && sw == 0) {
    					// 傾きの計算
    			dsnx(x, dx);
    			count++;
    					// 1次元最適化を行わない
    			if (opt_1 == 0) {
    						// 新しい点
    				for (i1 = 0; i1 < n; i1++)
    					x[i1] += step * 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];
    						}
    					}
    				}
    			}
    		}
    
    		if (sw == 0)
    			sw = -1;
    
    		return sw;
    	}
    
    	/****************************************************************/
    	/* 黄金分割法(与えられた方向での最小値)                         */
    	/*      a,b : 初期区間 a < b                                    */
    	/*      eps : 許容誤差                                          */
    	/*      val : 関数値                                            */
    	/*      ind : 計算状況                                          */
    	/*              >= 0 : 正常終了(収束回数)                       */
    	/*              = -1 : 収束せず                                 */
    	/*      max : 最大試行回数                                      */
    	/*      w : 位置                                                */
    	/*      dw : 傾きの成分                                         */
    	/*      return : 結果(w+y*dwのy)                              */
    	/****************************************************************/
    	double gold(double a, double b, double eps, double val[], int ind[], int max, double w[], double dw[])
    	{
    		double f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
    		int count = 0;
    					// 初期設定
    		ind[0] = -1;
    		x1     = b - tau * (b - a);
    		x2     = a + tau * (b - a);
    		f1     = snx(x1, w, dw);
    		f2     = snx(x2, w, dw);
    					// 計算
    		while (count < max && ind[0] < 0) {
    			count += 1;
    			if (f2 > f1) {
    				if (Math.abs(b-a) < eps) {
    					ind[0] = 0;
    					x      = x1;
    					val[0] = f1;
    				}
    				else {
    					b  = x2;
    					x2 = x1;
    					x1 = a + (1.0 - tau) * (b - a);
    					f2 = f1;
    					f1 = snx(x1, w, dw);
    				}
    			}
    			else {
    				if (Math.abs(b-a) < eps) {
    					ind[0] = 0;
    					x      = x2;
    					val[0] = f2;
    					f1     = f2;
    				}
    				else {
    					a  = x1;
    					x1 = x2;
    					x2 = b - (1.0 - tau) * (b - a);
    					f1 = f2;
    					f2 = snx(x2, w, dw);
    				}
    			}
    		}
    					// 収束した場合の処理
    		if (ind[0] == 0) {
    			ind[0] = count;
    			fa     = snx(a, w, dw);
    			fb     = snx(b, w, dw);
    			if (fb < fa) {
    				a  = b;
    				fa = fb;
    			}
    			if (fa < f1) {
    				x      = a;
    				val[0] = fa;
    			}
    		}
    
    		return x;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,JavaScript の仕様に適合した形で最小値を求めたい式を入力することによって,任意の関数の最小値を画面上で求めることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>最適化(最急降下法)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		str1 = "";
    		str2 = new Array();
    		function main()
    		{
    					// データの設定
    			let n    = parseInt(document.getElementById("n").value);
    			let max  = parseInt(document.getElementById("max").value);
    			let step = parseFloat(document.getElementById("step").value);
    			let x    = document.getElementById("x0").value.split(/ {1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				x[i1] = parseFloat(x[i1]);
    			let opt_1 = 1;
    			if (document.getElementById("line").options[1].selected)
    				opt_1 = 0;
    			str1 = document.getElementById("func").value + ";";
    			str2 = document.getElementById("dfunc").value.split(/\n{1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				str2[i1] = str2[i1] + ";";
    			let eps = 1.0e-10;
    			let y   = new Array();
    			let dx  = new Array();
    			for (let i1 = 0; i1 < n; i1++)
    				dx[i1] = 0.0;
    					// 実行
    			let sw = steep(opt_1, max, n, eps, step, y, x, dx, 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 = " 結果 =";
    				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]);
    		}
    
    		/******************************************************/
    		/* 最急降下法(steepest descent method)                */
    		/*      opt_1 : =0 : 1次元最適化を行わない           */
    		/*              =1 : 1次元最適化を行う               */
    		/*      max : 最大繰り返し回数                        */
    		/*      n : 次元                                      */
    		/*      eps : 収束判定条件                            */
    		/*      step : きざみ幅                               */
    		/*      y : 最小値                                    */
    		/*      x : x(初期値と答え)                           */
    		/*      dx : 関数の微分値                             */
    		/*      fn : 関数値を計算する関数                     */
    		/*      dfn : 関数の微分値を計算する関数              */
    		/*      return : >=0 : 正常終了(収束回数)             */
    		/*               =-1 : 収束せず                       */
    		/*               =-2 : 1次元最適化の区間が求まらない */
    		/*               =-3 : 黄金分割法が失敗               */
    		/******************************************************/
    		function steep(opt_1, max, n, eps, step, y, x, dx, fn, dfn)
    		{
    			let f1, f2, k, sp, count = 0, sw = 0
    			let sw1 = new Array();
    			let y1 = new Array();
    			let y2 = new Array();
    
    			y1[0] = fn(0.0, x, dx);
    
    			while (count < max && sw == 0) {
    					// 傾きの計算
    				dfn(n, x, dx);
    				count++;
    					// 1次元最適化を行わない
    				if (opt_1 == 0) {
    						// 新しい点
    					for (let i1 = 0; i1 < n; i1++)
    						x[i1] += step * 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, fn);
    							// 黄金分割法が失敗
    						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];
    							}
    						}
    					}
    				}
    			}
    
    			if (sw == 0)
    				sw = -1;
    
    			return sw;
    		}
    
    		/******************************************/
    		/* 黄金分割法(与えられた方向での最小値)   */
    		/*      a,b : 初期区間 a < b              */
    		/*      eps : 許容誤差                    */
    		/*      val : 関数値                      */
    		/*      ind : 計算状況                    */
    		/*              >= 0 : 正常終了(収束回数) */
    		/*              = -1 : 収束せず           */
    		/*      max : 最大試行回数                */
    		/*      w : 位置                          */
    		/*      dw : 傾きの成分                   */
    		/*      fn : 関数値を計算する関数         */
    		/*      return : 結果(w+y*dwのy)        */
    		/******************************************/
    		function gold(a, b, eps, val, ind, max, w, dw, fn)
    		{
    			let f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
    			let count = 0;
    					// 初期設定
    			ind[0] = -1;
    			x1     = b - tau * (b - a);
    			x2     = a + tau * (b - a);
    			f1     = fn(x1, w, dw);
    			f2     = fn(x2, w, dw);
    					// 計算
    			while (count < max && ind[0] < 0) {
    				count += 1;
    				if (f2 > f1) {
    					if (Math.abs(b-a) < eps) {
    						ind[0] = 0;
    						x      = x1;
    						val[0] = f1;
    					}
    					else {
    						b  = x2;
    						x2 = x1;
    						x1 = a + (1.0 - tau) * (b - a);
    						f2 = f1;
    						f1 = fn(x1, w, dw);
    					}
    				}
    				else {
    					if (Math.abs(b-a) < eps) {
    						ind[0] = 0;
    						x      = x2;
    						val[0] = f2;
    						f1     = f2;
    					}
    					else {
    						a  = x1;
    						x1 = x2;
    						x2 = b - (1.0 - tau) * (b - a);
    						f1 = f2;
    						f2 = fn(x2, w, dw);
    					}
    				}
    			}
    					// 収束した場合の処理
    			if (ind[0] == 0) {
    				ind[0] = count;
    				fa     = fn(a, w, dw);
    				fb     = fn(b, w, dw);
    				if (fb < fa) {
    					a  = b;
    					fa = fb;
    				}
    				if (fa < f1) {
    					x      = a;
    					val[0] = fa;
    				}
    			}
    
    			return x;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>最適化(最急降下法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,以下に示すような関数の最小値を求める場合に対する値が設定されています(他の問題を実行する場合は,それらを適切に修正してください).
    		<P STYLE="text-align:center"><IMG SRC="steep.gif"></P>
    		<DT>n 変数関数 f(<B>x</B>) に対する最急降下法においては,現在の点 <B>x</B> における傾き,
    		<P STYLE="text-align:center"><IMG SRC="steep1.gif"></P>
    		<DT>を求め,その方向にαだけ進んだ点の関数値 f( <B>x</B> + α<B>p</B>) を計算する,という処理を繰り返します.目的関数の箇所には,この例に対する f( <B>x</B> + α<B>p</B>) の計算式,また,目的関数の微分の箇所には,<B>p</B> の計算式が与えられています.ただし,プログラム内では,変数 x は,この例の x と y からなる配列,変数 dx は <B>p</B> に対応し,関数 f を x 及び y で微分したものからなる配列,また,k はαに対応する変数として表現されています.なお,目的関数及び目的関数の微分は,JavaScript の仕様に適合した形式で記述してあることに注意してください. 
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		変数の数(n):<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="2"> 
    		最大繰り返し回数:<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="10000"> 
    		刻み幅:<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.002"><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="yes"> 行う </OPTION><BR>
    							<OPTION VALUE="no"> 行わない </OPTION><BR>
    					  </SELECT><BR><BR>
    		目的関数: f(<B>x</B>+α<B>p</B>) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="70" VALUE="100.0 * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) + (1.0 - (x[0] + k * dx[0])) * (1.0 - (x[0] + k * dx[0]))"><BR><BR>
    		目的関数の微分: f'(x) = <TEXTAREA ID="dfunc" STYLE="font-size: 110%" COLS="70" ROWS="10">
    400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
    -200.0 * (x[1] - x[0] * x[0])
    		</TEXTAREA><BR><BR>
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
    		結果:<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="40" ROWS="2"></TEXTAREA>
    	<DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /********************************/
    /* 最急降下法による最小値の計算 */
    /*      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);
    
    	$str = trim(fgets(STDIN));
    	strtok($str, " ");
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$x[$i1] = floatval(strtok(" "));
    					// 実行
    	switch ($fun) {
    		case 1:
    			$sw = steep($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx1", "dsnx1");
    			break;
    		case 2:
    			$sw = steep($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx2", "dsnx2");
    			break;
    		case 3:
    			$sw = steep($opt_1, $max, $n, $eps, $step, $y, $x, $dx, "snx3", "dsnx3");
    			break;
    	}
    					// 結果の出力
    	if ($sw < 0) {
    		printf("   収束しませんでした!");
    		switch ($sw) {
    			case -1:
    				printf("(収束回数)\n");
    				break;
    			case -2:
    				printf("(1次元最適化の区間)\n");
    				break;
    			case -3:
    				printf("(黄金分割法)\n");
    				break;
    		}
    	}
    	else {
    		printf("   結果=");
    		for ($i1 = 0; $i1 < $n; $i1++)
    			printf("%f ", $x[$i1]);
    		printf(" 最小値=%f  回数=%d\n", $y, $sw);
    	}
    
    /********************************************************/
    /* 最急降下法(steepest descent method)                  */
    /*      opt_1 : =0 : 1次元最適化を行わない             */
    /*              =1 : 1次元最適化を行う                 */
    /*      max : 最大繰り返し回数                          */
    /*      n : 次元                                        */
    /*      eps : 収束判定条件                              */
    /*      step : きざみ幅                                 */
    /*      y : 最小値                                      */
    /*      x : x(初期値と答え)                             */
    /*      dx : 関数の微分値                               */
    /*      snx : 関数値を計算する関数名                    */
    /*      dsnx : 関数の微分を計算する関数名(符号を変える) */
    /*      return : >=0 : 正常終了(収束回数)               */
    /*               =-1 : 収束せず                         */
    /*               =-2 : 1次元最適化の区間が求まらない   */
    /*               =-3 : 黄金分割法が失敗                 */
    /********************************************************/
    function steep($opt_1, $max, $n, $eps, $step, &$y, &$x, $dx, $snx, $dsnx)
    {
    	$count = 0;
    	$sw    = 0;
    	$y1    = $snx(0.0, $x, $dx);
    
    	while ($count < $max && $sw == 0) {
    					// 傾きの計算
    		$dsnx($x, $dx);
    		$count++;
    					// 1次元最適化を行わない
    		if ($opt_1 == 0) {
    						// 新しい点
    			for ($i1 = 0; $i1 < $n; $i1++)
    				$x[$i1] += $step * $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;
    					}
    				}
    			}
    		}
    	}
    
    	if ($sw == 0)
    		$sw = -1;
    
    	return $sw;
    }
    
    /*********************************************/
    /* 与えられた点xから,dx方向へk*dxだけ進んだ */
    /* 点における関数値を計算する                */
    /*********************************************/
    					// 関数1
    function snx1($k, $x, $dx)
    {
    	$w  = array(2);
    
    	$w[0] = $x[0] + $k * $dx[0];
    	$w[1] = $x[1] + $k * $dx[1];
    	$x1   = $w[0] - 1.0;
    	$y1   = $w[1] - 2.0;
    	$y    = $x1 * $x1 + $y1 * $y1;
    
    	return $y;
    }
    					// 関数2
    function snx2($k, $x, $dx)
    {
    	$w  = array(2);
    
    	$w[0] = $x[0] + $k * $dx[0];
    	$w[1] = $x[1] + $k * $dx[1];
    	$x1   = $w[1] - $w[0] * $w[0];
    	$y1   = 1.0 - $w[0];
    	$y    = 100.0 * $x1 * $x1 + $y1 * $y1;
    
    	return $y;
    }
    					// 関数3
    function snx3($k, $x, $dx)
    {
    	$w = array(2);
    
    	$w[0] = $x[0] + $k * $dx[0];
    	$w[1] = $x[1] + $k * $dx[1];
    	$x1   = 1.5 - $w[0] * (1.0 - $w[1]);
    	$y1   = 2.25 - $w[0] * (1.0 - $w[1] * $w[1]);
    	$z1   = 2.625 - $w[0] * (1.0 - $w[1] * $w[1] * $w[1]);
    	$y    = $x1 * $x1 + $y1 * $y1 + $z1 * $z1;
    
    	return $y;
    }
    
    /********************/
    /* 微係数を計算する */
    /********************/
    					// 関数1
    function dsnx1($x, &$dx)
    {
    	$dx[0] = -2.0 * ($x[0] - 1.0);
    	$dx[1] = -2.0 * ($x[1] - 2.0);
    }
    					// 関数2
    function dsnx2($x, &$dx)
    {
    	$dx[0] = 400.0 * $x[0] * ($x[1] - $x[0] * $x[0]) + 2.0 * (1.0 - $x[0]);
    	$dx[1] = -200.0 * ($x[1] - $x[0] * $x[0]);
    }
    					// 関数3
    function dsnx3($x, &$dx)
    {
    	$dx[0] = 2.0 * (1.0 - $x[1]) * (1.5 - $x[0] * (1.0 - $x[1])) +
                2.0 * (1.0 - $x[1] * $x[1]) * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) +
                2.0 * (1.0 - $x[1] * $x[1] * $x[1]) * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
    	$dx[1] = -2.0 * $x[0] * (1.5 - $x[0] * (1.0 - $x[1])) -
                4.0 * $x[0] * $x[1] * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) -
                6.0 * $x[0] * $x[1] * $x[1] * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
    }
    
    /****************************************************************/
    /* 黄金分割法(与えられた方向での最小値)                         */
    /*      a,b : 初期区間 a < b                                    */
    /*      eps : 許容誤差                                          */
    /*      val : 関数値                                            */
    /*      ind : 計算状況                                          */
    /*              >= 0 : 正常終了(収束回数)                       */
    /*              = -1 : 収束せず                                 */
    /*      max : 最大試行回数                                      */
    /*      w : 位置                                                */
    /*      dw : 傾きの成分                                         */
    /*      snx : 関数値を計算する関数の名前                        */
    /*      return : 結果(w+y*dwのy)                              */
    /****************************************************************/
    function gold($a, $b, $eps, &$val, &$ind, $max, $w, $dw, $snx)
    {
    					// 初期設定
    	$tau   = (sqrt(5.0) - 1.0) / 2.0;
    	$x     = 0.0;
    	$count = 0;
    	$ind   = -1;
    	$x1    = $b - $tau * ($b - $a);
    	$x2    = $a + $tau * ($b - $a);
    	$f1    = $snx($x1, $w, $dw);
    	$f2    = $snx($x2, $w, $dw);
    					// 計算
    	while ($count < $max && $ind < 0) {
    		$count += 1;
    		if ($f2 > $f1) {
    			if (abs($b-$a) < $eps) {
    				$ind = 0;
    				$x   = $x1;
    				$val = $f1;
    			}
    			else {
    				$b  = $x2;
    				$x2 = $x1;
    				$x1 = $a + (1.0 - $tau) * ($b - $a);
    				$f2 = $f1;
    				$f1 = $snx($x1, $w, $dw);
    			}
    		}
    		else {
    			if (abs($b-$a) < $eps) {
    				$ind = 0;
    				$x   = $x2;
    				$val = $f2;
    				$f1  = $f2;
    			}
    			else {
    				$a = $x1;
    				$x1 = $x2;
    				$x2 = $b - (1.0 - $tau) * ($b - $a);
    				$f1 = $f2;
    				$f2 = $snx($x2, $w, $dw);
    			}
    		}
    	}
    					// 収束した場合の処理
    	if ($ind == 0) {
    		$ind = $count;
    		$fa  = $snx($a, $w, $dw);
    		$fb  = $snx($b, $w, $dw);
    		if ($fb < $fa) {
    			$a  = $b;
    			$fa = $fb;
    		}
    		if ($fa < $f1) {
    			$x   = $a;
    			$val = $fa;
    		}
    	}
    
    	return $x;
    }
    
    ?>
    			

  5. Ruby

    #*******************************/
    # 最急降下法による最小値の計算 */
    #      coded by Y.Suganuma     */
    #*******************************/
    
    #********************************************/
    # 与えられた点xから,dx方向へk*dxだけ進んだ */
    # 点における関数値及び微係数を計算する      */
    #********************************************/
    					# 関数1
    snx1 = Proc.new { |sw, k, x, dx|
    	if sw == 0
    		w    = Array.new(2)
    		w[0] = x[0] + k * dx[0]
    		w[1] = x[1] + k * dx[1]
    		x1   = w[0] - 1.0
    		y1   = w[1] - 2.0
    		x1 * x1 + y1 * y1
    	else
    		dx[0] = -2.0 * (x[0] - 1.0)
    		dx[1] = -2.0 * (x[1] - 2.0)
    	end
    }
    					# 関数2
    snx2 = Proc.new { |sw, k, x, dx|
    	if sw == 0
    		w    = Array.new(2)
    		w[0] = x[0] + k * dx[0]
    		w[1] = x[1] + k * dx[1]
    		x1   = w[1] - w[0] * w[0]
    		y1   = 1.0 - w[0]
    		100.0 * x1 * x1 + y1 * y1
    	else
    		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
    		dx[1] = -200.0 * (x[1] - x[0] * x[0])
    	end
    }
    					# 関数3
    snx3 = Proc.new { |sw, k, x, dx|
    	if sw == 0
    		w    = Array.new(2)
    		w[0] = x[0] + k * dx[0]
    		w[1] = x[1] + k * dx[1]
    		x1   = 1.5 - w[0] * (1.0 - w[1])
    		y1   = 2.25 - w[0] * (1.0 - w[1] * w[1])
    		z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
    		x1 * x1 + y1 * y1 + z1 * z1
    	else
    		dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                    2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                    2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
    		dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                    4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                    6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
    	end
    }
    
    #***************************************************************/
    # 黄金分割法(与えられた方向での最小値)                         */
    #      a,b : 初期区間 a < b                                    */
    #      eps : 許容誤差                                          */
    #      val : 関数値                                            */
    #      ind : 計算状況                                          */
    #              >= 0 : 正常終了(収束回数)                       */
    #              = -1 : 収束せず                                 */
    #      max : 最大試行回数                                      */
    #      w : 位置                                                */
    #      dw : 傾きの成分                                         */
    #      snx : 関数値を計算する関数の名前                        */
    #      return : 結果(w+y*dwのy)                              */
    #***************************************************************/
    def gold(a, b, eps, val, ind, max, w, dw, &snx)
    					# 初期設定
    	tau    = (Math.sqrt(5.0) - 1.0) / 2.0
    	x      = 0.0
    	count  = 0
    	ind[0] = -1
    	x1     = b - tau * (b - a)
    	x2     = a + tau * (b - a)
    	f1     = snx.call(0, x1, w, dw)
    	f2     = snx.call(0, x2, w, dw)
    					# 計算
    	while count < max && ind[0] < 0
    		count += 1
    		if f2 > f1
    			if (b-a).abs() < eps
    				ind[0] = 0
    				x      = x1
    				val[0] = f1
    			else
    				b  = x2
    				x2 = x1
    				x1 = a + (1.0 - tau) * (b - a)
    				f2 = f1
    				f1 = snx.call(0, x1, w, dw)
    			end
    		else
    			if (b-a).abs() < eps
    				ind[0] = 0
    				x      = x2
    				val[0] = f2
    				f1     = f2
    			else
    				a  = x1
    				x1 = x2
    				x2 = b - (1.0 - tau) * (b - a)
    				f1 = f2
    				f2 = snx.call(0, x2, w, dw)
    			end
    		end
    	end
    					# 収束した場合の処理
    	if ind[0] == 0
    		ind[0] = count
    		fa     = snx.call(0, a, w, dw)
    		fb     = snx.call(0, b, w, dw)
    		if fb < fa
    			a  = b
    			fa = fb
    		end
    		if fa < f1
    			x      = a
    			val[0] = fa
    		end
    	end
    
    	return x
    end
    
    #*******************************************************/
    # 最急降下法(steepest descent method)                  */
    #      opt_1 : =0 : 1次元最適化を行わない             */
    #              =1 : 1次元最適化を行う                 */
    #      max : 最大繰り返し回数                          */
    #      n : 次元                                        */
    #      eps : 収束判定条件                              */
    #      step : きざみ幅                                 */
    #      y : 最小値                                      */
    #      x : x(初期値と答え)                             */
    #      dx : 関数の微分値                               */
    #      snx : 関数値とその微分を計算する関数名          */
    #            (微分は,その符号を変える)               */
    #      return : >=0 : 正常終了(収束回数)               */
    #               =-1 : 収束せず                         */
    #               =-2 : 1次元最適化の区間が求まらない   */
    #               =-3 : 黄金分割法が失敗                 */
    #*******************************************************/
    def steep(opt_1, max, n, eps, step, y, x, dx, &snx)
    
    	count = 0
    	sw    = 0
    	y2    = Array.new(1)
    	sw1   = Array.new(1)
    
    	y1 = snx.call(0, 0.0, x, dx)
    
    	while count < max && sw == 0
    					# 傾きの計算
    		snx.call(1, 0.0, x, dx)
    		count += 1
    					# 1次元最適化を行わない
    		if opt_1 == 0
    						# 新しい点
    			for i1 in 0 ... n
    				x[i1] += step * dx[i1]
    			end
    						# 新しい関数値
    			y2[0] = snx.call(0, 0.0, x, dx)
    							# 関数値の変化が大きい
    			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)
    			if f2 > f1
    				sp = -step
    			end
    			for i1 in 0 ... max
    				f2 = snx.call(0, sp, x, dx)
    				if f2 > f1
    					sw1[0] = 1
    					break
    				else
    					sp *= 2.0
    					f1  = f2
    				end
    			end
    						# 区間が求まらない
    			if sw1[0] == 0
    				sw = -2
    						# 区間が求まった
    			else
    							# 黄金分割法
    				k = gold(0.0, sp, eps, y2, sw1, max, x, dx, &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
    	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(" ");
    for i1 in 0 ... n
    	x[i1]  = Float(a[i1+1])
    	dx[i1] = 0.0
    end
    				# 実行
    case fun
    	when 1
    		sw = steep(opt_1, max, n, eps, step, y, x, dx, &snx1)
    	when 2
    		sw = steep(opt_1, max, n, eps, step, y, x, dx, &snx2)
    	when 3
    		sw = steep(opt_1, max, n, eps, step, y, x, dx, &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
    
    ############################################
    # 最急降下法(steepest descent method)
    #      opt_1 : =0 : 1次元最適化を行わない
    #              =1 : 1次元最適化を行う
    #      max : 最大繰り返し回数
    #      n : 次元
    #      eps : 収束判定条件
    #      step : きざみ幅
    #      y : 最小値
    #      x : x(初期値と答え)
    #      dx : 関数の微分値
    #      snx : 関数値を計算する関数名
    #      dsnx : 関数の微分を計算する関数名(符号を変える)
    #      return : >=0 : 正常終了(収束回数)
    #               =-1 : 収束せず
    #               =-2 : 1次元最適化の区間が求まらない
    #               =-3 : 黄金分割法が失敗
    #      coded by Y.Suganuma
    ############################################
    
    def steep(opt_1, max, n, eps, step, y, x, dx, snx, dsnx) :
    
    	count = 0
    	sw    = 0
    	y1    = snx(0.0, x, dx)
    	y2    = np.empty(1, np.float)
    	sw1   = np.empty(1, np.int)
    
    	while count < max and sw == 0 :
    				# 傾きの計算
    		dsnx(x, dx)
    		count += 1
    				# 1次元最適化を行わない
    		if opt_1 == 0 :
    					# 新しい点
    			for i1 in range(0, n) :
    				x[i1] += step * 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]
    
    	if sw == 0 :
    		sw = -1
    
    	return sw
    
    ############################################
    # 最急降下法による最小値の計算
    #      coded by Y.Suganuma
    ############################################
    
    ############################################
    # 与えられた点xから,dx方向へk*dxだけ進んだ
    # 点における関数値を計算する
    ############################################
    			 # 関数1
    def snx1(k, x, dx) :
    	w    = np.empty(2, np.float)
    	w[0] = x[0] + k * dx[0]
    	w[1] = x[1] + k * dx[1]
    	x1   = w[0] - 1.0
    	y1   = w[1] - 2.0
    	y    = x1 * x1 + y1 * y1
    	return y
    			# 関数2
    def snx2(k, x, dx) :
    	w    = np.empty(2, np.float)
    	w[0] = x[0] + k * dx[0]
    	w[1] = x[1] + k * dx[1]
    	x1   = w[1] - w[0] * w[0]
    	y1   = 1.0 - w[0]
    	y    = 100.0 * x1 * x1 + y1 * y1
    	return y
    			# 関数3
    def snx3(k, x, dx) :
    	w    = np.empty(2, np.float)
    	w[0] = x[0] + k * dx[0]
    	w[1] = x[1] + k * dx[1]
    	x1   = 1.5 - w[0] * (1.0 - w[1])
    	y1   = 2.25 - w[0] * (1.0 - w[1] * w[1])
    	z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
    	y    = x1 * x1 + y1 * y1 + z1 * z1
    	return y
    
    ############################################
    # 微係数を計算する
    ############################################
    			# 関数1
    def dsnx1(x, dx) :
    	dx[0] = -2.0 * (x[0] - 1.0)
    	dx[1] = -2.0 * (x[1] - 2.0)
    			# 関数2
    def dsnx2(x, dx) :
    	dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
    	dx[1] = -200.0 * (x[1] - x[0] * x[0])
    			# 関数3
    def dsnx3(x, dx) :
    	dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) + 2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) + 2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
    	dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) - 4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) - 6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
    
    			# データの入力
    sw    = 0
    line  = sys.stdin.readline()
    s     = line.split()
    fun   = int(s[1])
    n     = int(s[3])
    max   = int(s[5])
    opt_1 = int(s[7])
    
    line  = sys.stdin.readline()
    s     = line.split()
    eps   = float(s[1])
    step  = float(s[3])
    
    x     = np.empty(n, np.float)
    dx    = np.empty(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 = steep(opt_1, max, n, eps, step, y, x, dx, snx1, dsnx1)
    elif fun == 2 :
    	sw = steep(opt_1, max, n, eps, step, y, x, dx, snx2, dsnx2)
    else :
    	sw = steep(opt_1, max, n, eps, step, y, x, dx, snx3, dsnx3)
    			# 結果の出力
    if sw < 0 :
    	print("   収束しませんでした!", end="")
    	if sw == -1 :
    		print("(収束回数)")
    	elif sw == -2 :
    		print("(1次元最適化の区間)")
    	else :
    		print("(黄金分割法)")
    else :
    	print("   結果=", end="")
    	for i1 in range(0, n) :
    		print(str(x[i1]) + " ", end="")
    	print(" 最小値=" + str(y[0]) + "  回数=" + str(sw))
    			

  7. C#

    /********************************/
    /* 最急降下法による最小値の計算 */
    /*      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[] 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]);
    					// 実行
    		int sw;
    		if (fun == 1)
    			sw = steep(opt_1, max, n, eps, step, ref y, x, dx, snx1, dsnx1);
    		else if (fun == 2)
    			sw = steep(opt_1, max, n, eps, step, ref y, x, dx, snx2, dsnx2);
    		else
    			sw = steep(opt_1, max, n, eps, step, ref y, x, dx, snx3, dsnx3);
    					// 結果の出力
    		if (sw < 0) {
    			Console.Write("   収束しませんでした!");
    			switch (sw) {
    				case -1:
    					Console.WriteLine("(収束回数)");
    					break;
    				case -2:
    					Console.WriteLine("(1次元最適化の区間)");
    					break;
    				case -3:
    					Console.WriteLine("(黄金分割法)");
    					break;
    			}
    		}
    		else {
    			Console.Write("   結果=");
    			for (int i1 = 0; i1 < n; i1++)
    				Console.Write(x[i1] + " ");
    			Console.WriteLine(" 最小値=" + y + "  回数=" + sw);
    		}
    	}
    			// 関数1の値の計算
    	double snx1(double k, double[] x, double[] dx)
    	{
    		double[] w = new double [2];
    		w[0]       = x[0] + k * dx[0];
    		w[1]       = x[1] + k * dx[1];
    		double x1  = w[0] - 1.0;
    		double y1  = w[1] - 2.0;
    		double y   = x1 * x1 + y1 * y1;
    		return y;
    	}
    			// 関数2の値の計算
    	double snx2(double k, double[] x, double[] dx)
    	{
    		double[] w = new double [2];
    		w[0]       = x[0] + k * dx[0];
    		w[1]       = x[1] + k * dx[1];
    		double x1  = w[1] - w[0] * w[0];
    		double y1  = 1.0 - w[0];
    		double y   = 100.0 * x1 * x1 + y1 * y1;
    		return y;
    	}
    			// 関数3の値の計算
    	double snx3(double k, double[] x, double[] dx)
    	{
    		double[] w = new double [2];
    		w[0]       = x[0] + k * dx[0];
    		w[1]       = x[1] + k * dx[1];
    		double x1  = 1.5 - w[0] * (1.0 - w[1]);
    		double y1  = 2.25 - w[0] * (1.0 - w[1] * w[1]);
    		double z1  = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
    		double y   = x1 * x1 + y1 * y1 + z1 * z1;
    		return y;
    	}
    					// 関数1の微分
    	int dsnx1(double[] x, double[] dx)
    	{
    		dx[0] = -2.0 * (x[0] - 1.0);
    		dx[1] = -2.0 * (x[1] - 2.0);
    		return 0;
    	}
    					// 関数2の微分
    	int dsnx2(double[] x, double[] dx)
    	{
    		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
    		dx[1] = -200.0 * (x[1] - x[0] * x[0]);
    		return 0;
    	}
    					// 関数3の微分
    	int dsnx3(double[] x, double[] dx)
    	{
    		dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                    2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                    2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    		dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                    4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                    6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
    		return 0;
    	}
    
    	/********************************************************/
    	/* 最急降下法(steepest descent method)                  */
    	/*      opt_1 : =0 : 1次元最適化を行わない             */
    	/*              =1 : 1次元最適化を行う                 */
    	/*      max : 最大繰り返し回数                          */
    	/*      n : 次元                                        */
    	/*      eps : 収束判定条件                              */
    	/*      step : きざみ幅                                 */
    	/*      y : 最小値                                      */
    	/*      x : x(初期値と答え)                             */
    	/*      dx : 関数の微分値                               */
    	/*      fn : 関数値を計算する関数                       */
    	/*      dfn : 関数の微分値を計算する関数                */
    	/*      return : >=0 : 正常終了(収束回数)               */
    	/*               =-1 : 収束せず                         */
    	/*               =-2 : 1次元最適化の区間が求まらない   */
    	/*               =-3 : 黄金分割法が失敗                 */
    	/********************************************************/
    	int steep(int opt_1, int max, int n, double eps, double step, ref double y,
    	          double[] x, double[] dx, Func<double, double[], double[], double> fn,
    	          Func<double[], double[], int> dfn)
    	{
    		int count = 0;
    		int sw    = 0;
    		double y1 = fn(0.0, x, dx);
    		double y2 = 0.0;
    
    		while (count < max && sw == 0) {
    					// 傾きの計算
    			dfn(x, dx);
    			count++;
    					// 1次元最適化を行わない
    			if (opt_1 == 0) {
    						// 新しい点
    				for (int i1 = 0; i1 < n; i1++)
    					x[i1] += step * dx[i1];
    						// 新しい関数値
    				y2 = fn(0.0, x, dx);
    							// 関数値の変化が大きい
    				if (Math.Abs(y2-y1) > eps)
    					y1 = y2;
    							// 収束(関数値の変化<eps)
    				else {
    					sw = count;
    					y  = y2;
    				}
    			}
    					// 1次元最適化を行う
    			else {
    						// 区間を決める
    				int 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;
    						}
    					}
    				}
    			}
    		}
    
    		if (sw == 0)
    			sw = -1;
    
    		return sw;
    	}
    
    	/******************************************/
    	/* 黄金分割法(与えられた方向での最小値)   */
    	/*      a,b : 初期区間 a < b              */
    	/*      eps : 許容誤差                    */
    	/*      val : 関数値                      */
    	/*      ind : 計算状況                    */
    	/*              >= 0 : 正常終了(収束回数) */
    	/*              = -1 : 収束せず           */
    	/*      max : 最大試行回数                */
    	/*      w : 位置                          */
    	/*      dw : 傾きの成分                   */
    	/*      fn : 関数値を計算する関数         */
    	/*      return : 結果(w+y*dwのy)        */
    	/******************************************/
    	double gold(double a, double b, double eps, ref double val, ref int ind,
    	            int max, double[] w, double[] dw,
    	            Func<double, double[], double[], double> fn)
    	{
    					// 初期設定
    		double tau = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0;
    		double x1  = b - tau * (b - a);
    		double x2  = a + tau * (b - a);
    		double f1  = fn(x1, w, dw);
    		double f2  = fn(x2, w, dw);
    		int count  = 0;
    		ind        = -1;
    					// 計算
    		while (count < max && ind < 0) {
    			count += 1;
    			if (f2 > f1) {
    				if (Math.Abs(b-a) < eps) {
    					ind = 0;
    					x   = x1;
    					val = f1;
    				}
    				else {
    					b  = x2;
    					x2 = x1;
    					x1 = a + (1.0 - tau) * (b - a);
    					f2 = f1;
    					f1 = fn(x1, w, dw);
    				}
    			}
    			else {
    				if (Math.Abs(b-a) < eps) {
    					ind = 0;
    					x   = x2;
    					val = f2;
    					f1  = f2;
    				}
    				else {
    					a  = x1;
    					x1 = x2;
    					x2 = b - (1.0 - tau) * (b - a);
    					f1 = f2;
    					f2 = fn(x2, w, dw);
    				}
    			}
    		}
    					// 収束した場合の処理
    		if (ind == 0) {
    			ind       = count;
    			double fa = fn(a, w, dw);
    			double fb = fn(b, w, dw);
    			if (fb < fa) {
    				a  = b;
    				fa = fb;
    			}
    			if (fa < f1) {
    				x   = a;
    				val = fa;
    			}
    		}
    
    		return x;
    	}
    }
    			

  8. VB

    ''''''''''''''''''''''''''''''''
    ' 最急降下法による最小値の計算 '
    '      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 x(n) As Double
    		Dim dx(n) As Double
    		Dim y As Double = 0.0
    		str             = MS.Split(Console.ReadLine().Trim())
    		For i1 As Integer = 0 To n-1
    			x(i1) = Double.Parse(str(i1+1))
    		Next
    					' 関数1の値の計算(ラムダ式)
    		Dim snx1 = Function(k0, x0, dx0) As Double
    			Dim w(2) As Double
    			w(0) = x0(0) + k0 * dx0(0)
    			w(1) = x0(1) + k0 * dx0(1)
    			Dim x1 As Double = w(0) - 1.0
    			Dim y1 As Double = w(1) - 2.0
    			Dim y0 As Double = x1 * x1 + y1 * y1
    			Return y0
    		End Function
    					' 関数2の値の計算(ラムダ式)
    		Dim snx2 = Function(k0, x0, dx0) As Double
    			Dim w(2) As Double
    			w(0) = x0(0) + k0 * dx0(0)
    			w(1) = x0(1) + k0 * dx0(1)
    			Dim x1 As Double = w(1) - w(0) * w(0)
    			Dim y1 As Double = 1.0 - w(0)
    			Dim y0 As Double = 100.0 * x1 * x1 + y1 * y1
    			Return y0
    		End Function
    					' 関数3の値の計算(ラムダ式)
    		Dim snx3 = Function(k0, x0, dx0) As Double
    			Dim w(2) As Double
    			w(0) = x0(0) + k0 * dx0(0)
    			w(1) = x0(1) + k0 * dx0(1)
    			Dim x1 As Double = 1.5 - w(0) * (1.0 - w(1))
    			Dim y1 As Double = 2.25 - w(0) * (1.0 - w(1) * w(1))
    			Dim z1 As Double = 2.625 - w(0) * (1.0 - w(1) * w(1) * w(1))
    			Dim y0 As Double = x1 * x1 + y1 * y1 + z1 * z1
    			Return y0
    		End Function
    					' 関数1の微分(ラムダ式)
    		Dim dsnx1 = Function(x0, dx0) As Integer
    			dx0(0) = -2.0 * (x0(0) - 1.0)
    			dx0(1) = -2.0 * (x0(1) - 2.0)
    			Return 0
    		End Function
    					' 関数2の微分(ラムダ式)
    		Dim dsnx2 = Function(x0, dx0) As Integer
    			dx0(0) = 400.0 * x0(0) * (x0(1) - x0(0) * x0(0)) + 2.0 * (1.0 - x0(0))
    			dx0(1) = -200.0 * (x0(1) - x0(0) * x0(0))
    			Return 0
    		End Function
    					' 関数3の微分(ラムダ式)
    		Dim dsnx3 = Function(x0, dx0) As Integer
    			dx0(0) = 2.0 * (1.0 - x0(1)) * (1.5 - x0(0) * (1.0 - x0(1))) +
    	                 2.0 * (1.0 - x0(1) * x0(1)) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) +
    	                 2.0 * (1.0 - x0(1) * x0(1) * x0(1)) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
    			dx0(1) = -2.0 * x0(0) * (1.5 - x0(0) * (1.0 - x0(1))) -
    	                 4.0 * x0(0) * x0(1) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) -
    	                 6.0 * x0(0) * x0(1) * x0(1) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
    			Return 0
    		End Function
    					' 実行
    		Dim sw As Integer
    		If fun = 1
    			sw = steep(opt_1, max, n, eps, step1, y, x, dx, snx1, dsnx1)
    		ElseIf fun = 2
    			sw = steep(opt_1, max, n, eps, step1, y, x, dx, snx2, dsnx2)
    		Else
    			sw = steep(opt_1, max, n, eps, step1, y, x, dx, snx3, dsnx3)
    		End If
    					' 結果の出力
    		If sw < 0
    			Console.Write("   収束しませんでした!")
    			If sw = -1
    				Console.WriteLine("(収束回数)")
    			ElseIf sw = -2
    				Console.WriteLine("(1次元最適化の区間)")
    			ElseIf sw = -3
    				Console.WriteLine("(黄金分割法)")
    			End If
    		Else
    			Console.Write("   結果=")
    			For i1 As Integer = 0 To n-1
    				Console.Write(x(i1) & " ")
    			Next
    			Console.WriteLine(" 最小値=" & y & "  回数=" & sw)
    		End If
    
    	End Sub
    
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' 最急降下法(steepest descent method)                  '
    	'      opt_1 : =0 : 1次元最適化を行わない             '
    	'              =1 : 1次元最適化を行う                 '
    	'      max : 最大繰り返し回数                          '
    	'      n : 次元                                        '
    	'      eps : 収束判定条件                              '
    	'      step1 : きざみ幅                                 '
    	'      y : 最小値                                     '
    	'      x : x(初期値と答え)                             '
    	'      dx : 関数の微分値                               '
    	'      fn : 関数値を計算する関数                       '
    	'      dfn : 関数の微分値を計算する関数                '
    	'      return : >=0 : 正常終了(収束回数)               '
    	'               =-1 : 収束せず                         '
    	'               =-2 : 1次元最適化の区間が求まらない   '
    	'               =-3 : 黄金分割法が失敗                 '
    	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function steep(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,
    	               fn As Func(Of Double, Double(), Double(), Double),
    	               dfn As Func(Of Double(), Double(), Integer))
    
    		Dim count As Integer = 0
    		Dim sw As Integer    = 0
    		Dim y1 As Double     = fn(0.0, x, dx)
    		Dim y2 As Double     = 0.0
    
    		Do While count < max and sw = 0
    					' 傾きの計算
    			dfn(x, dx)
    			count += 1
    					' 1次元最適化を行わない
    			If opt_1 = 0
    						' 新しい点
    				For i1 As Integer = 0 To n-1
    					x(i1) += step1 * 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
    						' 区間を決める
    				Dim sw1 As Integer = 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
    		Loop
    
    		If sw = 0
    			sw = -1
    		End If
    
    		Return sw
    
    	End Function
    
    	''''''''''''''''''''''''''''''''''''''''''
    	' 黄金分割法(与えられた方向での最小値)   '
    	'      a,b : 初期区間 a < b              '
    	'      eps : 許容誤差                    '
    	'      val : 関数値                      '
    	'      ind : 計算状況                    '
    	'              >= 0 : 正常終了(収束回数) '
    	'              = -1 : 収束せず           '
    	'      max : 最大試行回数                '
    	'      w : 位置                          '
    	'      dw : 傾きの成分                   '
    	'      fn : 関数値を計算する関数         '
    	'      return : 結果(w+y*dwのy)        '
    	''''''''''''''''''''''''''''''''''''''''''
    	Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
    	              ByRef ind As Integer, max As Integer, w() As Double, dw() As Double,
    	              fn As Func(Of Double, Double(), Double(), Double))
    
    					' 初期設定
    		Dim tau As Double    = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0
    		Dim x1 As Double     = b - tau * (b - a)
    		Dim x2 As Double     = a + tau * (b - a)
    		Dim f1 As Double     = fn(x1, w, dw)
    		Dim f2 As Double     = fn(x2, w, dw)
    		Dim count As Integer = 0
    		ind = -1
    					' 計算
    		Do While count < max and ind < 0
    			count += 1
    			If f2 > f1
    				If Math.Abs(b-a) < eps
    					ind = 0
    					x   = x1
    					val = f1
    				Else
    					b  = x2
    					x2 = x1
    					x1 = a + (1.0 - tau) * (b - a)
    					f2 = f1
    					f1 = fn(x1, w, dw)
    				End If
    			Else
    				If Math.Abs(b-a) < eps
    					ind = 0
    					x   = x2
    					val = f2
    					f1  = f2
    				Else
    					a  = x1
    					x1 = x2
    					x2 = b - (1.0 - tau) * (b - a)
    					f1 = f2
    					f2 = fn(x2, w, dw)
    				End If
    			End If
    		Loop
    					' 収束した場合の処理
    		If ind = 0
    			ind = count
    			Dim fa As Double = fn(a, w, dw)
    			Dim fb As Double = fn(b, w, dw)
    			If fb < fa
    				a  = b
    				fa = fb
    			End If
    			If fa < f1
    				x   = a
    				val = fa
    			End If
    		End If
    
    		Return x
    
    	End Function
    
    End Module
    			

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