情報学部 | 菅沼ホーム | 目次 | 索引 |
/*********************************************/ /* 多項式近似による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; }
/*********************************************/ /* 多項式近似による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; } }
<!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>
<?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; } ?>
#********************************************/ # 多項式近似による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])
# -*- 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]))
/*********************************************/ /* 多項式近似による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; } }
''''''''''''''''''''''''''''''''''''''''''''' ' 多項式近似による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
情報学部 | 菅沼ホーム | 目次 | 索引 |