grads.m

最优化理论与算法（第2版）这本书中的课后作业。用C++实现的一些具体算法。

function [p0,y0,err]=grads(F,G,p0,max1,delta,epsilon,show)

%input    F is the object function input as a string'F';
%        G=-(1/norm(grad F)*grad f;the search direction input as a string'G';
%        p0 is the initial starting point;
%        max1 is the maximun number of iterations;
%        delta is the tolerance for hmin in the single parameter
%              minimization in the search didection;
%        epsilon is the tolerance for the error in y0;
%        show; if show==1 the iterations are displayed;
%output  p0 is the point for the minimum;
%        y0 is the function value F(p0);
%        err is the error bound for y0;
%        p is a vector containing the iterations

if nargin==5,show=0;end
[mm,n]=size(p0);
maxj=10;big=1e8;h=1;
p=zeros(maxj,n+1);
len=norm(p0);
y0=feval(F,p0);
if (len>e4),h=len/1e4;end
err=1;cnt=0;cond=0;
p(cnt+1,:)=[p0 y0];
while(cnt<max1&cond~=5&(h>delta|err>epsilon))
    %computer search direction
    s=feval(G,p0);
    %start single parameter quadratic minimization
    p1=p0+h*s;
    p2=p0+2*h*s;
    y1=feval(F,p1);
    y2=feval(F,p2);
    cond=0;j=0;
    while(j<maxj&cond==0)
        len=norm(p0);
        if (y0<y1)
            p2=p1;
            y2=y1;
            h=h/2;
            p1=p0+h*s;
            y1=feval(F,p1);
        else
            if(y2<y1)
                p1=p2;
                y1=y2;
                h=2*h;
                p2=p0+2*h*s;
                y2=feval(F,p2);
            else
                cond=-1;
            end
        end
        j=j+1;
        if(h<delta),cond=1;end
        if(abs(h)>big|len>big),cond=5;end
    end
    if(cond==5)
        pmin=p1;
        ymin=y1;
    else
        d=4*y1-2*y0-2*y2;
        if(d<0)
            hmin=h*(4*y1-3*y0-y2)/d;
        else
            cond=4;
            hmin=h/3;
        end
        
        %construct the next point
        pmin=p0+hmin*s;
        ymin=feval(F,pmin);
        
        %determine magnitude of next h
        h0=abs(hmin);
        h1=abs(hmin-h);
        h2=abs(hmin-2*h);
        
        if(h0<h),h=h0;end
        if(h1<h),h=h1;end
        if(h2<h),h=h2;end
        if(h==0),h=hmin;end
        if(h<delta),cond=1;end
        
        %termination test for minimization
        e0=abs(y0-ymin);
        e1=abs(y1-ymin);
        e2=abs(y2-ymin);
        if(e0~=0&e0<err),err=e0;end
        if(e1~=0&e1<err),err=e1;end
        if(e2~=0&e0<err),err=e2;end
        if(e0==0&e1==0&e2==0),err=0;end
        if(err<epsilon),cond=2;end
        if(cond==2&h<delta),cond=3;end
    end
    cnt=cnt+1;
    p(cnt+1,:)=[pmin ymin];
    p0=pmin;
    y0=ymin;
end
if(show==1)
    disp(p);
end