用离散相似法仿真程序重现输出Y1的数据和曲线，并与四阶龙格一库塔法比较精度。

参考程序：

P=[1 0.01 1 0;0 0.085 1 0.17;1 0.01 1 0;0 0.051 1 0.15;
    1 0.0067 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.0044 0];
WIJ=[1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1;4 8 -1;5 4 1;6 5 1;6 7 -0.212;7 6 1;8 6 1;9 7 1];
n=9;
Y0=1;
Yt0=[0 0 0 0 0 0 0 0 0];
h=0.001;
T=0;
T0=0;
Tf=5;
Nout=7;
A=P(:,1);B=P(:,2);
C=P(:,3);D=P(:,4);
m=length(WIJ(:,1));
W0=zeros(n,1);W=zeros(n,n);
for k=1:m
    if (WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);
    else W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
    end;
end;
for i=1:n
    if(A(i)==0);
        FI(i)=1;
        FIM(i)=h*C(i)/B(i);
        FIJ(i)=h*h*C(i)/B(i)/2;
        FIC(i)=1;FID(i)=0;
        if(D(i)~=0);
            FID(i)=D(i)/B(i);
        else
        end;
    else
        FI(i)=exp(-h*A(i)/B(i));
        FIM(i)=(1-FI(i))*C(i)/A(i);
        FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);
        FIC(i)=1;FID(i)=0;
        if(D(i)~=0);
            FIM(i)=(1-FI(i))*D(i)/A(i);
            FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i);
            FIC(i)=C(i)/D(i)-A(i)/B(i);
            FID(i)=D(i)/B(i);
        else
        end;
    end;
end;
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;
t=T0:h:Tf;N=length(t);
for k=1:N-1
    Ub=Uk;
    Uk=W*Y+W0*Y0;
    Udot=(Uk-Ub)/h;
    Uf=2*Uk-Ub;
    X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
    Y=FIC'.*X+FID'.*Uf;
    y=[y,Y(Nout)];
end;
plot(t,y)

参考结果：

结果比较：

离散相似法(蓝色)，四阶龙格一库塔法(红色)。

────────────
　　北京信息科技大学自动化学院“控制系统数字仿真与CAD”实验内容，仅供参考

08.05.26最后更新