代码:
function c7fun7_4
t=[2002 2003 2004 2005 2006 2007];
%X10=[269751 2710156 270202 270202 272770 2777314 282783]/100
%此处多一个数据，需要把多的删去
X10=[269751 2710156 270202 272770 2777314 282783]/100;
X20=[17675 250172 265415 239894 255345 260253]/100;
X30=[566.9 756 761.3 681 590 1236];
X40=[812 1106 1108 1008 1177 1251];
global AA UU;
U(4)=0;A(4,4)=0;
kk1=c7fun73(X10);A(1,1)=-kk1(1);U(1)=kk1(2);
kk2=c7fun72(X20,[X10;X30]);
A(2,1)=kk2(2);
A(2,2)=-kk2(1);
A(2,3)=kk2(3);
kk3=c7fun72(X20,X30);A(3,3)=-kk3(1);A(3,2)=kk3(2);
kk4=c7fun73(X40);A(4,4)=kk4(1);U(4)=kk4(2);
disp('估计出的系统预测模型系数-微分方程组系数')
AA=A
UU=U(:)
%EDIT BY PSL@CSU
%QQ：547423688
%Email：[email protected]
%——求算GM（1，N)参数函数
function a=c7fun72(X10,Xn0)
[n,m]=size(Xn0);
%——1次累加生成
s1=0;s2([1:n],1)=0;
for jj=1:1:m
X11(jj)=s1+X10(jj); Xn1(:,jj)=s2+Xn0(:,jj); s1=X11(jj); s2=Xn1(:,jj);
end
%——建立矩阵B、yn
for ii=1:1:m-1
B(ii)=-(X11(ii)+X11(ii+1))/2;
end
XN1=Xn1';B=[B(:),XN1(2:m,1:n)];
yn=X10([2:m])'; a=inv((B'*B))*B'*yn;
%——求算GM（1，1）参数函数
function au=c7fun73(X0)
m=length(X0);
%——1次累加生成
s1=0;
for jj=1:1:m; X1(jj)=s1+X0(jj); s1=X1(jj); end
%——建立矩阵B、y
for ii=1:1:m-1; B(ii)=-(X1(ii)+X1(ii+1))/2; end
B=[B(:),ones(m-1,1)];
y=X0([2:m])';au=inv((B'*B))*B'*y;