黄金分割法（保留打卡matlab，有需要自己复制）

% function [s, phis,k,G,E]=zzzz(phi, a,b,delta, epsilon) 
% 输入:phi是目标函数，a, b是搜索区间的两个端点
%    delta,epsilon分别是自变量和函数值的容许误差
% 输出:s, phis分别是近似极小点和极小值，G是nx4矩阵，
%    其第k行分别是a，p,q, b的第k次迭代值[ak,pk,gk, bk], 
%    E=[ds, dphi],分别是s和phis的误差限
% t=(sqrt(5)-1)/2; 
% h=b-a;
% phia=feval(phi,a); phib=feval(phi,b);
% p=a+(1-t)*h;    q=a+t*h; 
% phip=feval(phi,p);phiq=feval(phi,q); 
% k=1;  G(k,:)=[a,p,q,b]; 
% while(abs(phib-phia)>epsilon)||(h>delta)
%       两端点函数值之差的绝对值是否小于预设的函数值误差或者
%       搜索长度h是否小于预设的自变量误差限Δ    （如果不满足终止条件，继续迭代）
%          if(phip<phiq) %如果phip更小，说明目标函数在区间可能存在局部极小值,
%              右端点b更新为q，
%              左端点p更新为原q,同时更新phiq。
%              根据更新后的端点计算新的黄金分割点p,并计算目标函数phip
%             b=q; phib=phiq; q=p; phiq=phip; 
%             h=b-a; p=a+(1-t)*h;phip=feval (phi,p); 
%          else %否则，如果phiq更小,左端点a更新为p,同时更新phia
%               右端点q更新为原p,同时更新phiq.然后根据更新后的端点计算新的黄金分割点q，
%               并计算目标函数phiq
% 
%             a=p; phia=phip; p=q; phip=phiq;
%             h=b-a; q=a+t*h; phiq=feval(phi,q);
%         end 
%         k=k+1; G(k,:)=[a,p,q,b];
%         增加迭代计数器k，并将当前迭代过程中的四个关键点存储到矩阵G的第k行
% end 
% 
% ds=abs(b-a); dphi=abs(phib-phia);   
% if(phip<=phiq) 
%    s=p; phis=phip; 
% else 
%    s=q; phis=phiq; 
% end 
% E=[ds, dphi];

% 
% function [s,phis,ds,dphi,S]=zzzz(phi,a,~,delta,epsilon) %输入:phi 是目标函数，a和b是搜索区间的端点% delta,epsilon是容许误差%输出:s是近似极小点，phis是对应的近似极小值；k是选代次数% ds是迭代终止时的步长，dphi是Iphi(s1)-phi(s)l;S是迭代向量
% s0=a; maxj=20; maxk=30; big=1e6; err=1;k=1; 
% S(k)=s0; cond=0; h=1; ds=0.00001; 
% if (abs(s0)>1e4),h=abs(s0)*(1e-4);end
% while (k<maxk && err>epsilon &&cond~=5) 
%     f1=(feval(phi,s0+ds)-feval(phi, s0-ds))/(2*ds);
%     if(f1>0),h=-abs(h);end 
%     s1=s0+h;     s2=s0+2*h;   bars=s0;
%     phi0=feval(phi,s0); phi1=feval(phi,s1); 
%     phi2=feval(phi,s2); barphi=phi0; cond=0;
%     j=0;%确定h使得phi1<phi0且phi1<phi2 
%     while(j<maxj&&abs(h)>delta&&cond==0) 
%         if (phi0<=phi1) 
%             s2=s1;phi2=phi1;h=0.5*h; 
%             s1=s0+h;phi1=feval(phi,s1); 
%         else if (phi2<phi1) 
%                 s1=s2;phi1=phi2;h=2*h; 
%                 s2=s0+2*h;phi2=feval(phi,s2); 
%             else 
%                 cond=-1; 
%              end 
%          end 
%          j=j+1; 
%          if(abs(h)>big||abs(s0)>big),cond=5;end 
%     end 
%     if (cond==5) 
%          bars=s1;barphi=feval(phi,s1);
%     else %二次插值求phis
%         d=2*(2*phi1-phi0-phi2); 
%         if(d<0)
%             barh=h*(4*phi1-3*phi0-phi2)/d;
%        else 
%             barh=h/3; cond=4; 
%        end 
%        bars=s0+barh;  barphi=feval(phi,bars); 
%        h=abs(h);  h0=abs(barh); 
%        h1=abs(barh-h);   h2=abs(barh-2*h); %确定下一次迭代的值
%        if(h0<h),h=h0; end 
%        if(h1<h),h=h1; end 
%        if(h2<h),h=h2; end 
%        if(h==0),h=barh; end 
%        if(h<delta),cond=1; end 
%        if(abs(h)>big||abs(bars)>big),cond=5; end 
%        err=abs(phi1-barphi); 
%        s0=bars; k=k+1;S(k)=s0; 
%      end 
%      if(cond==2&&h<delta),cond=3;end 
% end 
% s=s0;phis=feval(phi,s); 
% ds=h; dphi=err; 
% who
% whos


function mk=zzzz(xk, dk) 
beta=0.5;  sigma=0.2; 
m=0;mmax=20; 
while (m<=mmax) 
    if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) 
         mk=m; break; 
    end 
    m=m+1; 
end 
alpha=beta^mk; 
newxk=xk+alpha*dk; 
fk = fun(xk); 
newfk = fun(newxk); 
function f=fun(x) 
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2; %梯度
function gf=gfun(x) 
gf=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';
who
whos