系统仿真技术实验指导书(2024年春季实验1-6).docx

系统仿真技术试验指导书次新明,栩物编写中南高校信息科学与工程学院2014年5月18日书目试验一MAT1.AB中矩阵与多项式的基本运算1试验二MAT1.AB绘图吩咐2试验三MAT1.AB程序设计3试验四Matlab的符号计算与Simulink的运用5试验五MAT1.AB在限制系统分析中的应用8试验六连续系统数字仿真的基本算法15试验一MAT1.AB中矩阵与多项式的基本运算一、试验任务1 .了解MAT1.AB吩咐窗口和程序文件的调用。2 .熟识如下MAT1.AB的基本运算：矩阵的产生、数据的输入、相关元素的显示；矩阵的加法、乘法、左除、右除；特别矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；多项式的运算：多项式求根、多项式之间的乘除。二、基本吩咐训练1. eye(m)2. one(n)ones(m,n)3. zeros(m,n)4. rand(m,n)5. diag(v)6. AB、AB>inv（八）*B、B*inv（八）7. roots(p)8. poly9. conv、deconv10. A*B与A.*B的区分11. who与whos的运用12. dispsize（八）、Iength（八）的运用三、试验要求依据试验内容和相关吩咐进行试验，自拟输入元素，将上述各吩咐的输入和输出结果写成一、试验任务熟识MAT1.AB基本绘图吩咐，驾驭如下绘图方法：1 .坐标系的选择、图形的绘制；2 .图形注解（题目、标号、说明、分格线）的加入;3 .图形线型、符号、颜色的选取。二、基本吩咐训练1. plot2.loglog3.Semilogx4.semiIogy5.polar6.title7.xlabel8.ylabel9.text10.grid11.bar12.stairs13. contour三、试验举例1. t=0:pi/360:2*pi*22/3;x=93*cos(t)+36*cos(t*4.15);y=93*sin(t)+36*sin(t*4.15);plot(yzx),grid;2. t=0:0.05:100;x=t;y=2*t;z=sin(2*t);plot3(x,y,z,b:,)3. t=0:pi/20:2*pi;y=sin(x);stairs(xzy)4. th=pi/200:pi/200:2*pi,;r=cos(2*th);polar(th,r),grid5. th=0:pi/10:2*pi;x=exp(j*th);plot(real(x),imag(x),r*,);grid;四、试验要求在两种或两种以上坐标系绘制35个图形,要有颜色、图形种类、注解的不同一、试验任务1 .熟识MAT1.AB程序设计的方法和思路；2 .驾驭循环、分支语句的编写，学会运用lookfor、help吩咐。二、程序举例1 .计算l1000之内的斐波那契亚数列f=l,l;i=l;whilef(i)+f(i+l)<1000f(i+2)=f(i)+f(i+l);i=i+l;endf,i2 .m=3;n=4;fori=l:mforj=l:na(i,j)=l(i+j-l);endendformatrata3 .m=3;n=4;fori=l:mforj=l:na(izj)=l(i+j-l);endenda请比较程序2与程序3的区分4 .x=input('请输入X的值：')；ifx=10y=cos(x+l)+sqrt(x*x+l);elsey=x*sqrt(x+sqrt(x);endy5 .去掉多项式或数列开头的零项p=O001302009;fori=l:length(p),ifp(1)=0zp=p(2:length();end;end;P6 .建立MAT1.AB的函数文件，程序代码如下，以文件名ex2_4.m存盘functionf=ffibno(n)%ffibno计算斐波那契亚数列的函数文件%n可取随意自然数%程序如下f=l,U;i=l；whilef(i)+f(i+l)<nf(i+2)=f(i)+f(i+l);i=i+l;end输入完毕后在MAT1.AB的吩咐窗口输入ex2_4(200),得到运行结果。在MAT1.AB的吩咐窗口输入100kforffibno,得到结果：ex2_4.m:%ffibno计算斐波那契亚数列的函数文件在MAT1.AB的吩咐窗口输入helpex2_4,得到结果：ffibno计算斐波那契亚数列的函数文件n可取随意自然数程序如下三、程序设计题用一个MAT1.AB语言编写一个程序：输入一个自然数，推断它是否是素数，假如是，输出"Itisoneprime",假如不是，输出"Itisnotoneprime."。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关吩咐来限制是否接着进行素数的推断。假如操作者希望停止这种推断，则可以退出程序。假如所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入6,因为6不是素数。则程序中除了有“Itisnotoneprime”的结论外，还应有："6=2*3”的说明。四、试验要求1 .参照已知程序，改动程序中的参数和输入量，验证输出结果。2 .运用100kfOr、help吩咐，验证输出结果3 .试验结果写成试验报告三(全部试验完成后交试验报告)。试验四matlab的符号计算与Simulink的运用一、试验任务1 .驾驭MAT1.AB符号计算的特点和常用基本吩咐;2 .驾驭SIMU1.INK的运用。二、程序举例1. 求矩阵对应的行列式和特征根a=sym(,allal2;a21a22,);da=det（八）ea=eig（八）2. 求方程的解(包括精确解和肯定精度的解)rl=solve(,xa2-x-1')rv=vpa(rl)rv4=vpa(rl,4)rv20=vpa(rl,20)3. a=sym(,a,);b=sym(,b,);c=sym(,c,);d=sym(,d');%定义4个符号变量w=10;x=5;y=-8;Z=Il;A=a,b;c,dB=w,x;y,zdet（八）det(B)生定义4个数值变量生建立符号矩阵A务建立数值矩阵B%计算符号矩阵A的行列式%计算数值矩阵B的行列式4 .symsXy;s=(-7*xA2-8*yA2)*(-xa2+3*ya2);expand(三)名对S绽开collect(s,x)告对S按变量X合并同类项(无同类项)factorfans)%*ans分解因式5 .对方程AX=b求解A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=Iinsolve(A,b)皆调用IinSOIVe函数求解Ab考用另一种方法求解6 .对方程组求解alI*xl+al2*x2+al3*x3=bla21*xl+a22*x2÷a23*x3=b2a31*xl+a32*x2+a33*x3=b3symsallal2al3a21a22a23a31a32a33blb2b3;A=all,al2,al3;a21,a22za23;a31,a32,a33;b=bl;b2;b3;X=Iinsolve(A,b)传调用IinSolVe函数求的解XX=Ab告用左除运算求解7 .symsabtXyz;f=sqrt(l+exp(x);diff(f)f=x*cos(X);diff(f,x,2)diff(f,x,3)fl=a*cos(t);f2=b*sin(t);diff(f2)diff(fl)学未指定求导变量和阶数，按缺省规则处理考求f对X的二阶导数%求f对X的三阶导数行按参数方程求导公式求y对X的导数三、SIMU1.INK的运用1.在吩咐窗口中输入：SimUlink(回车)得到如下SimUlink模块:2.双击打开各模块，选择合适子模块构造限制系统，例如:3.双击各子模块可修改其参数，选择SimUEion菜单下的Start吩咐运行仿真，在示波器(Scope)中视察结果。四、试验要求1.符号计算部分：参照所给示例，自拟输入量，验证求矩阵的特征根、行列式及方程求解吩咐2.SimUIink部分：逐一熟识各模块的运用，对下面的随动系统模型进行仿真,试验报告中包含：连好的系统模型及用SCOPe观测的结果其中：R(三)为阶跃输入，C(三)为输出15+0.025.2G,(三)=Gl(三)=1105+0.00225(5+0.5)3.试验结果写成试验报告四备注：全部试验完成后，在试验报告的封面注明：课程、班级、姓名、学号等信息。按时交试验报告。试验五MAT1.AB在限制系统分析中的应用一、试验任务1 .驾驭MAT1.AB在限制系统时间响应分析中的应用;2 .驾驭MAT1.AB在系统根轨迹分析中的应用；3 .驾驭MAT1.AB限制系统频率分析中的应用；4 .驾驭MAT1.AB在限制系统稳定性分析中的应用二、基本吩咐1.step2.impulse3.initial4.Isim5.rlocfind6. bode7.margin8.nyquist9.Nichols10.cloop三、程序举例1 .求下面系统的单位阶跃响应G(三)=-5+5+4%程序如下：num=4;den=l,1,4;step(num,den)y,X,t=step(num,den);tp=spline(y,t,max(y)%计算峰值时间max(y)%计算峰值2 .求如下系统的单位阶跃响应-XUO1J1.-6-5FiJ1.2.Tw%程序如下：a=l,l;-6,-5;b=0;l;c=l,0;d=0;y,x=step(a,b,c,d);piot(y)3 .求下面系统的单位脉冲响应:%程序如下：num=4;den=l,1,4；impulse(num,den)4 .已知二阶系统的状态方程为:求系统的零输入响应和脉冲响应。%程序如下：a=O,l;-10,-2;b=O;l;c=l,O;d=O;x=l,0;subplot(1,2,1);initial(a,b,c,d,x)subplot(l,2,2);impulse(a,b,c,d)5：系统传递函数为:G(三)输入正弦信号时，视察输出信号的相位差。%程序如下：num=l;den=l,1；t=0:0.01:10;u=sin(2*t);holdonplot(t,u,T)lsim(num,den,u,t)6 .有一二阶系统，求出周期为4秒的方波的输出响应G(三)2s+5s+1s+2s÷3%程序如下:num=251;den=ll23;t=(0:.l:10);period=4;u=(rem(t,period)>=period./2);%Srem函数功能lsim(num,den,u,t);7 .已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性G(三)=k(s+2)(d+4$+3)2%程序如下：num=l2;denl=143;den=conv(den1.den1);figure(l)rlocus(num,den)k,p=rlocfind(num,den)figure(2)k=55;numl=k*l2;den=l43;den1=conv(den,den);num,den=cloop(numl,den1,-1);impulse(num,den)title('impulseresponse(k=55),)figure(3)k=56;numl=k*12;den=l43;den1=conv(den,den);num,den=cloop(numl,den1,-1);impulse(num,den)title('impulseresponse(k=56),)8 .作如下系统的bode图G(三)=-5+45÷115+7%程序如下：n=l,l;d=l,4,ll,7;bode(n,d)9 .系统传函如下G(三)='+I产(s+2)3求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应%程序如下：num=1;den=conv(12,conv(l2,12J);w=logspace(-l,2);t=0.5;ml,pl=bode(num,den,2);p1=p1-t*w*180/pi;n2,d2=pade(t,4);numt=conv(n2,num);dent=(conv(den,d2);m2,p2=bode(numt,dent,w);subplot(2,1,1);semilogx(w,20*log10(ml),w,20*log10(m2),'g-');gridon;title(,bodeplot')5xlabel(,frequency')jylabel(,gain');subplot(2,1,2);semilogx(w,p1,w,p2,g-,)igridon;XlabeI('frequency');ylabel('phase');10 .已知系统模型为、3.5G(三)=-s+2s+3s+2求它的幅值裕度和相角裕度%程序如下：n=3.5;d=l232;Gm,Pm,Wcg,Wcp=margin(n,d)11 .二阶系统为：GG)=3；s2+2ns+l令Wn=I,分别作出=2,1,0.707,0.5时的nyquist曲线。%程序如下：n=UJ；dl=l,4,l;d2=l,2,l;d3=l,1.414,1;d4=l,l,l;nyquist(n,dl);holdonnyquist(n,d2);nyquist(n,d3);nyquist(n,d4);12 .已知系统的开环传递函数为、100OG(三)=;G+3s+2)*(s+5)绘制系统的NyqUSit图，并探讨系统的稳定性%程序如下：G=tf(l000,conv(l,3,2,l,5);nyquist(G);axisfsquare')13 .分别由W的自动变量和人工变量作下列系统的nyquist曲线：G(三)=?s(s+l)%程序如下：n=l;d=l,1,0;nyquist(n,d);%自动变量n=l;d=l,1,0;w=0.5:0.1:3;nyquist(n,d,w);%人工变量14 .一多环系统，其结构图如下，运用NyqUiSt频率曲线推断系统的稳定性。(0.85S+1)(0.255+1)(0.0625s+1)%程序如下：kl=16.7/0.0125;zl=0;pl=-1.25-4-16;num1,den1=zp2tf(zl,pl,kl);num,den=cloop(numl,den1);z,p,k=tf2zp(num,den);pfigure(l)nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2);15 .已知系统为：G(三)=?s(s+l)作该系统的nichols曲线。%程序如下：n=l;d=l,l,0;ngrid(，new，);nichols(n,d);16 .已知系统的开环传递函数为：G(三)=s(s+l)(s+2)当k=2时，分别作nichols曲线和波特图。%程序如下：num=l;den=conv(conv(l0,ll),0.51);subplot(1,2,1);nichols(num,den);grid;%nichols曲线subplot(1,2,2);g=tf(num,den);bode(feedback(g,l,-l)；grid;%波特图17.系统的开环传递函数为:G(三)=s(s+l)(s+2)分别确定k=2和k=10时闭环系统的稳定性。%程序如下：dl=l,3,2,0;nl=2;ncl,dcl=cloop(nl,dl,-1);roots(dc1)d2=dl;n2=10;nc2,dc2=cloop(n2,d2,-l);roots(dc2)18.系统的状态方程为：试确定系统的稳定性。%程序如下：a=-4,-3,0;1,0,0;0,1,0;b=l;0;0;c=0,l,2;d=0;eig（八）%求特征根rank(ctrb(a,b)四、试验要求1 .参照已知程序，改动程序中的参数和输入量，验证输出结果。2 .驾驭相关吩咐的用法3 .试验结果写成试验报告五(全部试验完成后交试验报告)。试验六连续系统数字仿真的基本算法一、试验任务1 .理解欧拉法和龙格库塔法的基本思想；2 .理解数值积分算法的计算精度、速度、稳定性与步长的关系;二、程序举例1.¾h=0.2,试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。W)一今o)=l注：解析解：y(t)=+27%置自变量初值%置解析解和数值解的初值%步长%程序如下cleart(D=O;y(l)=l;y_euler(l)=l;y_rk2(l)=l;y_rk4(l)=l;h=0.2;%求解析解fork=l:5t(k+l)=t(k)+h;y(k+1)=sqrt(1+2*t(k+1);end%利用欧拉法求解fork=1:5y_euler(k+l)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);end%利用RK2法求解fork=l:5k1=y_rk2(k)-2*t(k)/y_rk2(k);k2=(y_rk2(k)+h*kl)-2*(t(k)+h)/(y_rk2(k)+h*kl);y_rk2(k+l)=y_rk2(k)+h*(k1+k2)2;end%利用RK4法求解fork=l:5k1=y_rk4(k)-2*t(k)/y_rk4(k);k2=(yJrk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2);k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2);k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+h*k3);y_rk4(k+l)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)6;end%输出结果disp(,时间解析解欧拉法RK2法RK4法)yt=t',y',y-euler,y_rk2;y_rk4'disp(yt)程序运行结果如下：时间解析解欧拉法RK2法RK4法01.00001.00001.00001.00000.20001.18321.20001.18671.18320.40001.34161.37331.34831.34170.60001.48321.53151.49371.48330.80001.61251.68111.62791.61251.00001.73211.82691.75421.73212.考虑如下二阶系统：y(f)+2R)V)+y(f)=0在010上的数字仿真解(已知：y(0)=100,y(0)=0),并将不同步长下的仿真结果与解析解进行精度比较。说明：已知该微分方程的解析解分别为：y(r)=100cosr(当R=O)/XInnTV51005T.百*ZiizD八y(t)=Q0e2cos/+-e2sinr(R=0.5)采纳RK4法进行计算,、选择状态变量：=y¼=j则有如下状态空间模型及初值条件x1=x2xl(0)=100x2=-x1-IRx2x2(0)=0F"采纳RK4法进行计算。程序如下:clearh=inputC请输入步长h-);M=round(10h);t(l)=0;y_0(l)=100;y_05(l)=100;对应于为R=O和R=0.5)xl(l)=100;x2(l)=0;y_rk4_0(l)=xl(l);y_rk4_05(l)=xl(l);%求解析解fork=l:Mt(k+l)=t(k)+h;y_O(k+1)=100*cos(t(k+1);%输入步长%置总计算步数%置自变量初值%置解析解的初始值（y_0和y_05分别%置状态向量初值%置数值解的初值y.05(k+1)=1OO*exp(-t(k÷1)2).*cos(sqrt(3)2*t(k+1)+100*sqrt(3)3*exp(-t(k+l)2).*sin(sqrt(3)2*t(k+l);end%利用RK4法求解%R=Ofork=l:Mkll=x2(k);kl2=-xl(k);k21=x2(k)+h*k12/2;k22=-(xl(k)+h*kl1/2);k3I=x2(k)+h*k222;k32=-(xl(k)+h*k21/2);k41=x2(k)+h*k32;k42=-(xl(k)+h*k31);xl(k+l)=xl(k)+h*(k11+2*k21+2*k31+k41)/6;x2(k+l)=x2(k)+h*(k12+2*k22+2*k32+k42)6;y_rk4_0(k+l)=xl(k+l);end%R=0.5fork=l:Mkll=x2(k);k12=-xl(k)-x2(k);k21=x2(k)+h*kl22;k22=-(xl(k)+h*k1l2)-(x2(k)+h*k122);k31=x2(k)+h*k222;k32=-(xl(k)+h*k212)-(x2(k)+h*k222);k41=x2(k)+h*k32;k42=-(xl(k)+h*k3l)-(x2(k)+h*k32);xl(k+l)=xl(k)+h*(k11+2*k21+2*k3l+k41)/6;x2(k+l)=x2(k)+h*(k12+2*k22+2*k32+k42)6;y_rk4_05(k+l)=xl(k+l);end%求出误差最大值err_0=max(abs(y_0-y_rk4_0);err_05=max(abs(y_05-y_rk4_05);%输出结果dispC最大误差(R=O)最大误差(R=O.5»err_max=err_0,err_05;disp(err_max)步长h取0.0001到0.5之间改变，并且将数值解干脆与解析解比较，求出最大误差数据列入下表中。步长方0.00010.00050.0010.0050.010.050.10.5最大误差5.4330XlO-IO1.6969XlO-IO1.0574XlO-IO4.1107×10-96.6029×10-84.1439×10-56.6602X10-44.2988XlO-I40.5最大误差2.7649XlO-Il6.8123×10-125.3753×10-124.0902XlO-IO6.5425×10-94.1365×10-66.7152×10-54.5976×10-2从上表中可以看出，当步长h=0.001时，总误差最小；当步长h小于0.001时，由于舍入误差变大而使总误差增加；当步长h大于0.001时，则由于截断误差的增加也使得总误差加大。另外，当系统的解改变激烈时.（如R=O）,误差对步长的改变较为敏感；当系统的解改变平稳时,步长的改变对误差的影响就要缓和得多。数值积分算法确定以后，在选择步长时，须要综合考虑。三、试验要求1 .参照已知程序，改动程序中的参数和输入量，验证输出结果。2 .分析不同参数取值对输出结果的影响。3 .试验结果写成试验报告六（全部试验完成后交试验报告）。