请matlab高手帮忙看一下小弟编的BFGS算法程序

cwd198318

这是我编的一个BFGS程序，但是运行时好像是死循环，用Ctrl+C中断程序后出错信息为：
Error in ==> num2cell at 33
        c{i} = a(i);

Error in ==> sym.subs at 68
VaR = num2cell(sym([ '[' vars ']' ]));

Error in ==> cwdBFGS at 29
fz0=subs(f,{x1,x2},{x0(1,1),x0(2,1)});%计算函数值

Error in ==> cwdBFGS at 120
              p0=-g0;

Error in ==> cwdBFGS at 115
        if fz>=fz0   %进入算法第五步

小弟虚心求教，请高手帮忙指点一下。

%BFGS算法求解
%求梯度
syms x1 x2 t   
f=21.5+x1*sin(4.0*pi*x1)+x2*sin(20.0*pi*x2); %目标函数，
fx=diff(f,x1);  %对x求梯度
fy=diff(f,x2);  %对y求梯度
fxy=[fx;fy];
%初始化赋值
dim=2;    %维数
e=0.01;    % H的终止准则限


%算法第一步
x0=[1;5]; %给初值
fz0=subs(f,{x1,x2},{x0(1,1),x0(2,1)});%计算函数值
g0=subs(fxy,{x1,x2},{x0(1,1),x0(2,1)});%计算梯度值


%算法第二步
k=0;  %迭代次数
H0=[1,0;0,1];  %迭代矩阵H初值为单位矩阵
p0=-g0;
      
%算法第三步,一维搜索
m=1; %标志量
iter=0;
while (m)
    iter=iter+1;
    x=x0+t*p0;
    ft=subs(f,{x1,x2},{x(1,1),x(2,1)});
   
    %调用进退法算法，确定范围
    d0=0;%初始点
    h=1;%初始步长
    s=2; %加倍系数
    y0=subs(ft,t,d0);
    n=0;   
    p=1;%标志量   
    while (p)
        d1=d0+h;
        y1=subs(ft,t,d1);
        if y1<y0
            h=s*h;
            d=d0;
            d0=d1;
            y0=y1;
            n=n+1;
        elseif n==0;
                 h=-h;
                 d=d1;
        else
                p=0;
                break;
        end
    end   
    a=min(d,d1);%得到搜索区间
    b=max(d,d1);
    %0.618算法确定t的值
    r=0.618;
    h1=a+(1-r)*(b-a);%两个边界点
    u1=a+r*(b-a);
    y1=subs(ft,t,h1);%边界点的值
    y2=subs(ft,t,u1);
    while (abs(h1-u1)>0.01)
        if y1<y2
            a=u1;
            u1=h1;
            y2=y1;
            h1=a+(1-r)*(b-a);
            y1=subs(ft,t,h1);
        else
            a=h1;
            h1=u1;
            y1=y2;
            u1=a+r*(b-a);
            y2=subs(ft,t,u1);
        end
    end
    t=(h1+u1)/2; %得到最优值
    %一维搜索后使新的迭代点不超过x的取值范围,限幅处理
    x=x0+t*p0;
    x=double(x);
    if x(1,1)>12.1  
        x(1,1)=12.1;
    end
    if x(1,1)<-3.0
        x(1,1)=-3.0;
    end
    if x(2,1)>5.8
        x(2,1)=5.8;
    end
    if x(2,1)<4.1
        x(2,1)=4.1;
    end
    %计算f(k+1)和g(k+1)的值
    fz=subs(f,{x1,x2},{x(1,1),x(2,1)});%计算新的值
    g=subs(fxy,{x1,x2},{x(1,1),x(2,1)});% 计算新的梯度
   
    %算法第四步
    while(norm(g(:,1),2)>=e) %判断H是否满足终止条件,计算2－范数
        if fz>=fz0   %进入算法第五步
              x0=x0;    %迭代后函数值没有下降，重新迭代
              fz0=fz0;
              g0=g0;
              H0=[1,0;0,1]; %迭代矩阵H初值为单位矩阵
              p0=-g0;
              m=1;
              k=0;
              continue; %转算法第三步
        elseif k==dim  %进入算法第六步
                  x0=x;
                  fz0=fz;
                  g0=g;
                  H0=[1,0;0,1]; %迭代矩阵H初值为单位矩阵
                  p0=-g0;
                  m=1;
                  k=0;
                  continue; %转算法第三步
        else                %进入算法第七步
              y=g-g0;
              s=x-x0;
              H=H0+[(1+y'*H0*y/(s'*y))*s*s'-H0*y*s'-s*y'*H0]/(s'*y);
              p=-H*g;
              H0=H;
              p0=p;
              fz0=fz;
              g0=g;
              x0=x;
              k=k+1;
              continue; %转算法第三步
        end
    end
    m=0;
end
xpso=x %输出x,f的值，BFGS算法结束
fpso=fz        
iter

射天狼

网上可以搜索到bggs的matlab源程序的.