%共轭梯度法求函数的极小值
function [minf,optPoint]=conjGM(objFun)

%极小值minf

%最优解optPoint
syms w x y z;
v=[w;x;y;z];%变量向量,对应于目标函数objFun中的变量

f=objFun(v);
gradf=gradient(f);%函数的梯度

epsilon=1e-3;%最小误差
x1=[0 0 0 0]';%初始点

yj=x1;
dj=-subs(gradf,v,yj);%搜索方向

k=0;%统计循环次数
while norm(subs(gradf,v,yj))>=epsilon
    t=formula(objFun(yj+x*dj));
    lambdaf=@(x)(t);%关于lambda(步长)函数
    [~, lambda]=newton(lambdaf,0);%一维搜索最优的步长
    yj1=yj+lambda*dj;%新搜索到的点
    gj=subs(gradf,v,yj);%上一点的梯度
    gj1=subs(gradf,v,yj1);%新搜索到点的梯度  
    beta=gj1'*gj1/(gj'*gj);
    dj1=-gj1+beta*dj;
    yj=yj1;
    %如果与新的负梯度方向的夹角大于某个值，则选择负梯度作为新的搜索方向
    if acos(dot(-gj1,dj1)/(norm(-gj1)*norm(dj1)))*180/pi>75
        dj=-gj1;
    else
        dj=dj1;
    end   
    k=k+1;
    disp([yj' subs(f,v,yj)  norm(subs(gradf,v,yj))])
end
disp(k)
minf=subs(f,v,yj);%函数的最小值
optPoint=yj;
end

%%牛顿法一维搜索
function [optV,optPoint]=newton(objFun,x0)
%目标函数objFun为一元函数，x0为起始点
epsilon=1e-5;%设定最小误差
syms x;
f=objFun(x);
df=diff(f,'x');%一阶导数
dff=diff(f,2);%二阶导数
xk=x0;
t=subs(df,x,xk);
while  abs(t)>=epsilon
    xk=x0-subs(df,x,x0)/subs(dff,x,x0);
    t=subs(df,x,xk);
    x0=xk;
end
optV=objFun(xk);
optPoint=xk;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%示例运行结果如下：

>> f=@(x)((x(1)-2).^4+(x(2)+3).^4+(x(3)-5).^4+(x(4)-7).^4+20);
>> [minf,optPoint]=conjGM(f)
      0.19899     -0.67159       3.1092       8.5317       78.199       63.509
       1.5455      -3.5854       4.7173       7.8538       20.698       2.6443
       1.7232      -3.2802       4.7766       6.7818       20.017      0.13657
       2.0485      -2.9399       4.9468       6.9138           20     0.002808
       2.0485      -2.9399       4.9468       6.9138           20     0.002808
       2.0431      -2.9504       4.9544       6.9456           20   0.00094904
     6
minf =
           20
optPoint =
       2.0431
      -2.9504
       4.9544
       6.9456

%%%%%%%%%%%%%%%

注意：

1）语句“v=[w;x;y;z];”中向量v的维数一定要同目标函数objFun中变量（向量）的维数一致。根据实际情况调整。所给实例中x为4维向量，故而v的维数也为4.

2）算法运行时间同初始点的选择及运算精度有关。