情報学部 菅沼ホーム 目次 索引

最適化（黄金分割法）

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

　　プログラムは，f(x) = x4 + 3x3 + 2x2 + 1 の最小値を黄金分割法で求めた例です．

  1. C++ 

    /*********************************************/
/* 黄金分割法によるy=x^4+3x^3+2x^2+1の最小値 */
/*      coded by Y.Suganuma                  */
/*********************************************/
#include <stdio.h>

double snx(double);
double gold(double, double, double, double *, int *, int, double (*)(double));

int main()
{
	double a, b, eps, val, x;
	int ind, max;

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

	x = gold(a, b, eps, &val, &ind, max, snx);

	printf("x %f val %f ind %d\n", x, val, ind);

	return 0;
}

/****************/
/* 関数値の計算 */
/****************/
double snx(double x)
{
	double f;

	f  = x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0;

	return f;
}

/******************************************/
/* 黄金分割法(関数の最小値)               */
/*      a,b : 初期区間 a < b              */
/*      eps : 許容誤差                    */
/*      val : 間数値                      */
/*      ind : 計算状況                    */
/*              >= 0 : 正常終了(収束回数) */
/*              = -1 : 収束せず           */
/*      max : 最大試行回数                */
/*      fun : 関数値を計算する関数の名前  */
/*      return : 結果                     */
/******************************************/
#include <math.h>

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

	tau   = (sqrt(5.0) - 1.0) / 2.0;
	count = 0;
	*ind  = -1;
	x1    = b - tau * (b - a);
	x2    = a + tau * (b - a);
	f1    = fun(x1);
	f2    = fun(x2);

	while (count < max && *ind < 0) {
		count += 1;
		if (f2 > f1) {
			if (fabs(b-a) < eps && fabs(b-a) < eps*fabs(b)) {
				*ind = 0;
				x    = x1;
				*val = f1;
			}
			else {
				b  = x2;
				x2 = x1;
				x1 = a + (1.0 - tau) * (b - a);
				f2 = f1;
				f1 = fun(x1);
			}
		}
		else {
			if (fabs(b-a) < eps && fabs(b-a) < eps*fabs(b)) {
				*ind = 0;
				x    = x2;
				*val = f2;
				f1   = f2;
			}
			else {
				a  = x1;
				x1 = x2;
				x2 = b - (1.0 - tau) * (b - a);
				f1 = f2;
				f2 = fun(x2);
			}
		}
	}

	if (*ind == 0) {
		*ind = count;
		fa   = fun(a);
		fb   = fun(b);
		if (fb < fa) {
			a  = b;
			fa = fb;
		}
		if (fa < f1) {
			x    = a;
			*val = fa;
		}
	}

	return x;
}

  2. Java 

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

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

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

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

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

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

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

	abstract double snx(double x);

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

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

		while (count < max && ind[0] < 0) {
			count += 1;
			if (f2 > f1) {
				if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
					ind[0] = 0;
					x      = x1;
					val[0] = f1;
				}
				else {
					b  = x2;
					x2 = x1;
					x1 = a + (1.0 - tau) * (b - a);
					f2 = f1;
					f1 = snx(x1);
				}
			}
			else {
				if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
					ind[0] = 0;
					x      = x2;
					val[0] = f2;
					f1     = f2;
				}
				else {
					a  = x1;
					x1 = x2;
					x2 = b - (1.0 - tau) * (b - a);
					f1 = f2;
					f2 = snx(x2);
				}
			}
		}

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

		return x;
	}
}

  3. JavaScript

    　　ここをクリックすると，JavaScript の仕様に適合した形で最小値を求めたい式を入力することによって，任意の関数の最小値を画面上で求めることができます． 
    <!DOCTYPE HTML>

<HTML>

<HEAD>

	<TITLE>最適化（黄金分割法）</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		str = "";
		function main()
		{
			str = document.getElementById("func").value + ";";
					// データの設定
			let a   = parseFloat(document.getElementById("a").value);
			let b   = parseFloat(document.getElementById("b").value);
			let max = parseInt(document.getElementById("max").value);
			let eps = 1.0e-10;
			let val = new Array();
			let ind = new Array();
					// 実行
			let x = gold(a, b, eps, val, ind, max, snx);
					// 結果
			if (ind < 0)
				document.getElementById("res").value = "解を求めることができませんでした！";
			else {
				let c = 100000;
				let s1 = Math.round(x * c) / c;
				let s2 = Math.round(val[0] * c) / c;
				document.getElementById("res").value = " x = " + s1 + "，f(x) = " + s2 + "，ind = " + ind[0];
			}
		}

		/****************/
		/* 関数値の計算 */
		/****************/
		function snx(x)
		{
			let y = eval(str);
			return y;
		}

		/*******************************************/
		/* 黄金分割法(関数の最小値)                */
		/*      a,b : 初期区間 a < b               */
		/*      eps : 許容誤差                     */
		/*      val : 間数値                       */
		/*      ind : 計算状況                     */
		/*              >= 0 : 正常終了(収束回数)  */
		/*              = -1 : 収束せず            */
		/*      max : 最大試行回数                 */
		/*      return : 結果                      */
		/*******************************************/
		function gold(a, b, eps, val, ind, max)
		{
			let tau   = (Math.sqrt(5.0) - 1.0) / 2.0;
			let count = 0;
			ind[0]    = -1;
			let x     = 0.0;
			let x1    = b - tau * (b - a);
			let x2    = a + tau * (b - a);
			let f1    = snx(x1);
			let f2    = snx(x2);

			while (count < max && ind[0] < 0) {
				count += 1;
				if (f2 > f1) {
					if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
						ind[0] = 0;
						x      = x1;
						val[0] = f1;
					}
					else {
						b  = x2;
						x2 = x1;
						x1 = a + (1.0 - tau) * (b - a);
						f2 = f1;
						f1 = snx(x1);
					}
				}
				else {
					if (Math.abs(b-a) < eps && Math.abs(b-a) < eps*Math.abs(b)) {
						ind[0] = 0;
						x      = x2;
						val[0] = f2;
						f1     = f2;
					}
					else {
						a  = x1;
						x1 = x2;
						x2 = b - (1.0 - tau) * (b - a);
						f1 = f2;
						f2 = snx(x2);
					}
				}
			}

			if (ind[0] == 0) {
				ind[0] = count;
				let fa = snx(a);
				let fb = snx(b);
				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>　　テキストフィールドには，例として，以下に示すような関数の [-2, -1] 区間における最小値を求める場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．なお，目的関数は，JavaScript の仕様に適合した形式で記述してあることに注意してください． 
		<P STYLE="text-align:center"><IMG SRC="gold.gif"></P>
	</DL>

	<DIV STYLE="text-align:center">
		初期値：<INPUT ID="a" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="-2">，
				<INPUT ID="b" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="-1">　
		最大繰り返し回数：<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="100"><BR><BR>
		目的関数：　f(x) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="50" VALUE="x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0">　　
		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
		結果：<INPUT ID="res" STYLE="font-size: 100%" TYPE="text" SIZE="40">
	</DIV>

</BODY>

</HTML>

  4. PHP 

    <?php

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

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

	$x = gold($a, $b, $eps, $val, $ind, $max, "snx");

	printf("x %f val %f ind %d\n", $x, $val, $ind);

/****************/
/* 関数値の計算 */
/****************/
function snx($x)
{
	return $x * $x * $x * $x + 3.0 * $x * $x * $x + 2.0 * $x * $x + 1.0;
}

/******************************************/
/* 黄金分割法(関数の最小値)               */
/*      a,b : 初期区間 a < b              */
/*      eps : 許容誤差                    */
/*      val : 間数値                      */
/*      ind : 計算状況                    */
/*              >= 0 : 正常終了(収束回数) */
/*              = -1 : 収束せず           */
/*      max : 最大試行回数                */
/*      fun : 関数値を計算する関数の名前  */
/*      return : 結果                     */
/******************************************/
function gold($a, $b, $eps, &$val, &$ind, $max, $fun)
{
	$x     = 0.0;
	$tau   = (sqrt(5.0) - 1.0) / 2.0;
	$count = 0;
	$ind   = -1;
	$x1    = $b - $tau * ($b - $a);
	$x2    = $a + $tau * ($b - $a);
	$f1    = $fun($x1);
	$f2    = $fun($x2);

	while ($count < $max && $ind < 0) {
		$count += 1;
		if ($f2 > $f1) {
			if (abs($b-$a) < $eps && abs($b-$a) < $eps*abs($b)) {
				$ind = 0;
				$x   = $x1;
				$val = $f1;
			}
			else {
				$b  = $x2;
				$x2 = $x1;
				$x1 = $a + (1.0 - $tau) * ($b - $a);
				$f2 = $f1;
				$f1 = $fun($x1);
			}
		}
		else {
			if (abs($b-$a) < $eps && abs($b-$a) < $eps*abs($b)) {
				$ind = 0;
				$x   = $x2;
				$val = $f2;
				$f1  = $f2;
			}
			else {
				$a  = $x1;
				$x1 = $x2;
				$x2 = $b - (1.0 - $tau) * ($b - $a);
				$f1 = $f2;
				$f2 = $fun($x2);
			}
		}
	}

	if ($ind == 0) {
		$ind = $count;
		$fa  = $fun($a);
		$fb  = $fun($b);
		if ($fb < $fa) {
			$a  = $b;
			$fa = $fb;
		}
		if ($fa < $f1) {
			$x   = $a;
			$val = $fa;
		}
	}

	return $x;
}

?>

  5. Ruby 

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

#***************/
# 関数値の計算 */
#***************/
snx = Proc.new { |x|
	x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0
}

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

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

	while count < max && ind[0] < 0
		count += 1
		if f2 > f1
			if (b-a).abs() < eps && (b-a).abs() < eps*b.abs()
				ind[0] = 0
				x      = x1
				val[0] = f1
			else
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = fun.call(x1)
			end
		else
			if (b-a).abs() < eps && (b-a).abs() < eps*b.abs()
				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 = fun.call(x2)
			end
		end
	end

	if ind[0] == 0
		ind[0] = count
		fa     = fun.call(a)
		fb     = fun.call(b)
		if fb < fa
			a  = b
			fa = fb
		end
		if fa < f1
			x      = a
			val[0] = fa
		end
	end

	return x
end

ind = Array.new(1)
val = Array.new(1)
a   = -2.0
b   = -1.0
eps = 1.0e-10
max = 100

x = gold(a, b, eps, val, ind, max, &snx)

printf("x %f val %f ind %d\n", x, val[0], ind[0])

  6. Python 

    # -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *

############################################
# 黄金分割法(関数の最小値)
#      a,b : 初期区間 a < b
#      eps : 許容誤差
#      val : 間数値
#      ind : 計算状況
#              >= 0 : 正常終了(収束回数)
#              = -1 : 収束せず
#      max : 最大試行回数
#      fun : 関数値を計算する関数の名前
#      return : 結果
#      coded by Y.Suganuma
############################################

def gold(a, b, eps, val, ind, max, fun) :

	x      = 0.0
	tau    = (sqrt(5.0) - 1.0) / 2.0
	count  = 0
	ind[0] = -1
	x1     = b - tau * (b - a)
	x2     = a + tau * (b - a)
	f1     = fun(x1)
	f2     = fun(x2)

	while count < max and ind[0] < 0 :
		count += 1
		if f2 > f1 :
			if (abs(b-a) < eps) and (abs(b-a) < eps*abs(b)) :
				ind[0] = 0
				x      = x1
				val[0] = f1
			else :
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = fun(x1)
		else :
			if (abs(b-a) < eps) and (abs(b-a) < eps*abs(b)) :
				ind[0] = 0
				x      = x2
				val[0] = f2
				f1     = f2
			else :
				a  = x1
				x1 = x2
				x2 = b - (1.0 - tau) * (b - a)
				f1 = f2
				f2 = fun(x2)

	if ind[0] == 0 :
		ind[0] = count
		fa     = fun(a)
		fb     = fun(b)
		if fb < fa :
			a  = b
			fa = fb
		if fa < f1 :
			x      = a
			val[0] = fa

	return x

############################################
# 黄金分割法によるy=x^4+3x^3+2x^2+1の最小値
#      coded by Y.Suganuma
############################################

			# 関数値の計算
def snx(x) :
	return x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0
			# 設定と実行
a   = -2.0
b   = -1.0
eps = 1.0e-10
max = 100
val = np.empty(1, np.float)
ind = np.empty(1, np.int)

x = gold(a, b, eps, val, ind, max, snx)

print("x " + str(x) + " val " + str(val[0]) + " ind " + str(ind[0]))

  7. C# 

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

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	public Test1()
	{

		double a   = -2.0;
		double b   = -1.0;
		double eps = 1.0e-10;
		double val = 0.0;
		int max    = 100;
		int ind    = 0;

		double x = gold(a, b, eps, ref val, ref ind, max, snx);

		Console.WriteLine("x " + x + " val " + val + " ind " + ind);
	}
			// 関数値の計算
	double snx(double x)
	{
		return x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0;
	}

	/************************************************/
	/* 黄金分割法(関数の最小値)                     */
	/*      a,b : 初期区間 a < b                    */
	/*      eps : 許容誤差                          */
	/*      val : 間数値                            */
	/*      ind : 計算状況                          */
	/*              >= 0 : 正常終了(収束回数)       */
	/*              = -1 : 収束せず                 */
	/*      max : 最大試行回数                      */
	/*      fn : 関数値を計算する関数               */
	/*      return : 結果                           */
	/************************************************/
	double gold(double a, double b, double eps, ref double val, ref int ind, int max,
	            Func<double, double> fn)
	{
		ind        = -1;
		int count  = 0;
		double tau = (Math.Sqrt(5.0) - 1.0) / 2.0;
		double x1  = b - tau * (b - a);
		double x2  = a + tau * (b - a);
		double f1  = fn(x1);
		double f2  = fn(x2);
		double x   = 0.0;

		while (count < max && ind < 0) {
			count += 1;
			if (f2 > f1) {
				if (Math.Abs(b-a) < eps && Math.Abs(b-a) < eps*Math.Abs(b)) {
					ind = 0;
					x   = x1;
					val = f1;
				}
				else {
					b  = x2;
					x2 = x1;
					x1 = a + (1.0 - tau) * (b - a);
					f2 = f1;
					f1 = fn(x1);
				}
			}
			else {
				if (Math.Abs(b-a) < eps && Math.Abs(b-a) < eps*Math.Abs(b)) {
					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);
				}
			}
		}

		if (ind == 0) {
			ind       = count;
			double fa = fn(a);
			double fb = fn(b);
			if (fb < fa) {
				a  = b;
				fa = fb;
			}
			if (fa < f1) {
				x   = a;
				val = fa;
			}
		}

		return x;
	}
}

  8. VB 

    '''''''''''''''''''''''''''''''''''''''''''''
' 黄金分割法によるy=x^4+3x^3+2x^2+1の最小値 '
'        coded by Y.Suganuma                '
'''''''''''''''''''''''''''''''''''''''''''''
Module Test

	Sub Main()

		Dim a As Double    = -2.0
		Dim b As Double    = -1.0
		Dim eps As Double  = 1.0e-10
		Dim val As Double  = 0.0
		Dim max As Integer = 100
		Dim ind As Integer = 0
			' 関数値の計算（ラムダ式）
		Dim snx = Function(v) As Double 
			Return v * v * v * v + 3.0 * v * v * v + 2.0 * v * v + 1.0
		End Function
			' 実行
		Dim x As Double = gold(a, b, eps, val, ind, max, snx)

		Console.WriteLine("x " & x & " val " & val & " ind " & ind)

	End Sub

	''''''''''''''''''''''''''''''''''''''''''''''''
	' 黄金分割法(関数の最小値)                     '
	'      a,b : 初期区間 a < b                    '
	'      eps : 許容誤差                          '
	'      val : 間数値                            '
	'      ind : 計算状況                          '
	'              >= 0 : 正常終了(収束回数)       '
	'              = -1 : 収束せず                 '
	'      max : 最大試行回数                      '
	'      fn : 関数値を計算する関数               '
	'      return : 結果                           '
	''''''''''''''''''''''''''''''''''''''''''''''''
	Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
	              ByRef ind As Integer, max As Integer, fn As Func(Of Double, Double))

		ind = -1
		Dim count As Integer  = 0
		Dim tau As Double = (Math.Sqrt(5.0) - 1.0) / 2.0
		Dim x1 As Double  = b - tau * (b - a)
		Dim x2 As Double  = a + tau * (b - a)
		Dim f1 As Double  = fn(x1)
		Dim f2 As Double  = fn(x2)
		Dim x As Double   = 0.0

		Do while count < max and ind < 0
			count += 1
			If f2 > f1
				If Math.Abs(b-a) < eps and Math.Abs(b-a) < eps*Math.Abs(b)
					ind = 0
					x   = x1
					val = f1
				Else
					b  = x2
					x2 = x1
					x1 = a + (1.0 - tau) * (b - a)
					f2 = f1
					f1 = fn(x1)
				End If
			Else
				If Math.Abs(b-a) < eps and Math.Abs(b-a) < eps*Math.Abs(b)
					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)
				End If
			End If
		Loop

		If ind = 0
			ind = count
			Dim fa As Double = fn(a)
			Dim fb As Double = fn(b)
			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

情報学部 菅沼ホーム 目次 索引