2023年中南大学自动化胡杨系统仿真实验报告完整版

中南大学系统仿真实验报告指导老师：实验者:学号:专业班级:完毕时间:实验一ＭATLAＢ中矩阵与多项式的基本运算基本命令训练：eye（m)取n＝3,程序如下：>>ｅyｅ(3）ans＝10001000１结论：eye(n)用于产生n×n维的单位矩阵，在这里n取3,故产生３×3维单位矩阵。one(n)、ones(m,n)对oｎes(n）取n=5,对ones(ｍ,n)取m=３,n=５，程序如下:>>ｏneｓ(5)ans=11１１１１1１１111１1111111１1111>>ones(３,5)anｓ＝111１11111111111结论：oｎeｓ(n)用于产生n×n维的全1矩阵，在这里ｎ取5,故产生5行5列全1矩阵。ones(m，ｎ)用于产生ｍ×n维的全１矩阵,在这里ｍ取3,n取5,故产生3行5列的全1矩阵。zｅrｏs(m，n)取m=３，n=2,程序如下:＞>zeros(３,2)ａns=00０0０0结论：zeros(m,n)用于产生m×ｎ维全0矩阵，在这里ｍ取3,n取2，故产生3行2列全0矩阵。4．ｒaｎd(m,ｎ）取m=3,n＝4，程序如下:>>rａｎｄ(3,4)ans=0．95０10.４8６00．４56５０．４4４70.23１10.89130．01850.61540.6０6８０.7621０.８214０.791９结论:rand（m,ｎ)用于产生ｍ×n维平均分布的随机矩阵，在这里m取3，ｎ取４,故产生了3行4列的随机矩阵5．diag(v)先创建３×3的魔方矩阵ｖ,在进行diag（v)运算,程序如下：>>ｖ＝maｇｉc(3)diag(v)v=816357492ans=8５2结论:diag(v）用于得到矩阵v的对角元素6．A\B、A/B、inv（Ａ)*B、B*inv(A)先创建A、B两个矩阵,在进行运算，程序如下:>>A=[1,2;３，4]；＞>Ｂ=［5,6;7，8]；>>a=A\Bb=A/Ｂｃ=ｉｎv（Ａ)＊Bd=Ｂ＊inv（A)ａ=-３-４45b=3．00０0－2．000０2.００00-1.００00c=-3.00０0－４.0０004.00００5.00０0ｄ=-1.000０２.0000－2.00０03.0０00结论:’/’表达矩阵右除,’\’表达矩阵左除，iｎv(Ａ）表达求A的逆矩阵，由实验结果可知，矩阵左除与右除结果不同样，矩阵左乘与右乘结果也不同样,A\Ｂ是求ＡＸ＝Ｂ的解，A／Ｂ是求XＢ=Ａ的解。所以编程求解的时候要注意区分他们的区别。7、roｏts（ｐ)>>syｍsx;>＞a=3*x.^3+2＊x+5;>>p=［3，0,２,5］>>roots(ｐ)p＝３０2５aｎs=0.5０00+1.1902i0.5000-1.1９0２ｉ-1.０000结论:ｒｏotｓ(p)函数用于求多项式的根,以向量形式输入多项式的系数,相应降幂排列，然后调用函数,即可求得相应多项式的根。8、ｐoｌｙ>>A＝[1,2；3,4]；＞>ｐoly（A)ans=１.００00-5.0000－２．000０结论：poｌy(A)用于求矩阵A的特性多项式的系数9.conv、dｅｃonv>>A=［1,2];>>B=[3,4］;>＞a=conｖ（Ａ,B)ｂ=ｄeconv(Ａ，B）a=3108b=0．3３33结论：使用conｖ函数对多项式进行乘法运算,其使用格式为c=conv(a,b),其中a和b为两个多项式的系数向量，c为相乘所生成的多项式的系数向量。使用ｄeｃonv(a,b)完毕除法运算。A*B与Ａ.*Ｂ的区别＞>A=[１,２];>>B＝[5,６]';>>ａ=A＊BA=[1,2］;B＝[5，6];b=A.*Ｂa=１7ｂ＝512结论:Ａ.*B称为“点乘”、“位乘“，即为两个行列数相同的矩阵，相应位置一一相乘，得到的结果依位置相应到结果矩阵中，而A*B为矩阵乘法,规定前者A的列数与后者B行数相应。1１.ｗhｏ与whｏs的使用>>A＝[1,２；3,4];>>ｗhowhｏsYoｕrvａriablｅsａre:ANameＳizeBytesClａssAttｒibuｔesA2x2３２ｄouｂle结论：who给出变量的名称清单;而whos给出所有变量的具体信息。ｄiｓｐ、ｓiｚe(ａ)、lｅngtｈ(a)的使用＞>a=＇helｌｏｗorlｄ＇;>>dｉsp（a)a=[1，2，3,4]；B=sｉｚe(ａ)C=length(a)heｌｌoworｌdB=14C＝4结论：disp函数的作用是直接将内容输出在Matlab命令窗口中,size(a)表达矩阵每个维度的长度，ｌength（ａ)表达矩阵a的最大的长度。实验二ＭAＴＬAＢ绘图命令基本命令训练１．ｐlot2.logloｇ3．sｅmiｌoｇx４.ｓeｍilogy５.ｐolar6.titlｅ7.xlａｂel8.yｌabel９．teｘt１0．gｒid11.ｂaｒ1２.ｓｔａiｒｓ13.ｃｏntoｕｒ1.＞＞t=[0:ｐｉ/３60：２*ｐｉ*22/3];x＝93*cｏｓ(t）+36*ｃoｓ(t＊４.15）;ｙ＝9３*sin(ｔ)＋36*ｓｉn(t*4.15);plot(y,x）,grid;实验结果为：>>t=[０:pi/3６0:２*pi*22/3];x=93*ｃｏs（t)+３6*cｏs(ｔ＊4.15);y=９3*siｎ（t）+３6＊sin(t＊4．15）;plot（y，ｘ）实验结果为:实验结论：plot（)用于绘制二维曲线，grid用于切换有无网格的状态。2.t=0:0.05:１00;ｘ=t;y=２＊t;z=siｎ（2*t）；pｌot３（x,y,z,＇b:')实验结果为：实验结论：ｐlot3(x,y,z)用于绘制三维曲线，ｂ表达设立曲线的颜色为蓝色,：表达曲线线型为点线，格式为ploｔ3(函数参数，函数参数，’曲线参数设立’)3．t=０:pi/20:２*pi；y=ｓｉn（ｘ);staｉrs(x,y)实验结果为：实验结论：staiｒs(x,y)表达绘制出的二维曲线为阶梯图。４．ｔh＝[pｉ/200:pi／200:2＊pi]';r=cos（2*th）；pｏlａr(ｔh,r),grid实验结果为：实验结论:pｏlar()用于绘制二维曲线的极坐标图。５．ｔh=[0:pi/1０:2*pi];x=ｅｘp（j＊ｔh);ｐlｏt(ｒeal(x),imaｇ(ｘ),＇ｒ*');grid；实验结果为：实验结论:r表达设立曲线颜色为红色，*表达曲线的数据点形为星号。６、>＞x＝0:１０0０;＞＞y=0:100０；>＞lｏglｏg(ｘ，y）;ｔｉtle('Loglｏg');ｇriｄon;实验结果为:实验结论：loｇloｇ（)用于绘制横纵轴均为对数刻度的图形,title（)用于为图形添加标题,本例为添加Logｌog作为标题。7、＞>x=0：1000;>>y=０:1000;>>semiloｇx(ｘ,y);titｌｅ('Lｏｇｌｏg');gridon;实验结果为：将semilogx换成semｉｌｏgy程序如下：>＞x=0:1000;>>y＝０:1０00;>>semiｌogy(x,y);title('Ｌogloｇ');griｄoｎ;实验结果为:实验结论：sｅmilogx（）用于绘制半对数图，其中x轴坐标为对数,若换成semilｏgy则表达y轴坐标为对数。8、>>x=０：1000;>>ｙ＝0：1０00;>>plｏｔ(x,ｙ）；＞>x=0：１000；>>y＝0：1０00;>>pｌｏt(x,y);grｉdoｎ;xlabｅl('\fontsize{20}\iｔx＼rｍ＇)；ｙlaｂeｌ（'\fｏｎtｓize{20}ｙ＇)；text(500,500,'中点')实验结果为:实验结论:xlａbｅｌ和ylａbｅｌ分别表达给x轴和y轴添加标注，ｔext(x,y,’ｓtrinｇ’）用于给图形坐标(ｘ，ｙ）处书写注释，本程序给x轴和y轴分别标注x，y，,在(50０，５00）坐标处注释“中点”。9、>>t=0:pi／1００:2*pｉ;＞＞alpha=3;>>y＝ｓin(aｌｐha*ｔ);＞>ｂar(ｔ,y);gridon;实验结果为：实验结论：bar(ｘ,y)用于绘制二维条形图。10、>>x=－8:0．5：8;>＞y=－8:0.5:8;>>［ｘx，yy]=mｅshgｒｉd（x,y)；＞>c=sqrt(xx.^２+ｙy.^2）+eps;>＞z=siｎ（c)./c;>>ｃontｏur(xｘ,yｙ,z)实验结果为：实验结论：ｃoｎtour(x,y,z)用于绘制等高线。补充实验:多窗口绘制图形suｂpｌｏt(）>>subploｔ(2,2,１)；t=[０:pi/2０0:2＊pｉ]；y=ｓin(t);ｐloｔ(t，y）;ｓubplot(2,2,２);t=［0:pi／200:2*ｐi];y=coｓ（t);ploｔ(t,ｙ）；ｓuｂplot(2,2,4);t=[0:pi／２00:2*ｐi];ｙ=t;plot（t,y);实验结果为:实验结论:本实验测试sｕｂplｏt（)函数,由实验结果可知，ｓubplot（）函数中某一个未编写并不会影响整个函数的运营，只是未编写的那个部分不显示,其他的照常显示，比如编写了subpｌｏt(2,２,1)，sｕbplot(2,2，２),ｓｕｂplot(２,2,4)，但是未编写subplot(2，2,3),那么结果只显示subpｌoｔ(2,2，1),ｓuｂplｏt(2,2,2),subpｌｏｔ(2,2,４)中的结果,并且顺序按原位置,而ｓubplｏt(2，２,３）的不会显示。实验三MATLAB程序设计１.计算1~1000之内的斐波那契亚数列＞＞ｆ=[１，1];＞>i=1；>>wｈilef(i）+f(i＋1)<１000f（i+2)=f(i）＋f（ｉ＋１)；i=ｉ＋１；end>>ｆ，ｉf=Ｃolｕmｎs1tｈｒoｕｇh1011235813２13455Cｏlumnｓ11thｒough１６8９１4４233377６１09８7i=152．>>m=３;>＞n＝4;＞>ｆoｒｉ＝１:mforj=１:na（i,j）＝1/(i+ｊ－1）;endeｎd>>fｏrmaｔraｔ>>aa=１1/21/３1/41/２1/３1/41/51/31／41/51/63.>>m=3；n=4;fori=１:mfoｒj=1：ｎa(ｉ,ｊ）=1/（ｉ＋ｊ-1);enｄｅndaa＝10．50.333３3０.250．50.3３3330.2５0.20.333330．２５０.20.16667实验２用了formatrat，结果为分数表达，实验3没有则用小数表达。4、>>x=iｎpｕt(＇请输入x的值:');iｆｘ<=１0;ｙ=coｓ(x+1)+sqrt(x*x＋1)；elseifx>15y＝x*ｓqｒt(x+sqrｔ(ｘ）);elsey＝x；endy请输入ｘ的值:10ｙ＝1０.054>>x＝iｎput（'请输入x的值：＇)；ifｘ＜=１0;y=ｃos(x+1)＋ｓqrｔ(x*x＋１);eｌseifx>15y=x＊sqrｔ(x+sqrt(x))；ｅｌsey=ｘ;ｅｎdy请输入ｘ的值:1１y=１15、去掉多项式或数列开头的零项.＞>ｐ=[00０1３０2009];fori＝1:leｎgｔh(ｐ),iｆp(1)==0,p＝p(2:length(ｐ))；ｅnd;end;ｐp=1302００96、建立ＭＡTLAB的函数文献，程序代码如下,以文献名eｘ２_4.m存盘点击File-New-Ｆｕncｔion建立文献,文献名为ex2＿４．ｍ,结果如下:在ＭＡTLAB的命令窗口输入ex２_4(200)，得到运营结果:>＞eｘ2_４(200)anｓ=112３58１３２13455８９144在ＭATLAB的命令窗口输入lookforffibno,得到结果：>>loｏｋｆorfｆiｂnoex2_４－ffibno计算斐波那契亚数列的函数文献在ＭATLAB的命令窗口输入helｐex2_４，得到结果：＞＞ｈeｌｐeｘ２_4ｆfiｂno计算斐波那契亚数列的函数文献ｎ可取任意自然数程序如下程序设计题functiｏnsuｓhun=inｐut(＇请输入一个数n：＇）；if（n==1)fpｒｉｎtf('１既不是素数也不是合数\n'）；ｄiｓp(＇是否继续？');disｐ('1.是；2．否')；b=input('')；ifb＝=１sｕsｈｕ；eｌｓeｄisp('谢谢使用！');ｂｒeak;endendif(n=＝２)fpｒiｎtf（'2是素数\n＇)；ｄｉｓp(＇是否继续？'）；dｉsp('１.是；2.否＇)；b=iｎput('');ifｂ==１sｕｓhu;elｓｅdisp(＇谢谢使用！');breａｋ;ｅｎdenｄiｆ(n＝=3)fｐrｉｎtf（'3是素数\ｎ'）;disｐ（＇是否继续？');dｉｓp('1.是;2.否')；b=input(＇');ifb=＝1suｓhｕ;eｌsedｉsp('谢谢使用!'）;brｅak;enｄeｎdif（ｎ>３)fori=２:(n-１）ifmod(ｎ，ｉ)＝=0ａ＝ｎ/i；t=0;ｆpｒinｔf('%d＝%d*%ｄ\n'，n,a,ｉ);elｓeｔ=1；eｎｄｅnｄｉf（ｔ==1)ｃｌearｔ;ｆprintf（'%d不是素数\n',n);ｅlseｆprintｆ('％d是素数\n＇,ｎ);cleart；enddisp('是否继续？');disp（'１.是；２．否');ｂ=input（''）；ifｂ==１susｈu;elsedｉsｐ('谢谢使用!');eｎdend运营结果为：>＞sushu请输入一个数n:11既不是素数也不是合数是否继续？1．是;2.否1请输入一个数ｎ:２２是素数是否继续?1.是;2．否１请输入一个数ｎ:3８３8＝１9*2３8=２＊193８不是素数是否继续？1.是;2．否１请输入一个数n:2323不是素数是否继续？1.是；2.否2谢谢使用！实验四ＭATLAB的符号计算与SIＭULIＮK的使用程序举例求矩阵相应的行列式和特性根>>a＝syｍ(＇[１00１]');>>ｄa=deｔ(ａ)ea=eig（a)ｄa＝1ea=１实验结论:det(）函数用于计算矩阵对于的行列式的值,ｅig()函数用于计算矩阵的特性值和特性向量２.求方程的解（涉及精确解和一定精度的解)＞>r１=solve('x^２-x-1')rv=ｖｐａ(r1)ｒv4=vpa(r1,４）rv20＝vpa(r1,20)r1＝1/2－5^(1／2)/2５^(1/2)/２+1/２rv=-0.１.６83436５6ｒv4＝－0.６181．６18rv20=-0.6８9484821．６８9484８２实验结论：vpa（s,n)称为变精度算法函数，表达将ｓ表达为n位有效数的符号对象3、>>a＝sym('ａ');ｂ=sym（＇ｂ');c=syｍ('ｃ'）;ｄ=ｓyｍ('ｄ＇);%定义4个符号变量ｗ＝10；x=5;y=-８;z=11;%定义4个数值变量Ａ=[a,b；c，d]%建立符号矩阵AB=[w,x;y,z]％建立数值矩阵Bｄet(Ａ)%计算符号矩阵Ａ的行列式det（Ｂ)%计算数值矩阵B的行列式A=[ａ,b］［c,d]B=105-811ans＝ａ＊d－b*ｃanｓ＝15０4、>>sｙmsxy;ｓ=（-7*x^2-8*ｙ^2)*(-x^2+３＊y＾2);expanｄ(s)%对s展开colleｃt(s，x)％对s按变量x合并同类项(无同类项)factor(ａns)%对ans分解因式ａnｓ＝7*x^4-１3＊x^2*y＾2-2４*y^4ａnｓ＝7*x^4+（-１3*y^2)*x^2-２4*y^4ans=(7*ｘ^２＋8*y^2)*(x＾２－3＊y^2)实验结论:expaｎd函数用于多项式的展开运算,collect函数用于符号表达式的展开运算和合并同类项，factoｒ用于对函数进行因式分解。5、对方程AX=b求解>>A=[34，８,4;３，34,3；3，６,8];b＝[4;6;2］；X＝lｉｎsｏlve(A,b)%调用liｎsolve函数求解A＼ｂ%用另一种方法求解Ｘ=0.0674820.16１370．1036７ａｎs=0.0６7482０.161３70.1036７实验结论:运用linsoｖｅ（Ａ,b)与A＼b结果同样,都用于对ＡX＝b进行求解6、对方程组求解>＞ａ11＊x1+ａ１2*x2＋a13＊x3=b1a21*x1+a２2*x2+a2３＊ｘ3=b2a31＊x１＋a32*x2+a33＊x3=ｂ3A=［a11,a12,a１３;a２1,ａ２2,ａ23；a３1,a32,a33]；b=[b1；b2;ｂ3］;ＸX=A\bXＸ=(a１2*ａ2３*ｂ3－ａ１2＊b2*a33+ａ13*a32*b2-a1３*ａ2２＊b3＋b1*ａ22＊a33-b1＊ａ32＊a23)／(a11*ａ22*a33-a１1*a３２*ａ2３－a21*a12*a33＋a3２*ａ2１*a13-a22＊a３1*a13+ａ３1*a12＊ａ23)-(a１1*a2３*b３-a11*b2＊a３3-ａ２1*a13＊ｂ３-a23*a３１*ｂ1+b2＊ａ31*ａ1３+a２1*b1*ａ3３)／(ａ11*a2２*ａ33-a１１*a32＊a2３-a２1*a１2*a3３＋ａ３2*a21＊a13－a2２*a31*a1３+a31*a12＊a23）(ａ3２*a２1＊b１-a11*a32*b2+ａ11*a22*b３-a2２＊a3１＊b1-a21*a１2*b３+ａ３1*a1２＊ｂ2)/(a11*a２２*a33－a11*ａ3２＊ａ２３－ａ2１*ａ12*a3３+a32＊a21＊a13-a22＊a３１＊a1３+a31*a12＊a23)7、>＞ｓymsａbｔxyz;f=sqrt(1+ｅｘｐ（ｘ)）;diff（f)％未指定求导变量和阶数,按缺省规则解决f=x*cｏs（ｘ);diff(f,x,2)%求f对x的二阶导数diｆf(f,x,3）%求f对ｘ的三阶导数f1＝ａ*coｓ(t）;f2=b*sｉｎ(ｔ)；diff(ｆ2)/ｄiｆf（f1）％按参数方程求导公式求ｙ对x的导数ans=exｐ(x)/(2＊(ｅxｐ（ｘ)＋１)＾(1/2)）ans=-2＊sin(x)－x＊ｃoｓ(x）aｎｓ=ｘ＊ｓｉn(x）-3*coｓ(x）aｎｓ=－(b＊coｓ（ｔ)）/(ａ*sin(ｔ))SIMULINK的使用选择合适子模块构造控制系统如下：选择Sｉmｕlatiｏn菜单下的ｓtart命令运营仿真,在示波器(Scope）中观测结果：R（s)为阶跃输入,C(s)为输出实验结果为:实验五ＭAＴLＡB在控制系统分析中的应用基本命令1．ｓtep2.impulｓe3.iniｔial4．lsim5.rlocfind6．bｏｄe７.mａrgin8.nyqｕｉst9.Nichｏｌs1０.clｏｏp1、求下面系统的单位阶跃响应＞>ｎｕm＝［４]；ｄｅｎ＝[1,1,4];sｔｅｐ(ｎum,den)[y,x,t］=steｐ(ｎum,den);tｐ=spline（y,t,mａx(y))max(y)ｔp=1.6062ａns=1.４４41实验结论:ｓtep(nuｍ,ｄeｎ)用于绘制系统阶跃响应曲线,spｌine(ｙ,t,maｘ(y）)函数是由y,ｔ的值计算mａx（y）相应的函数值t，在本例中即峰值时间，而max(y)用于计算峰值。2、求如下系统的单位阶跃响应>>a＝［0,1；-6,－5];b＝[0;1];ｃ＝［１，0]；d＝0；[y,x］=step(ａ,b,c，d）;pｌｏt(y)程序结果为：实验结论：steｐ(a,b,c,d）用于绘制系统阶跃响应曲线３、求下面系统的单位脉冲响应%程序如下：>>num＝[4］;ｄen=[1,１,４];impulse(nuｍ,den)实验结果:实验结论:iｍｐulse(nuｍ,ｄeｎ)函数用于求取系统单位脉冲响应，其用法基本同step函数。4、已知二阶系统的状态方程为：ﻩ求系统的零输入响应和脉冲响应。%程序如下:>>a=［０,1;-１０,-2]；b=[0;1］;c=[１,0］；d=[0];ｘ０=［1,0];ｓｕｂplot(1，2,1);ｉniｔｉal(a,ｂ,ｃ,d，x0)sｕｂpｌot（1,2,2）;impulse(ａ,b,c,ｄ)实验结果为:实验结论:initial（a,b，ｃ,d,x０）用于求解零输入响应系统,x０为初始条件,impuｌｓe（a，b,c，d)用于求单位脉冲响应。5、系统传递函数为:输入正弦信号时,观测输出信号的相位差。%程序如下:>>nｕm=［１];ｄｅn=[1,1］；t＝0:0．01:1０;ｕ=sin（2*ｔ)；holdonplot(t，u,'r＇)lsｉｍ（num,ｄｅn,u,t)程序结果为：实验结论：lsiｍ(num,dｅn,u,t)用于求系统对任意输入ｕ的响应。6、有一二阶系统,求出周期为4秒的方波的输出响应％程序如下：num=[2５1]；den=[12３];t=(0:.１：10)；ｐeriod=4;ｕ＝(rｅm（t,pｅrｉod)＞=pｅriod./2);％看rem函数功能lsｉm(num,dｅn,u，ｔ);实验结果为:7.已知开环系统传递函数,绘制系统的根轨迹，并分析其稳定性％程序如下：num=[１２]；ｄen1=［143];den=cｏnv（den1,deｎ１);fiｇｕrｅ(1)rlocus(ｎum,deｎ)［k,p]=rloｃfｉｎd(nｕｍ,den)ｆigｕre(2)k＝５５;num１＝k*[12];deｎ=［１43];deｎ1=cｏnｖ(deｎ,den）;[num,deｎ］=cloop(nｕm1,ｄen1,-1）;impulｓe(ｎｕm,ｄen）titlｅ('impulsereｓpoｎse（k=５5)')figurｅ（３)ｋ＝56;num1＝k＊［１2];den=[143];den1=conv(ｄen,den);[nｕm,den]=cｌoop(num1,dｅn1,-1）;imｐulｓe(ｎum，ｄeｎ)tｉtｌe('impulｓｅresponse(ｋ=５6)'）实验结果：Ｓelecｔaｐointinthegraｐhiｃｓwinｄｏwselecｔed_point=-3.３886+2.3６02ik=２3．5611p=-５.0868-0.43５５＋２.28３１ｉ-０．435５-２.28３１i-2．０423实验结论：dｅｎ＝conv(A,B)用于多项式A,Ｂ以系数行向量表达，进行相乘,rｌｏcus(nｕm,deｎ）用于绘制指定系统的根轨迹。分析得系统是不稳定的。函数rlocｆinｄ用于计算给定一组根的根轨迹增益8、作如下系统的bode图%程序如下：n=[１,1］;d=[1,4,１１,7］;boｄe(n，d)实验结果为:实验结论：bｏde（n,d）用于绘制Bode图9.系统传递函数如下求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表达的系统频率响应%程序如下：nuｍ=[1];dｅn=conv(［1２],conｖ([12］,［１2])）;ｗ＝logｓpace(-1，２);t=0.5；[m1,p1］=ｂode（num,den,2);p1＝ｐ１-t*w'*18０/pｉ；[ｎ2，d2]=pａde（t,4);ｎuｍｔ=conｖ(n２,num);dent=（conｖ（den,d2））;[m２,p2］＝bｏｄe(numt，dent，w);subｐｌot(2,1,1);semｉloｇx(w，20*log10(m１），w,2０*log１0（ｍ2）,＇g－-')；gridon；ｔitle('ｂodeｐlot');xlabel(＇fｒeｑuency＇);yｌabel('gaｉｎ'）;ｓuｂｐloｔ(２,１,2);seｍilogx(w,p1，w,p2,'g--')；griｄon;ｘlａbel('ｆrequenｃy')；ylａbel（'phasｅ');实验结果为：10.已知系统模型为求它的幅值裕度和相角裕度%程序如下:n=[３．５]；d＝[１232];［Gm,Ｐm，Wcg,Wcｐ]＝ｍargin(n,ｄ)实验结果为:Gm=1.1433Ｐm=7.１６88Wcg＝1.7３23Ｗｃp=1．6５41实验结论:[Gm,Ｐm,Wcg,Wｃp]=margiｎ(n,d)给出系统相对稳定参数,返回参数分别为幅值裕度、相角裕度、幅值穿越频率、相角穿越频率。１1、二阶系统为:令ｗn=1，分别作出ξ=2,1，0.7０７,０．5时的ｎyquｉｓt曲线。%程序如下:n=[１];ｄ1=[１,4，1];d２=［１，2,1];ｄ3=[1,1．４14，１];d４=［1,1，１]；nyquiｓｔ(n，d1)；holdonnyquisｔ（n,ｄ２);ｎｙquist(ｎ,d３);ｎyquist(n,d4)；实验结果为:实验结论:nyqｕiｓt(sys,w)函数用于绘制系统Ｎｙｑｕist图,由用户指定选取频率范围。１２、已知系统的开环传递函数为绘制系统的Ｎyqusｉt图，并讨论系统的稳定性%程序如下:G=tf(100０,coｎv(［1,3,2],[１，５］));nyquisｔ(G);aｘｉs('sqｕare')实验结果为:Ｗａrniｎｇ:ThisplottypeｄoesnotsｕppｏrｔthisｏptionfortheAXISｃｏｍmａnd.＞IｎD:＼MATLＡＢ6p5p1\tooｌboｘ＼ｃontroｌ＼ｃｔrlｇｕis\@ctrluiｓ\@axesｇroup\adｄbypass．matｌｉnｅ1１4InD:\MATLＡB６ｐ5p1\ｔooｌbｏx\maｔlaｂ\graph2d\ｐｒivate\mwbyｐａss.patｌiｎe24ＩnD:＼MＡTLAB6p5ｐ1＼ｔoｏlｂｏｘ\mａｔｌab\graph2ｄ＼aｘｉs．matline7613.分别由w的自动变量和人工变量作下列系统的ｎyquｉsｔ曲线:%程序如下：ｎ＝［1];d=[1,1,0]；nyquist(n,d）;％自动变量n=[1]；d＝[１,1，0];ｗ=[0．5:0.1：3];nyquist(n,d,w);%人工变量实验结果为：实验结论:nｙquisｔ(n，ｄ，ｗ)用于求指定范围ｗ的奈氏值14、一多环系统,其结构图如下,使用Ｎyｑuisｔ频率曲线判断系统的稳定性。%程序如下：ｋ1=１６.7/0.０12５;z１=[0];ﻫp1=［－1.25-４-16］;[ｎum1,ｄｅn１]＝zp２tf(ｚ1，p1,k１）；

[ｎum,dｅｎ]=cloｏp(ｎｕｍ１,den1);ﻫ[z，p,k］＝tf2zp(ｎuｍ，dｅn);p

figure(1)ﻫｎyqｕiｓt(ｎum，den)

fｉgure（2）

[nｕm2,deｎ2]=clooｐ(num，den）；ﻫimｐulse(num２，den２）；实验结果为：p=-1０.5９７+36.２１5i-１0．5９７-36.2１5i－0．05６１87实验结论：[num1,ｄen1]=zp2tf（z1,p1,ｋ1）用于将系统函数的零极点转化为系统函数一般形式的系数,ｃlｏop(）函数用于将系统通过正负反馈连接成闭环系统,分析得系统稳定15.已知系统为:作该系统的nicｈolｓ曲线。%程序如下:n=[１]；d=[1,1,0];ngｒid(‘new’）；ｎiｃｈｏls(n，d);16、已知系统的开环传递函数为:当k=2时,分别作nｉchols曲线和波特图。%程序如下：ｎum=1;dｅｎ=conv(conｖ(［10],[１1］）,[0．51]);subｐloｔ(1，2,1);nichoｌs（ｎum,ｄen)；gｒid;％niｃhｏls曲线ｓubplot（１,2,2);g=tf(nuｍ，dｅn）；ｂode(feedbａｃｋ(g,1，-1));griｄ；％波特图17.系统的开环传递函数为:分别拟定ｋ=2和k=1０时闭环系统的稳定性。％程序如下：d1＝[1，3，2,0]；n１＝[2];［nc1,dc1]=clooｐ(n1，d1,-1);rootｓ(dc1)d2=d１；n2=[1０];[nc2,dc2]=cｌooｐ（n２,d2，-1）；ｒｏots(dc2)实验结果为:ans=-２．5214－0.23931+0．85７87ｉ-0.２３931-0.85787ｉans=－３.308９0．1544５＋1.7316i0.1５445-1．731６i1８、系统的状态方程为:试拟定系统的稳定性。%程序如下：a=［-4,-3,０；１,0，0;0，１,0];b=[1;0；0］;ｃ=[0,1,２];d＝0;eig(ａ)%求特性根raｎk（ctｒｂ(ａ,ｂ))实验结果为:anｓ=0－１-3anｓ=3系统可控性辨别矩阵的秩为3，满秩，系统可通过状态反馈配置极点使系统稳定。实验六连续系统数字仿真的基本算法程序举例１.取h=0.２,试分别用欧拉法、RK２法和RＫ4法求解微分方程的数值解，并比较计算精度。注:解析解:％程序如下clｅａｒt(１)=０；%置自变量初值y(1）＝1;y_eulｅr(1）=1;ｙ_rk2(１)=1;y_ｒｋ4（1)=1;%置解析解和数值解的初值h=0.2;%步长%求解析解ｆorｋ=1:５ｔ(k+1)＝ｔ（k）＋ｈ；ｙ(k＋1)=sｑrt（1+2*t(k+1)）;enｄ%运用欧拉法求解forｋ=１:5ｙ_eulｅｒ（k+1)=ｙ_eｕlｅr(k）+h*(y_eｕler（k）-2＊t(k)/ｙ＿euｌer(k)）;end%运用RK2法求解fｏrk=１：5k1＝ｙ＿rｋ2(k）-2＊t（k）／y_rk2(k);ｋ2=(ｙ_rk2(k)+h*k1)-2*(t(k)+h）／(ｙ_ｒk2（k)+h*k1）;y＿ｒk2(k+1)=y_ｒk2(k)+h＊(ｋ1+ｋ2)/2;end%运用RK4法求解fork=１:5ｋ1=ｙ_rk4(ｋ)-2*t(k)/y＿ｒk4（k);k2＝(y_rｋ4(k)＋h*k１/2)-2*(t(k)+h/2）/（ｙ_ｒk4(k）+h*k1／2);k3=(ｙ_rk4(k）＋ｈ*k2／2)－2*(t(k)＋h/2）/(y＿rk4(ｋ）+h＊k2/2);k4=(y＿rk4（k)＋ｈ*k3)－２*(t（k)＋h)／(y_rk4(k）+ｈ*k3)；y_rk4(k+1)=ｙ_ｒk4(k)+h＊（k１＋2*k2+２*k3+k4)／6;end％输出结果disp('时间解析解欧拉法RK2法RK4法')yt=［t',y＇,y_eｕｌeｒ',y_ｒk2',y_ｒｋ4＇］；disp(yｔ)程序运营结果如下:时间解析解欧拉法RK2法ＲK4法０１．000０1.0０0０１．０00０1.０00０0．２０23１．18321.20２31.186７1.１８320.40001．34１61.3７331.３483１.３４１7０.60001.48321.531５1.４９371.4833０．8000１.61251.6８111.6２７9１.61２51.００001.73211.82６91.７5421．７3212、考虑如下二阶系统:在上的数字仿真解（已知:，),并将不同步长下的仿真结果与解析解进行精