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

最適化(多項式近似法)

    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 approx(double, double, double, double *, int *, int, double (*)(double));
    
    int main()
    {
    	double eps, val, x, d;
    	int ind, max;
    
    	x   = -3.0;
    	d   = 0.1;
    	eps = 1.0e-7;
    	max = 100;
    
    	x = approx(x, d, 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;
    }
    
    /******************************************/
    /* 多項式近似(関数の最小値)               */
    /*      x0  : 初期値                      */
    /*      d0  : 初期ステップ                */
    /*      eps : 許容誤差                    */
    /*      val : 関数値                      */
    /*      ind : 計算状況                    */
    /*              >= 0 : 正常終了(収束回数) */
    /*              = -1 : 収束せず           */
    /*      max : 最大試行回数                */
    /*      fun : 関数値を計算する関数の名前  */
    /*      return : 結果                     */
    /******************************************/
    #include <math.h>
    
    double approx(double x0, double d0, double eps, double *val, int *ind, int max, double (*fun)(double))
    {
    	double f[4], x[4], xx = 0.0, d, dl;
    	int i1, k = 0, count = 0, min, sw;
    
    	d    = d0;
    	x[1] = x0;
    	f[1] = fun(x0);
    	*ind = -1;
    
    	while (count < max && *ind < 0) {
    		x[3] = x[1] + d;
    		f[3] = fun(x[3]);
    		while (k < max && f[3] <= f[1]) {
    			k++;
    			d *= 2.0;
    			x[0] = x[1];
    			f[0] = f[1];
    			x[1] = x[3];
    			f[1] = f[3];
    			x[3] = x[1] + d;
    			f[3] = fun(x[3]);
    		}
    					// 初期値が不適当
    		if (k >= max)
    			count = max;
    		else {
    					// 3点の選択
    			sw = 0;
    			if (k > 0) {
    				x[2] = x[3] - 0.5 * d;
    				f[2] = fun(x[2]);
    				min  = -1;
    				for (i1 = 0; i1 < 4; i1++) {
    					if (min < 0 || f[i1] < f[min])
    						min = i1;
    				}
    				if (min >= 2) {
    					for (i1 = 0; i1 < 3; i1++) {
    						x[i1] = x[i1+1];
    						f[i1] = f[i1+1];
    					}
    				}
    				sw = 1;
    			}
    			else {
    				x[0] = x[1] - d0;
    				f[0] = fun(x[0]);
    				if (f[0] > f[1]) {
    					x[2] = x[3];
    					f[2] = f[3];
    					sw = 1;
    				}
    				else {
    					x[1] = x[0];
    					f[1] = f[0];
    					d0   = -d0;
    					d    = 2.0 * d0;
    					k    = 1;
    				}
    			}
    					// 収束?
    			if (sw > 0) {
    				count++;
    				dl = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2]);
    				xx   = x[1] - dl;
    				*val = fun(xx);
    				if (fabs(dl) < eps)
    					*ind = count;
    				else {
    					k  = 0;
    					d0 = 0.5 * d;
    					d  = d0;
    					if (*val < f[1]) {
    						x[1] = xx;
    						f[1] = *val;
    					}
    				}
    			}
    		}
    	}
    
    	return xx;
    }
    			

  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 eps, x, d, val[] = new double [1];
    		int max, ind[] = new int [1];
    
    		x   = -2.0;
    		d   = 0.1;
    		eps = 1.0e-10;
    		max = 100;
    
    		Kansu kn = new Kansu();
    		x = kn.approx(x, d, eps, val, ind, max);
    
    		System.out.println("x " + x + " val " + val[0] + " ind " + ind[0]);
    	}
    }
    
    /****************/
    /* 関数値の計算 */
    /****************/
    class Kansu extends Approximation
    {
    	double snx(double x)
    	{
    		double y = x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0;
    		return y;
    	}
    }
    
    abstract class Approximation {
    
    	/************************************************/
    	/* 多項式近似(関数の最小値)                     */
    	/*      x0  : 初期値                            */
    	/*      d0  : 初期ステップ                      */
    	/*      eps : 許容誤差                          */
    	/*      val : 間数値                            */
    	/*      ind : 計算状況                          */
    	/*              >= 0 : 正常終了(収束回数)       */
    	/*              = -1 : 収束せず                 */
    	/*      max : 最大試行回数                      */
    	/*      return : 結果                           */
    	/************************************************/
    
    	abstract double snx(double x);
    
    	double approx(double x0, double d0, double eps, double val[], int ind[], int max)
    	{
    		double f[] = new double [4];
    		double x[] = new double [4];
    		double xx = 0.0, d, dl;
    		int i1, k = 0, count = 0, min, sw = 0;
    
    		d      = d0;
    		x[1]   = x0;
    		f[1]   = snx(x0);
    		ind[0] = -1;
    
    		while (count < max && ind[0] < 0) {
    			x[3] = x[1] + d;
    			f[3] = snx(x[3]);
    			while (k < max && f[3] <= f[1]) {
    				k++;
    				d *= 2.0;
    				x[0] = x[1];
    				f[0] = f[1];
    				x[1] = x[3];
    				f[1] = f[3];
    				x[3] = x[1] + d;
    				f[3] = snx(x[3]);
    			}
    						// 初期値が不適当
    			if (k >= max)
    				count = max;
    			else {
    						// 3点の選択
    				sw = 0;
    				if (k > 0) {
    					x[2] = x[3] - 0.5 * d;
    					f[2] = snx(x[2]);
    					min  = -1;
    					for (i1 = 0; i1 < 4; i1++) {
    						if (min < 0 || f[i1] < f[min])
    							min = i1;
    					}
    					if (min >= 2) {
    						for (i1 = 0; i1 < 3; i1++) {
    							x[i1] = x[i1+1];
    							f[i1] = f[i1+1];
    						}
    					}
    					sw = 1;
    				}
    				else {
    					x[0] = x[1] - d0;
    					f[0] = snx(x[0]);
    					if (f[0] > f[1]) {
    						x[2] = x[3];
    						f[2] = f[3];
    						sw = 1;
    					}
    					else {
    						x[1] = x[0];
    						f[1] = f[0];
    						d0   = -d0;
    						d    = 2.0 * d0;
    						k    = 1;
    					}
    				}
    						// 収束?
    				if (sw > 0) {
    					count++;
    					dl = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2]);
    					xx     = x[1] - dl;
    					val[0] = snx(xx);
    					if (Math.abs(dl) < eps)
    						ind[0] = count;
    					else {
    						k  = 0;
    						d0 = 0.5 * d;
    						d  = d0;
    						if (val[0] < f[1]) {
    							x[1] = xx;
    							f[1] = val[0];
    						}
    					}
    				}
    			}
    		}
    
    		return xx;
    	}
    }
    			

  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 x   = parseFloat(document.getElementById("x0").value);
    			let d   = parseFloat(document.getElementById("d0").value);
    			let max = parseInt(document.getElementById("max").value);
    			let eps = 1.0e-10;
    			let val = new Array();
    			let ind = new Array();
    					// 実行
    			x = approx(x, d, 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;
    		}
    
    		/*******************************************/
    		/* 多項式近似(関数の最小値)                */
    		/*      x0  : 初期値                       */
    		/*      d0  : 初期ステップ                 */
    		/*      eps : 許容誤差                     */
    		/*      val : 間数値                       */
    		/*      ind : 計算状況                     */
    		/*              >= 0 : 正常終了(収束回数)  */
    		/*              = -1 : 収束せず            */
    		/*      max : 最大試行回数                 */
    		/*      fn : 関数値を計算する関数          */
    		/*      return : 結果                      */
    		/*******************************************/
    		function approx(x0, d0, eps, val, ind, max, fn)
    		{
    			let f = new Array();
    			let x = new Array();
    			let xx = 0.0;
    			let dl;
    			let k = 0;
    			let count = 0;
    			let min;
    			let sw = 0;
    
    			let d  = d0;
    			x[1]   = x0;
    			f[1]   = fn(x0);
    			ind[0] = -1;
    
    			while (count < max && ind[0] < 0) {
    				x[3] = x[1] + d;
    				f[3] = fn(x[3]);
    				while (k < max && f[3] <= f[1]) {
    					k++;
    					d *= 2.0;
    					x[0] = x[1];
    					f[0] = f[1];
    					x[1] = x[3];
    					f[1] = f[3];
    					x[3] = x[1] + d;
    					f[3] = fn(x[3]);
    				}
    						// 初期値が不適当
    				if (k >= max)
    					count = max;
    				else {
    						// 3点の選択
    					sw = 0;
    					if (k > 0) {
    						x[2] = x[3] - 0.5 * d;
    						f[2] = fn(x[2]);
    						min  = -1;
    						for (i1 = 0; i1 < 4; i1++) {
    							if (min < 0 || f[i1] < f[min])
    								min = i1;
    						}
    						if (min >= 2) {
    							for (i1 = 0; i1 < 3; i1++) {
    								x[i1] = x[i1+1];
    								f[i1] = f[i1+1];
    							}
    						}
    						sw = 1;
    					}
    					else {
    						x[0] = x[1] - d0;
    						f[0] = fn(x[0]);
    						if (f[0] > f[1]) {
    							x[2] = x[3];
    							f[2] = f[3];
    							sw = 1;
    						}
    						else {
    							x[1] = x[0];
    							f[1] = f[0];
    							d0   = -d0;
    							d    = 2.0 * d0;
    							k    = 1;
    						}
    					}
    						// 収束?
    					if (sw > 0) {
    						count++;
    						dl = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2]);
    						xx     = x[1] - dl;
    						val[0] = fn(xx);
    						if (Math.abs(dl) < eps)
    							ind[0] = count;
    						else {
    							k  = 0;
    							d0 = 0.5 * d;
    							d  = d0;
    							if (val[0] < f[1]) {
    								x[1] = xx;
    								f[1] = val[0];
    							}
    						}
    					}
    				}
    			}
    
    			return xx;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>最適化(多項式近似法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドには,例として,以下に示すような関数の -2 の近傍における最小値を求める場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.なお,目的関数は,JavaScript の仕様に適合した形式で記述してあることに注意してください.
    		<P STYLE="text-align:center"><IMG SRC="approx.gif"></P>
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		初期値:<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="-2"> 
    		初期ステップ幅:<INPUT ID="d0" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="0.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="50">
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /*********************************************/
    /* 多項式近似によるy=x^4+3x^3+2x^2+1の最小値 */
    /*      coded by Y.Suganuma                  */
    /*********************************************/
    
    	$x   = -3.0;
    	$d   = 0.1;
    	$eps = 1.0e-7;
    	$max = 100;
    
    	$x = approx($x, $d, $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;
    }
    
    /******************************************/
    /* 多項式近似(関数の最小値)               */
    /*      x0  : 初期値                      */
    /*      d0  : 初期ステップ                */
    /*      eps : 許容誤差                    */
    /*      val : 関数値                      */
    /*      ind : 計算状況                    */
    /*              >= 0 : 正常終了(収束回数) */
    /*              = -1 : 収束せず           */
    /*      max : 最大試行回数                */
    /*      fun : 関数値を計算する関数の名前  */
    /*      return : 結果                     */
    /******************************************/
    function approx($x0, $d0, $eps, &$val, &$ind, $max, $fun)
    {
    	$xx    = 0.0;
    	$k     = 0;
    	$count = 0;
    	$d     = $d0;
    	$x     = array(4);
    	$x[1]  = $x0;
    	$f     = array(4);
    	$f[1]  = $fun($x0);
    	$ind   = -1;
    
    	while ($count < $max && $ind < 0) {
    		$x[3] = $x[1] + $d;
    		$f[3] = $fun($x[3]);
    		while ($k < $max && $f[3] <= $f[1]) {
    			$k++;
    			$d *= 2.0;
    			$x[0] = $x[1];
    			$f[0] = $f[1];
    			$x[1] = $x[3];
    			$f[1] = $f[3];
    			$x[3] = $x[1] + $d;
    			$f[3] = $fun($x[3]);
    		}
    					// 初期値が不適当
    		if ($k >= $max)
    			$count = $max;
    		else {
    					// 3点の選択
    			$sw = 0;
    			if ($k > 0) {
    				$x[2] = $x[3] - 0.5 * $d;
    				$f[2] = $fun($x[2]);
    				$min  = -1;
    				for ($i1 = 0; $i1 < 4; $i1++) {
    					if ($min < 0 || $f[$i1] < $f[$min])
    						$min = $i1;
    				}
    				if ($min >= 2) {
    					for ($i1 = 0; $i1 < 3; $i1++) {
    						$x[$i1] = $x[$i1+1];
    						$f[$i1] = $f[$i1+1];
    					}
    				}
    				$sw = 1;
    			}
    			else {
    				$x[0] = $x[1] - $d0;
    				$f[0] = $fun($x[0]);
    				if ($f[0] > $f[1]) {
    					$x[2] = $x[3];
    					$f[2] = $f[3];
    					$sw   = 1;
    				}
    				else {
    					$x[1] = $x[0];
    					$f[1] = $f[0];
    					$d0   = -$d0;
    					$d    = 2.0 * $d0;
    					$k    = 1;
    				}
    			}
    					// 収束?
    			if ($sw > 0) {
    				$count++;
    				$dl  = 0.5 * $d * ($f[2] - $f[0]) / ($f[0] - 2.0 * $f[1] + $f[2]);
    				$xx  = $x[1] - $dl;
    				$val = $fun($xx);
    				if (abs($dl) < $eps)
    					$ind = $count;
    				else {
    					$k  = 0;
    					$d0 = 0.5 * $d;
    					$d  = $d0;
    					if ($val < $f[1]) {
    						$x[1] = $xx;
    						$f[1] = $val;
    					}
    				}
    			}
    		}
    	}
    
    	return $xx;
    }
    
    ?>
    			

  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
    }
    
    #*****************************************/
    # 多項式近似(関数の最小値)               */
    #      x0  : 初期値                      */
    #      d0  : 初期ステップ                */
    #      eps : 許容誤差                    */
    #      val : 関数値                      */
    #      ind : 計算状況                    */
    #              >= 0 : 正常終了(収束回数) */
    #              = -1 : 収束せず           */
    #      max : 最大試行回数                */
    #      fun : 関数値を計算する関数の名前  */
    #      return : 結果                     */
    #*****************************************/
    def approx(x0, d0, eps, val, ind, max, &fun)
    
    	f      = Array.new(4)
    	x      = Array.new(4)
    	xx     = 0.0
    	d      = d0
    	x[1]   = x0
    	f[1]   = fun.call(x0)
    	k      = 0
    	count  = 0
    	ind[0] = -1
    
    	while count < max && ind[0] < 0
    		x[3] = x[1] + d
    		f[3] = fun.call(x[3])
    		while k < max && f[3] <= f[1]
    			k   += 1
    			d   *= 2.0
    			x[0] = x[1]
    			f[0] = f[1]
    			x[1] = x[3]
    			f[1] = f[3]
    			x[3] = x[1] + d
    			f[3] = fun.call(x[3])
    		end
    					# 初期値が不適当
    		if k >= max
    			count = max
    		else
    					# 3点の選択
    			sw = 0
    			if k > 0
    				x[2] = x[3] - 0.5 * d
    				f[2] = fun.call(x[2])
    				min  = -1
    				for i1 in 0 ... 4
    					if min < 0 || f[i1] < f[min]
    						min = i1
    					end
    				end
    				if min >= 2
    					for i1 in 0 ... 3
    						x[i1] = x[i1+1]
    						f[i1] = f[i1+1]
    					end
    				end
    				sw = 1
    			else
    				x[0] = x[1] - d0
    				f[0] = fun.call(x[0])
    				if f[0] > f[1]
    					x[2] = x[3]
    					f[2] = f[3]
    					sw   = 1
    				else
    					x[1] = x[0]
    					f[1] = f[0]
    					d0   = -d0
    					d    = 2.0 * d0
    					k    = 1
    				end
    			end
    					# 収束?
    			if sw > 0
    				count += 1
    				dl     = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2])
    				xx     = x[1] - dl
    				val[0] = fun.call(xx)
    				if dl.abs() < eps
    					ind[0] = count
    				else
    					k  = 0
    					d0 = 0.5 * d
    					d  = d0
    					if val[0] < f[1]
    						x[1] = xx
    						f[1] = val[0]
    					end
    				end
    			end
    		end
    	end
    
    	return xx
    end
    
    ind = Array.new(1)
    val = Array.new(1)
    x   = -3.0
    d   = 0.1
    eps = 1.0e-7
    max = 100
    
    x = approx(x, d, 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 *
    
    ##########################################
    # 多項式近似(関数の最小値)
    #      x0  : 初期値
    #      d0  : 初期ステップ
    #      eps : 許容誤差
    #      val : 間数値
    #      ind : 計算状況
    #              >= 0 : 正常終了(収束回数)
    #              = -1 : 収束せず
    #      max : 最大試行回数
    #      fun : 関数値を計算する関数の名前
    #      return : 結果
    #      coded by Y.Suganuma
    ##########################################
    
    def approx(x0, d0, eps, val, ind, max, fun) :
    
    	x      = np.empty(4, np.float)
    	f      = np.empty(4, np.float)
    	xx     = 0.0
    	k      = 0
    	count  = 0
    	d      = d0
    	x[1]   = x0
    	f[1]   = fun(x0)
    	ind[0] = -1
    
    	while count < max and ind[0] < 0 :
    		x[3] = x[1] + d
    		f[3] = fun(x[3])
    		while k < max and f[3] <= f[1] :
    			k += 1
    			d *= 2.0
    			x[0] = x[1]
    			f[0] = f[1]
    			x[1] = x[3]
    			f[1] = f[3]
    			x[3] = x[1] + d
    			f[3] = fun(x[3])
    					# 初期値が不適当
    		if k >= max :
    			count = max
    		else :
    					# 3点の選択
    			sw = 0
    			if k > 0 :
    				x[2] = x[3] - 0.5 * d
    				f[2] = fun(x[2])
    				min  = -1
    				for i1 in range(0, 4) :
    					if min < 0 or f[i1] < f[min] :
    						min = i1
    				if min >= 2 :
    					for i1 in range(0, 3) :
    						x[i1] = x[i1+1]
    						f[i1] = f[i1+1]
    				sw = 1
    			else :
    				x[0] = x[1] - d0
    				f[0] = fun(x[0])
    				if f[0] > f[1] :
    					x[2] = x[3]
    					f[2] = f[3]
    					sw = 1
    				else :
    					x[1] = x[0]
    					f[1] = f[0]
    					d0   = -d0
    					d    = 2.0 * d0
    					k    = 1
    					# 収束?
    			if sw > 0 :
    				count += 1
    				dl     = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2])
    				xx     = x[1] - dl
    				val[0] = fun(xx)
    				if abs(dl) < eps :
    					ind[0] = count
    				else :
    					k  = 0
    					d0 = 0.5 * d
    					d  = d0
    					if val[0] < f[1] :
    						x[1] = xx
    						f[1] = val[0]
    
    	return xx
    
    ############################################
    # 多項式近似による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
    
    			# 設定と実行
    x   = -3.0
    d   = 0.1
    eps = 1.0e-7
    max = 100
    val = np.empty(1, np.float)
    ind = np.empty(1, np.int)
    
    x = approx(x, d, 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 val = 0.0;
    		double x   = -2.0;
    		double d   = 0.1;
    		double eps = 1.0e-10;
    		int max    = 100;
    		int ind    = 0;
    
    		x = approx(x, d, 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;
    	}
    
    	/************************************************/
    	/* 多項式近似(関数の最小値)                     */
    	/*      x0  : 初期値                            */
    	/*      d0  : 初期ステップ                      */
    	/*      eps : 許容誤差                          */
    	/*      val : 間数値                            */
    	/*      ind : 計算状況                          */
    	/*              >= 0 : 正常終了(収束回数)       */
    	/*              = -1 : 収束せず                 */
    	/*      max : 最大試行回数                      */
    	/*      fn : 関数値を計算する関数               */
    	/*      return : 結果                           */
    	/************************************************/
    	double approx(double x0, double d0, double eps, ref double val, ref int ind, int max,
    	              Func<double, double> fn)
    	{
    		double[] f = new double [4];
    		f[1]       = fn(x0);
    		double[] x = new double [4];
    		x[1]       = x0;
    		double xx  = 0.0;
    		double d   = d0;
    		int count  = 0;
    		ind        = -1;
    
    		while (count < max && ind < 0) {
    			x[3]  = x[1] + d;
    			f[3]  = fn(x[3]);
    			int k = 0;
    			while (k < max && f[3] <= f[1]) {
    				k++;
    				d   *= 2.0;
    				x[0] = x[1];
    				f[0] = f[1];
    				x[1] = x[3];
    				f[1] = f[3];
    				x[3] = x[1] + d;
    				f[3] = fn(x[3]);
    			}
    						// 初期値が不適当
    			if (k >= max)
    				count = max;
    			else {
    						// 3点の選択
    				int sw = 0;
    				if (k > 0) {
    					x[2] = x[3] - 0.5 * d;
    					f[2] = fn(x[2]);
    					int min  = -1;
    					for (int i1 = 0; i1 < 4; i1++) {
    						if (min < 0 || f[i1] < f[min])
    							min = i1;
    					}
    					if (min >= 2) {
    						for (int i1 = 0; i1 < 3; i1++) {
    							x[i1] = x[i1+1];
    							f[i1] = f[i1+1];
    						}
    					}
    					sw = 1;
    				}
    				else {
    					x[0] = x[1] - d0;
    					f[0] = fn(x[0]);
    					if (f[0] > f[1]) {
    						x[2] = x[3];
    						f[2] = f[3];
    						sw   = 1;
    					}
    					else {
    						x[1] = x[0];
    						f[1] = f[0];
    						d0   = -d0;
    						d    = 2.0 * d0;
    						k    = 1;
    					}
    				}
    						// 収束?
    				if (sw > 0) {
    					count++;
    					double dl = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2]);
    					xx  = x[1] - dl;
    					val = fn(xx);
    					if (Math.Abs(dl) < eps)
    						ind = count;
    					else {
    						k  = 0;
    						d0 = 0.5 * d;
    						d  = d0;
    						if (val < f[1]) {
    							x[1] = xx;
    							f[1] = val;
    						}
    					}
    				}
    			}
    		}
    
    		return xx;
    	}
    }
    			

  8. VB

    '''''''''''''''''''''''''''''''''''''''''''''
    ' 多項式近似によるy=x^4+3x^3+2x^2+1の最小値 '
    '        coded by Y.Suganuma                '
    '''''''''''''''''''''''''''''''''''''''''''''
    Module Test
    
    	Sub Main()
    
    		Dim val As Double  = 0.0
    		Dim x As Double    = -2.0
    		Dim d As Double    = 0.1
    		Dim eps As Double  = 1.0e-10
    		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
    			' 実行
    		x = approx(x, d, eps, val, ind, max, snx)
    
    		Console.WriteLine("x " & x & " val " & val & " ind " & ind)
    
    	End Sub
    
    	''''''''''''''''''''''''''''''''''''''''''''''''
    	' 多項式近似(関数の最小値)                     '
    	'      x0  : 初期値                            '
    	'      d0  : 初期ステップ                      '
    	'      eps : 許容誤差                          '
    	'      val : 間数値                            '
    	'      ind : 計算状況                          '
    	'              >= 0 : 正常終了(収束回数)       '
    	'              = -1 : 収束せず                 '
    	'      max : 最大試行回数                      '
    	'      fn : 関数値を計算する関数               '
    	'      return : 結果                           '
    	''''''''''''''''''''''''''''''''''''''''''''''''
    	Function approx(x0 As Double, d0 As Double, eps As Double, Byref val As Double,
    	                ByRef ind As Integer, max As Integer, fn As Func(Of Double, Double))
    
    		Dim f(4) As Double
    		f(1) = fn(x0)
    		Dim x(4) As Double
    		x(1) = x0
    		Dim xx As Double     = 0.0
    		Dim d As Double      = d0
    		Dim count As Integer = 0
    		ind = -1
    
    		Do While count < max and ind < 0
    			x(3) = x(1) + d
    			f(3) = fn(x(3))
    			Dim k As Integer = 0
    			Do While k < max and f(3) <= f(1)
    				k   += 1
    				d   *= 2.0
    				x(0) = x(1)
    				f(0) = f(1)
    				x(1) = x(3)
    				f(1) = f(3)
    				x(3) = x(1) + d
    				f(3) = fn(x(3))
    			Loop
    						' 初期値が不適当
    			If k >= max
    				count = max
    			Else
    						' 3点の選択
    				Dim sw As Integer = 0
    				If k > 0
    					x(2) = x(3) - 0.5 * d
    					f(2) = fn(x(2))
    					Dim min As Integer = 0
    					For i1 As Integer = 0 To 3
    						If min <= 0 or f(i1) < f(min)
    							min = i1
    						End If
    						i1 += 1
    					Next
    					If min >= 2
    						For i1 As Integer = 0 To 2
    							x(i1) = x(i1+1)
    							f(i1) = f(i1+1)
    						Next
    					End If
    					sw = 1
    				Else
    					x(0) = x(1) - d0
    					f(0) = fn(x(0))
    					If f(0) > f(1)
    						x(2) = x(3)
    						f(2) = f(3)
    						sw   = 1
    					Else
    						x(1) = x(0)
    						f(1) = f(0)
    						d0   = -d0
    						d    = 2.0 * d0
    						k    = 1
    					End If
    				End If
    						' 収束?
    				If sw > 0
    					count += 1
    					Dim dl As Double = 0.5 * d * (f(2) - f(0)) / (f(0) - 2.0 * f(1) + f(2))
    					xx  = x(1) - dl
    					val = fn(xx)
    					If Math.Abs(dl) < eps
    						ind = count
    					Else
    						k  = 0
    						d0 = 0.5 * d
    						d  = d0
    						If val < f(1)
    							x(1) = xx
    							f(1) = val
    						End If
    					End If
    				End If
    			End If
    		Loop
    
    		Return xx
    
    	End Function
    
    End Module
    			

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