matlab代码：
%在workspace里输入x =[5999,5903,5848,5700,7884];gm1(x)   %其中5999,5903,5848,5700,7884可以换成已知的历史数据，无论几个都可以。function  gmcal=gm1(x)sizexd2 = size(x,2);%求数组长度k=0;for y1=x    k=k+1;    if k>1        x1(k)=x1(k-1)+x(k);        %累加生成        z1(k-1)=-0.5*(x1(k)+x1(k-1));           %z1维数减1，用于计算B        yn1(k-1)=x(k);    else        x1(k)=x(k);    endend%x1,z1,k,yn1sizez1=size(z1,2);%size(yn1);z2 = z1';z3 = ones(1,sizez1)';YN = yn1';   %转置%YNB=[z2 z3];au0=inv(B'*B)*B'*YN;au = au0';%B,au0,auafor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,ua %输出预测的  a u 和 u/a的值constant1 = x(1)-ua;afor1 = -afor;x1t1 = 'x1(t+1)';estr = 'exp';tstr = 't';leftbra = '(';rightbra = ')';%constant1,afor1,x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)%输出时间响应方程,也就是最终要求的灰色模型