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

实验四离散相似法数字仿真 实验四离散相似法数字仿真 PAGE 实验四离散相似法数字仿真 实验四离散相似法数字仿真 【实验目的】 掌握离散相似法仿真方法 【实验内容】 用离散相似法仿真程序（参考课本sp4-4m）重现下图输出Y1的数据和曲线，并与四阶龙格一库塔法比较精度。 【实验的程序】 P=[1 1 0;0 1 ;1 1 0;0 1 ; 1 70 0;1 0;0 1 130 0;1 0;1 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 ;7 6 1;8 6 1;9 7 1]; n=9; Y0=1; Yt0=[0 0 0 0 0 0 0 0 0]; h=; 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);%若为积分比例，修正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)/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.*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(2) [T,y];%输出 plot(T,y,r);%绘图 axis([0 5 0 350]) 实验的结果如下： 运用离散化的方法得到的结果 运用龙格-库塔方法得到的结果 【实验总结：】