%DFP法
clear;
syms x1 x2;
fun=fun1;
x0=[0;0];%設(shè)置初值點(diǎn)
x=x0;
Dfx=[diff(fun,x1);diff(fun,x2)];
g=subs(Dfx,[x1 x2],x);
fx=subs(fun,[x1 x2],x);
d=-g;
e=1e-6;%設(shè)置精度
num=0;
alpha=1;
sigma=0.5;
rho=0.5;
%開始用Armijo型線搜索方法求步長(zhǎng)，并取σ=0.5，ρ=0.5
while(subs(fun,[x1 x2],x+alpha*d)>(subs(fun,[x1 x2],x)+sigma*alpha*g'*d))
    alpha=rho*alpha;
end
x_next=x+alpha*d;
num=num+1;
fx=subs(fun,[x1 x2],x);
g_next=subs(Dfx,[x1 x2],x_next);
g_fan=sqrt(g'*g);
k=1;
H=eye(2);
while(g_fan>e)
    if(k==2)
        H_next=eye(2);
        k=1;
    else
        p=x_next-x;
        q=g_next-g;
        H_next=H+(p*p')/(p'*q)-(H*q*q'*H)/(q'*H*q);
        k=k+1;
    end
    H=H_next;
    g=g_next;
    x=x_next;
    d=-H*g;
    alpha=1;
    sigma=0.5;
    rho=0.5;
    %開始用Armijo型線搜索方法求步長(zhǎng)，并取σ=0.5，ρ=0.5
    while(subs(fun,[x1 x2],x+alpha*d)>(subs(fun,[x1 x2],x)+sigma*alpha*g'*d))
        alpha=rho*alpha;
    end
    x_next=x+alpha*d;
    g_next=subs(Dfx,[x1 x2],x_next);
    g_fan=sqrt(g_next'*g_next);
    num=num+1;
    fx=subs(fun,[x1 x2],x);
end
x=x_next;
fx=subs(fun,[x1 x2],x);
fprintf('\nDFP法\n結(jié)果：\n優(yōu)化精度為ε=%d\n初指點(diǎn)x0=[%d %d]\n最優(yōu)點(diǎn)x=[ %d %d ] \n最優(yōu)值fx=%d \n梯度模||g(x)||=%d \n迭代次數(shù)k=%d\n',e,x0(1),x0(2),x(1),x(2),fx,g_fan,num);
[x,y]=meshgrid(-3:0.1:3,-3:0.1:3);
z=(x.^2-y+1).^2+(x-cos(pi/2*y)).^2;
meshc(x,y,z);
xlabel('x1');
ylabel('x2');