拟牛顿 DFP matlab

function sevnn
x=[1,0]';
[x,val]=dfp('fun','gfun',x)
end

function f=fun(x) 
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2; 
end
 
function g=gfun(x) 
g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]'; 
end

function He=Hess(x) 
He=[1200*x(1)^2-400*x(2)+2, -400*x(1); -400*x(1), 200 ];
end

function [x,val]=dfp(fun,gfun,x0) 
%DFP算法
maxk=100;
rho=0.55;sigma=0.4; 
epsilon=1e-5; 
k=0; n=length(x0); 
Hk=inv(feval('Hess',x0)); %Hk=eye(n); 
while(k<maxk) 
    gk=feval(gfun,x0); %计算梯度 
    if(norm(gk)<epsilon), break; end %检验终止准则
    dk=-Hk*gk; 
    m=0; mk=0; 
    while(m<15) % 用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk) 
            mk=m; break; 
        end
        m=m+1; 
    end
    %DFP校正 
    x=x0+rho^mk*dk; 
    sk=x-x0; 
    yk=feval(gfun,x)-gk; 
    
    Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk); 
    
    k=k+1; 
    x0=x;
end
val=feval(fun,x0);
end