最適化（黄金分割法）

黄金分割法

黄金分割法を使って関数の極小値を求めます。（後でもう少し説明を加えます）

close all
axis([1,5,1,10])
hold on

f=@(x) 3*(x-3).^2+2    %対象の関数

n=10; %ループ回数


%極小値を含む区間を決める
%(この方法は囲い込み法と呼ばれるらしい)

x(1)=10;  %初期値
h=1;    %ステップサイズ
x(2)=x(1)+h;
if f(x(2))>f(x(1))    %f(x(2))>f(x(1))ならx(1)とx(2)を入れ替え、ステップhの方向を変える
    change=x(2);
    x(2)=x(1);
    x(1)=change;
    h=-h;
end

for i=1:n
    h=2*h;
    x(3)=x(2)+h;
    if f(x(3))<=f(x(2))
        x(1)=x(2);
        x(2)=x(3);
    else
        break;
    end
end

if h>0
    a=x(1);
    b=x(3);
else
    a=x(3);
    b=x(1);
end

fprintf('[a b]=[%f %f]\n',a,b)  %決定された区間を出力


%黄金分割法
p=(5^(1/2)-1)/2;
c=a+(b-a).*p.^2;
d=a+(b-a).*p;

for i=1:n
    if f(c)<=f(d)
        b=d;
        d=c;
        c=a+(b-a).*p.^2;
        minx(i)=c;     %分割の様子をプロットしたいので結果を配列に格納しておきます
        miny(i)=f(c);
    else
        a=c;
        c=d;
        d=a+(b-a).*p;
        minx(i)=d;
        miny(i)=f(d);
    end
end

plot(minx,miny,'r')     %途中経過を赤ライン、最終結果を黒丸でプロット
plot(minx(n),miny(n),'ko')
fprintf('minimum f(%f) = %f\n',minx(n),miny(n)) %最終的な結果を表示


%関数のプロット
x=-10:0.001:10;
plot(x,f(x))

結果

EDU>> golden_section

f = 

    @(x)3*(x-3).^2+2

[a b]=[-4.000000 8.000000]
minimum f(3.026138) = 2.002050

黄金分割で使用した点を順に赤ラインでつないでいます。