黃金分割法求極小值

    功能：用黃金分割法求f(x)在區間[a,b]上的近似極小值。當且僅當f(x)在[a,b]上為單峰時次方法適用

----------------------------------------------------------------------

function [S,E,G]=golden(f,a,b,delta,epsilon)

% input  - f is the object function input as a string 'f'

%        - a and b are the end points of the interval

%        - delta is the tolerance for the abscissas

%        - epsilon is the tolerance for the ordinates

% output - S=(p,yp) contains the abscissa p and the ordinate yp of  
%          the minimum

%        - E=(dp,dy) contains the error bounds for p and yp

%        -G is an n*4 matrix : the kth row contains

%        -[ak ck dk bk]: the values of a,c,d and b at the kth interation

r1= (sqrt(5)-1)/2;

r2=r1^2;

h=b-a;

ya=feval(f,a);

yb=feval(f,b);

c=a+r2*h;

d=a+r1*h;

yc=feval(f,c);

yd=feval(f,d);

k=1;

A(k)=a;B(k)=b;C(k)=c;D(k)=d;

while (abs(yb-ya)>epsilon) | (h>delta)

      k=k+1;

      if (yc<yd)

         b=d;
   
         yb=yd;

         d=c;

         yd=yc;

         h=b-a;

         c=a+r2*h;

         yc=feval(f,c);

      else 

         a=c;

         ya=yc;

         c=d;

         yc=yd;

         h=b-a;

         d=a+r1*h;

         ya=feval(f,d);

      end

      A(k)=a;B(k)=b;C(k)=c;D(k)=d;

end

dp=abs(b-a);

dy=abs(yb-ya);

p=a;

yp=ya;

if (yb<ya)

   p=b;

   yp=yb;

end;

G=[A' C' D' B'];

S=[p yp];

E=[dp dy];