实验四离散相似法数字仿真

1、实验四离散相似法数字仿真【实验目的】掌握离散相似法仿真方法【实验内容】用离散相似法仿真程序（参考课本sp4-4m）重现下图输出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

2、=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); FI

3、J(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; %积分环节系数求取完毕 if(D(i)=0);%若为积分比例，修正d 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); %若为惯性比例环节修正 c d FIM(i)=(1-FI(i)*D(i)/

4、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);%建立时间序列向量 n维for k=1:N-1 Ub=Uk;%保存当前一次输入值 Uk=W*Y+W0*Y0;%保存当前uk Udot=(Uk-Ub)/h;%求当前uk导数 Uf=2*Uk-Ub;%求下一步输入uf X=FI.*X+FIM.*

5、Uk+FIJ.*Udot;%点运算求取状态和输出 Y=FIC.*X+FID.*Uf; y=y,Y(Nout);end;figure(1)plot(t,y,g)hold以下运用龙格-库塔法实现：Q=B-D*w;%求出各阵Qn=inv(Q);R=C*w-A;V1=C*w0;Ab=Qn*R;b1=Qn*V1;Y=yt0;y=Y(nout);t=T0;N=round(Tf-T0)/h);%龙哥库塔法for i=1:N k1=Ab*Y+b1*y0; k2=Ab*(Y+h*k1/2)+b1*y0; k3=Ab*(Y+h*k2/2)+b1*y0; k4=Ab*(Y+h*k3)+b1*y0; Y=Y+h*(k1+2*k2+2*k3+k4)/6; y=y,Y(nout); T=T,t+h; t=t+h;end;figure(