系統(tǒng)仿真技術(shù)實(shí)驗(yàn)報(bào)告

1、系統(tǒng)仿真技術(shù)實(shí)驗(yàn)報(bào)告學(xué) 院：信息科學(xué)與工程學(xué)院專 業(yè)：電氣工程及其自動(dòng)化班 級(jí)：XXXXX 姓 名：XXXXX 學(xué) 號(hào)：XXXXXXXXXXX 目 錄一、 實(shí)驗(yàn)一（MATLAB中矩陣與多項(xiàng)式的基本運(yùn)算）.41.1 實(shí)驗(yàn)題目.41.2 實(shí)驗(yàn)擴(kuò)展.131.3 實(shí)驗(yàn)心得.20二、 實(shí)驗(yàn)二（MATLAB繪圖命令）.212.1 實(shí)驗(yàn)題目.212.2 實(shí)驗(yàn)擴(kuò)展.262.3 實(shí)驗(yàn)心得.29三、 實(shí)驗(yàn)三（MATLAB程序設(shè)計(jì)）.303.1 實(shí)驗(yàn)題目.303.2 實(shí)驗(yàn)擴(kuò)展.373.3 實(shí)驗(yàn)心得.38四、 實(shí)驗(yàn)四（MATLAB的符號(hào)計(jì)算與SIMULINK的使用）.394.1 實(shí)驗(yàn)題目.394.2 實(shí)驗(yàn)擴(kuò)展.484

2、.3 實(shí)驗(yàn)心得.49五、實(shí)驗(yàn)五（MATLAB在控制系統(tǒng)分析中的應(yīng)用）.505.1 實(shí)驗(yàn)題目.505.2 實(shí)驗(yàn)心得.71六、 實(shí)驗(yàn)六（連續(xù)系統(tǒng)數(shù)字仿真的基本算法）.726.1實(shí)驗(yàn)題目.72七、總結(jié).78 實(shí)驗(yàn)一：MATLAB中矩陣與多項(xiàng)式的基本運(yùn)算1.1 實(shí)驗(yàn)題目1.1.1 實(shí)驗(yàn)任務(wù)： 1了解MATLAB命令窗口和程序文件的調(diào)用。2熟悉如下MATLAB的基本運(yùn)算： 矩陣的產(chǎn)生、數(shù)據(jù)的輸入、相關(guān)元素的顯示； 矩陣的加法、乘法、左除、右除； 特殊矩陣：單位矩陣、“1”矩陣、“0”矩陣、對(duì)角陣、隨機(jī)矩陣的產(chǎn)生和運(yùn)算； 多項(xiàng)式的運(yùn)算：多項(xiàng)式求根、多項(xiàng)式之間的乘除。1.1.2 基本命令訓(xùn)練：1eye(m)

3、2one(n)、ones(m,n) 3zeros(m,n)實(shí)驗(yàn)結(jié)果：4rand(m,n) 5diag(v)分析：rand(m,n):產(chǎn)生m*n的隨機(jī)陣diag（v）：以v的元素作為對(duì)角陣的元素產(chǎn)生對(duì)角陣實(shí)驗(yàn)結(jié)果：6AB 、A/B、 inv(A)*B 、B*inv(A)分析：X=AB是AX=B的解，且AB=A(-1)BX=A/B是XB=A的解，且A/B=AB(-1)Inv(A)*B:A的逆矩陣與左乘BB*inv(A)：B左乘A的逆矩陣實(shí)驗(yàn)結(jié)果：7roots(p)8Poly分析：roots(p):多項(xiàng)式p的根poly(C):用C的元素為根構(gòu)造多項(xiàng)式實(shí)驗(yàn)結(jié)果：9conv 、deconv分析：conv

4、:進(jìn)行多項(xiàng)式的乘法deconv:進(jìn)行多項(xiàng)式的除法實(shí)驗(yàn)結(jié)果：10A*B 與 A.*B的區(qū)別分析：A*B：進(jìn)行矩陣運(yùn)算A.*B：進(jìn)行數(shù)組運(yùn)算實(shí)驗(yàn)結(jié)果：11who與whos的使用分析：who:列出當(dāng)前工作空間中所有變量名字whos:列出當(dāng)前工作空間中所有變量，以及它們的名字、尺寸、所占字節(jié)數(shù)、屬性等信息。實(shí)驗(yàn)結(jié)果：12disp、size(a)、length(a)的使用分析：disp:屏幕輸出函數(shù)size(a):求矩陣的大小(獲取數(shù)組的行數(shù)和列數(shù))length(a):數(shù)組長度(即行數(shù)或列數(shù)中較大者)實(shí)驗(yàn)結(jié)果：12 實(shí)驗(yàn)擴(kuò)展1.2.1 向量及其運(yùn)算1）、向量產(chǎn)生實(shí)驗(yàn)結(jié)果：2）、向量基本運(yùn)算分析：dot:

5、維數(shù)相同的兩個(gè)向量的點(diǎn)積，結(jié)果是一個(gè)標(biāo)量cross:兩個(gè)向量的叉積就是一個(gè)過兩個(gè)相交向量的交點(diǎn)且垂直于兩個(gè)向量所在平面的向量，結(jié)果是一個(gè)矢量實(shí)驗(yàn)結(jié)果：1.2.2 矩陣及其運(yùn)算1）、矩陣產(chǎn)生通過冒號(hào)可以截取指定矩陣中的部分2）、矩陣基本運(yùn)算det(A):返回矩陣A的行列式的值eig(A):返回矩陣A的特征值inv(A):求A的逆矩陣fliplr(A):矩陣A關(guān)于中軸的反轉(zhuǎn)矩陣(左右反轉(zhuǎn))flipud(A):上下反轉(zhuǎn)rot90(A):矩陣A逆時(shí)針旋轉(zhuǎn)90度diag(A):抽取主對(duì)角線元素構(gòu)成向量rem(A,x):A矩陣除以模數(shù)x后的余數(shù)1.2.3 數(shù)組及其運(yùn)算實(shí)驗(yàn)結(jié)果：1.2.4 多項(xiàng)式及其運(yùn)算分

6、析：polyval(P,X):計(jì)算多項(xiàng)式P在X點(diǎn)的值（若X是向量或矩陣，則采用數(shù)組運(yùn)算（點(diǎn)運(yùn)算）polyvalm(P,X):以方陣X為自變量，計(jì)算多項(xiàng)式的值，采用矩陣運(yùn)算1.3、 實(shí)驗(yàn)心得MATLAB中數(shù)值運(yùn)算的部分較為繁多，但和之前學(xué)過的C語言等也有相同的部分。但MATLAB的語言更為精簡，比較方便。MATLAB是一款工程性質(zhì)十分強(qiáng)的軟件，只有大量的應(yīng)用實(shí)踐，才能夠?qū)ζ渲械暮瘮?shù)運(yùn)用自如，才能理解其中的內(nèi)涵。實(shí)驗(yàn)二 MATLAB繪圖命令2.1 實(shí)驗(yàn)題目21.1實(shí)驗(yàn)任務(wù)熟悉MATLAB基本繪圖命令，掌握如下繪圖方法：1坐標(biāo)系的選擇、圖形的繪制；2圖形注解（題目、標(biāo)號(hào)、說明、分格線）的加入；3圖形

7、線型、符號(hào)、顏色的選取。 基本命令訓(xùn)練1plot 2loglog 3semilogx 4semilogy5polar 6title 7xlabel 8ylabel9text 10grid 11bar 12stairs13contour實(shí)驗(yàn)結(jié)果：2.2 實(shí)驗(yàn)擴(kuò)展 2.2.1 圖形的修飾 1. 曲線顏色: b(blue),c(cyan),g(green),k(black),m(magenta),r(red),w(white),y(yellow)2.曲線線型：-（實(shí)線，默認(rèn)）， ：（點(diǎn)線）， -.(點(diǎn)畫線), -(虛線) 3.數(shù)據(jù)點(diǎn)形：.(實(shí)點(diǎn)), +(十字形), 。(圓圈) , *(星號(hào)), x(

8、叉號(hào)), s(正方形), d(棱形), h(六角形), p(五角形), >(右三角), <(左三角) 實(shí)驗(yàn)舉例：2.2.2 三維圖形的繪制仿真舉例：t=0:pi/30:10*pi; %設(shè)定時(shí)間范圍plot3(cos(t),sin(t),t,'-b','Linewidth',4)%繪制三維圖形，并作修飾grid on %加網(wǎng)格axis square %坐標(biāo)為方形figure(2) %新建圖形窗口2stem3(cos(t),sin(t),t,'-.k') %繪制三維火柴桿圖實(shí)驗(yàn)結(jié)果：2.3 實(shí)驗(yàn)心得用戶應(yīng)用MATLAB，只需指定繪圖方式，

9、并提供充足的繪圖數(shù)據(jù)，可以得到所需的圖形，也可對(duì)圖形加以修飾。二維圖形的繪制基本函數(shù)為plot。較為復(fù)雜的三維圖形分為三維曲線，三維網(wǎng)線圖和三維曲面圖。對(duì)符號(hào)函數(shù)或是無法求出函數(shù)參數(shù)之間顯示關(guān)系的隱函數(shù)，可以使用符號(hào)函數(shù)繪制圖形的函數(shù)。實(shí)驗(yàn)三 MATLAB程序設(shè)計(jì)3.1 實(shí)驗(yàn)題目實(shí)驗(yàn)任務(wù)1熟悉MATLAB程序設(shè)計(jì)的方法和思路；2掌握循環(huán)、分支語句的編寫，學(xué)會(huì)使用look for、help命令。31.2 程序訓(xùn)練1計(jì)算11000之內(nèi)的斐波那契亞數(shù)列 f=1,1;i=1;while f(i)+f(i+1)<1000 f(i+2)=f(i)+f(i+1); i=i+1;endf,i實(shí)驗(yàn)結(jié)果：2

10、 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endendformat rata實(shí)驗(yàn)結(jié)果：3 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1);endenda實(shí)驗(yàn)結(jié)果：程序2和程序3的區(qū)別：format rat 小數(shù)分?jǐn)?shù)表示。4 x=input('請(qǐng)輸入x的值:'); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y實(shí)驗(yàn)結(jié)果：5去掉多項(xiàng)式或數(shù)列開頭的零項(xiàng)p=0 0 0 1 3 0 2 0 0 9;for i=1:le

11、ngth(p),if p(1)=0,p=p(2:length(p); end;end;p實(shí)驗(yàn)結(jié)果：6 建立MATLAB的函數(shù)文件，程序代碼如下，以文件名ex2_4.m存盤function f=ffibno(n)%ffibno 計(jì)算斐波那契亞數(shù)列的函數(shù)文件%n可取任意自然數(shù)%程序如下f=1,1;i=1;while f(i)+f(i+1)<n f(i+2)=f(i)+f(i+1); i=i+1;end輸入完畢后在MATLAB的命令窗口輸入ex2_4(200)，得到運(yùn)行結(jié)果。在MATLAB的命令窗口輸入lookfor ffibno，得到結(jié)果：ex2_4.m: %ffibno 計(jì)算斐波那契亞數(shù)列

12、的函數(shù)文件在MATLAB的命令窗口輸入help ex2_4，得到結(jié)果： ffibno 計(jì)算斐波那契亞數(shù)列的函數(shù)文件 n可取任意自然數(shù) 程序如下實(shí)驗(yàn)結(jié)果：3.1.3 程序設(shè)計(jì)用一個(gè)MATLAB語言編寫一個(gè)程序：輸入一個(gè)自然數(shù)，判斷它是否是素?cái)?shù)，如果是，輸出“It is one prime”，如果不是，輸出“It is not one prime.”。要求通過調(diào)用子函數(shù)實(shí)現(xiàn)。最好能具有如下功能：設(shè)計(jì)較好的人機(jī)對(duì)話界面，程序中含有提示性的輸入輸出語句。能實(shí)現(xiàn)循環(huán)操作，由操作者輸入相關(guān)命令來控制是否繼續(xù)進(jìn)行素?cái)?shù)的判斷。如果操作者希望停止這種判斷，則可以退出程序。如果所輸入的自然數(shù)是一個(gè)合數(shù)，除了給出其

13、不是素?cái)?shù)的結(jié)論外，還應(yīng)給出至少一種其因數(shù)分解形式。例：輸入 6， 因?yàn)?不是素?cái)?shù)。則程序中除了有“It is not one prime”的結(jié)論外，還應(yīng)有：“6=2*3”的說明。實(shí)驗(yàn)結(jié)果：3.2 實(shí)驗(yàn)擴(kuò)展3.2.1 選擇結(jié)構(gòu)舉例：g=input('please input grade:','s'); %由用戶交互輸入成績等級(jí)switch(g) %判斷用戶輸入case'A','a'disp('85100');case'B','b'disp('7084');case

14、9;C','c'disp('6069');case'D','d'disp('<60');otherwisedisp('輸入錯(cuò)誤!')end實(shí)驗(yàn)結(jié)果：3.2.2 循環(huán)結(jié)構(gòu)舉例：A=magic(3)sum=zeros(3,1);for n=A %給定循環(huán)表達(dá)式，循環(huán)次數(shù)為A的列數(shù)，第一次取出A的第一列賦給n,依此類推共三次sum=sum+n;endsum' %輸出結(jié)果，以行向量的形式給出實(shí)驗(yàn)結(jié)果：3.3 實(shí)驗(yàn)心得MATLAB程序可以再任何文本編輯器中編輯完成，但使用MATLAB自身

15、的編輯器是個(gè)不錯(cuò)的選擇，因?yàn)樗嬗姓{(diào)試功能。MATLAB程序的流程結(jié)構(gòu)有條件轉(zhuǎn)移、循環(huán)、開關(guān)結(jié)構(gòu)以及試探結(jié)構(gòu)。MATLAB函數(shù)的編寫調(diào)試要求盡量按規(guī)范編寫程序，以便調(diào)試和維護(hù)，增強(qiáng)可用性。實(shí)驗(yàn)四 MATLAB的符號(hào)計(jì)算與SIMULINK的使用4.1 實(shí)驗(yàn)題目 4.1.1 實(shí)驗(yàn)任務(wù)1掌握MATLAB符號(hào)計(jì)算的特點(diǎn)和常用基本命令；2掌握SIMULINK的使用。4.1.2 程序舉例1求矩陣對(duì)應(yīng)的行列式和特征根a=sym('a11 a12;a21 a22');da=det(a)ea=eig(a)實(shí)驗(yàn)結(jié)果：2求方程的解（包括精確解和一定精度的解）r1=solve('x2-x-1&

16、#39;)rv=vpa(r1)rv4=vpa(r1,4)rv20=vpa(r1,20)實(shí)驗(yàn)結(jié)果：3 a=sym('a');b=sym('b');c=sym('c');d=sym('d'); %定義4個(gè)符號(hào)變量w=10;x=5;y=-8;z=11; %定義4個(gè)數(shù)值變量A=a,b;c,d %建立符號(hào)矩陣AB=w,x;y,z %建立數(shù)值矩陣Bdet(A) %計(jì)算符號(hào)矩陣A的行列式det(B) %計(jì)算數(shù)值矩陣B的行列式實(shí)驗(yàn)結(jié)果：4 syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %對(duì)s展開coll

17、ect(s,x) %對(duì)s按變量x合并同類項(xiàng)(無同類項(xiàng))factor(ans) % 對(duì)ans分解因式實(shí)驗(yàn)結(jié)果：5 對(duì)方程 AX=b求解A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=linsolve(A,b) %調(diào)用linsolve函數(shù)求解Ab %用另一種方法求解實(shí)驗(yàn)結(jié)果：6 對(duì)方程組求解a11*x1+a12*x2+a13*x3=b1a21*x1+a22*x2+a23*x3=b2a31*x1+a32*x2+a33*x3=b3syms a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3;A=a11,a12,a13;a21,a22,a23;a31

18、,a32,a33;b=b1;b2;b3;X=linsolve(A,b) %調(diào)用linsolve函數(shù)求的解XX=Ab %用左除運(yùn)算求解實(shí)驗(yàn)結(jié)果：7syms a b t x y z;f=sqrt(1+exp(x);diff(f) %未指定求導(dǎo)變量和階數(shù)，按缺省規(guī)則處理f=x*cos(x);diff(f,x,2) %求f對(duì)x的二階導(dǎo)數(shù)diff(f,x,3) %求f對(duì)x的三階導(dǎo)數(shù)f1=a*cos(t);f2=b*sin(t);diff(f2)/diff(f1) %按參數(shù)方程求導(dǎo)公式求y對(duì)x的導(dǎo)數(shù)實(shí)驗(yàn)結(jié)果：4.1.3 simulink的使用1在命令窗口中輸入：simulink（回車）得到如下simuli

19、nk模塊：2雙擊打開各模塊，選擇合適子模塊構(gòu)造控制系統(tǒng)，例如：3雙擊各子模塊可修改其參數(shù)，選擇Simulation菜單下的start命令運(yùn)行仿真，在示波器（Scope）中觀察結(jié)果。4.2 實(shí)驗(yàn)擴(kuò)展程序舉例：Simulink部分：逐一熟悉各模塊的使用，對(duì)下面的隨動(dòng)系統(tǒng)模型進(jìn)行仿真，實(shí)驗(yàn)報(bào)告中包含：連好的系統(tǒng)模型及用Scope觀測(cè)的結(jié)果其中：R(s)為階躍輸入，C(s)為輸出示波器顯示結(jié)果：4.4 實(shí)驗(yàn)心得Simulink是MATLAB中一個(gè)用于實(shí)現(xiàn)動(dòng)態(tài)系統(tǒng)建模、仿真與分析的仿真軟件集成環(huán)境軟件工具包。為便于仿真，simulink模塊庫提供了豐富的功能模塊。這些模塊分屬于不同功能、不同類別的子模塊

20、庫。熟悉這些模塊所在位置以及它們的功能和使用方法需要強(qiáng)加練習(xí)。實(shí)驗(yàn)五 MATLAB在控制系統(tǒng)分析中的應(yīng)用5.1 實(shí)驗(yàn)題目實(shí)驗(yàn)任務(wù)1掌握MATLAB在控制系統(tǒng)時(shí)間響應(yīng)分析中的應(yīng)用；2掌握MATLAB在系統(tǒng)根軌跡分析中的應(yīng)用； 3. 掌握MATLAB控制系統(tǒng)頻率分析中的應(yīng)用； 4. 掌握MATLAB在控制系統(tǒng)穩(wěn)定性分析中的應(yīng)用5.1.2 基本命令1. step 2. impulse 3. initial 4. lsim 5. rlocfind 6. bode 7. margin 8. nyquist 9. Nichols 10. cloop5.1.3 程序練習(xí)1. 求下面系統(tǒng)的單位階躍響應(yīng)%程序如

21、下：num=4 ; den=1 , 1 , 4 ;step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %計(jì)算峰值時(shí)間max(y) %計(jì)算峰值實(shí)驗(yàn)結(jié)果：2. 求如下系統(tǒng)的單位階躍響應(yīng)%程序如下：a=0,1;-6,-5;b=0;1;c=1,0;d=0;y,x=step(a,b,c,d);plot(y) 實(shí)驗(yàn)結(jié)果：3. 求下面系統(tǒng)的單位脈沖響應(yīng)：%程序如下：num=4 ; den=1 , 1 ,4 ;impulse(num,den)實(shí)驗(yàn)結(jié)果：4. 已知二階系統(tǒng)的狀態(tài)方程為：求系統(tǒng)的零輸入響應(yīng)和脈沖響應(yīng)。 %程序如下：a

22、=0 , 1 ; -10 , -2 ; b=0 ; 1 ; c=1 , 0 ; d=0 ; x0=1 ,0; subplot(1 , 2 , 1) ; initial(a , b , c ,d,x0) subplot(1 , 2 , 2) ; impulse(a , b , c , d)實(shí)驗(yàn)結(jié)果：5：系統(tǒng)傳遞函數(shù)為：輸入正弦信號(hào)時(shí)，觀察輸出信號(hào)的相位差。 %程序如下： num=1 ; den=1 , 1 ; t=0 : 0.01 : 10 ; u=sin(2*t) ; hold on plot(t , u , r) lsim(num , den , u , t)實(shí)驗(yàn)結(jié)果:6. 有一二階系統(tǒng)，求

23、出周期為4秒的方波的輸出響應(yīng) %程序如下： num=2 5 1;den=1 2 3;t=(0:.1:10);period=4;u=(rem(t,period)>=period./2);%看rem函數(shù)功能lsim(num,den,u,t);實(shí)驗(yàn)結(jié)果：7. 已知開環(huán)系統(tǒng)傳遞函數(shù)，繪制系統(tǒng)的根軌跡，并分析其穩(wěn)定性 %程序如下： num=1 2;den1=1 4 3;den=conv(den1,den1);figure(1)rlocus(num,den)k,p= rlocfind(num,den) figure(2)k=55;num1=k*1 2;den=1 4 3;den1=conv(den,

24、den);num,den=cloop(num1,den1,-1);impulse(num,den)title('impulse response (k=55) ')figure(3)k=56;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title('impulse response(k=56)')實(shí)驗(yàn)結(jié)果：8. 作如下系統(tǒng)的bode圖 %程序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d)實(shí)驗(yàn)結(jié)果：9

25、. 系統(tǒng)傳函如下 求有理傳函的頻率響應(yīng)，然后在同一張圖上繪出以四階伯德近似表示的系統(tǒng)頻率響應(yīng) %程序如下： num=1;den=conv(1 2,conv(1 2,1 2); w=logspace(-1,2); t=0.5;m1,p1=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(m1),w,20*log10(m2),'g-')

26、;grid on ; title('bode plot');xlabel('frequency');ylabel('gain');subplot(2,1,2);semilogx(w,p1,w,p2,'g-');grid on;xlabel('frequency');ylabel('phase');實(shí)驗(yàn)結(jié)果：10. 已知系統(tǒng)模型為 求它的幅值裕度和相角裕度 %程序如下： n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)實(shí)驗(yàn)結(jié)果：11. 二階系統(tǒng)為：令wn=1,分別

27、作出=2 , 1 , 0.707 , 0.5時(shí)的nyquist曲線。 %程序如下：n=1 ; d1=1 , 4 , 1 ; d2=1 , 2 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1; nyquist(n,d1) ; hold on nyquist(n,d2) ; nyquist(n,d3) ; nyquist(n,d4) ;實(shí)驗(yàn)結(jié)果：12. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為 繪制系統(tǒng)的Nyqusit圖，并討論系統(tǒng)的穩(wěn)定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis('square')實(shí)驗(yàn)結(jié)果：13. 分別由

28、w的自動(dòng)變量和人工變量作下列系統(tǒng)的nyquist曲線： %程序如下： n=1 ; d=1 , 1 ,0 ;nyquist(n ,d) ; %自動(dòng)變量n=1 ; d=1 , 1 ,0 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工變量實(shí)驗(yàn)結(jié)果：14. 一多環(huán)系統(tǒng)，其結(jié)構(gòu)圖如下，使用Nyquist頻率曲線判斷系統(tǒng)的穩(wěn)定性。 %程序如下： k1=16.7/0.0125;z1=0;p1=-1.25 -4 -16;num1,den1=zp2tf(z1,p1,k1);num,den=cloop(num1,den1);z,p,k=tf2zp(num,den);p fi

29、gure(1)nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2);實(shí)驗(yàn)結(jié)果：15. 已知系統(tǒng)為：作該系統(tǒng)的nichols曲線。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ;實(shí)驗(yàn)結(jié)果：16. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為：當(dāng)k=2時(shí)，分別作nichols曲線和波特圖。 %程序如下：num=1;den=conv(conv(1 0,1 1),0.5 1);subplot(1,2,1);nichols(num,den);grid; % nichols曲線

30、subplot(1,2,2);g=tf(num,den);bode(feedback(g,1,-1);grid; %波特圖實(shí)驗(yàn)結(jié)果：17. 系統(tǒng)的開環(huán)傳遞函數(shù)為：分別確定k=2和k=10時(shí)閉環(huán)系統(tǒng)的穩(wěn)定性。%程序如下：d1=1 , 3 , 2 , 0 ; n1=2 ;nc1 , dc1=cloop(n1 , d1 ,-1) ;roots(dc1)d2=d1 ; n2=10 ;nc2 , dc2=cloop(n2 , d2,-1) ; roots(dc2)實(shí)驗(yàn)結(jié)果：18. 系統(tǒng)的狀態(tài)方程為： 試確定系統(tǒng)的穩(wěn)定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;

31、0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根rank(ctrb(a,b)實(shí)驗(yàn)結(jié)果：5.2 實(shí)驗(yàn)心得 通過本次試驗(yàn)，熟悉了解了MATLAB中SIMULINK的應(yīng)用，特別是在自動(dòng)控制系統(tǒng)建模中，它非常的簡便易操作，相比輸入語句的建模來說，它形象而快速，不僅便于理解和分析，而且大大提高了仿真速度。同時(shí)，也鞏固了MATLAB的畫圖功能，同樣為之后對(duì)控制系統(tǒng)的仿真打下基礎(chǔ)，對(duì)于相關(guān)傳遞函數(shù)圖形的繪制有很大的幫助。實(shí)驗(yàn)六 連續(xù)系統(tǒng)數(shù)字仿真的基本算法6.1 實(shí)驗(yàn)題目6.1.1 實(shí)驗(yàn)任務(wù)1理解歐拉法和龍格-庫塔法的基本思想；2理解數(shù)值積分算法的計(jì)算精度、速度、穩(wěn)定性與步長的關(guān)系。6

32、.1.2 程序練習(xí)1. 取h=0.2，試分別用歐拉法、RK2法和RK4法求解微分方程的數(shù)值解，并比較計(jì)算精度。 注：解析解： %程序如下cleart(1)=0; % 置自變量初值y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和數(shù)值解的初值h=0.2; % 步長% 求解析解for k=1:5 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1);end% 利用歐拉法求解for k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);end% 利用RK

33、2法求解for k=1:5 k1=y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;end% 利用RK4法求解for k=1:5 k1=y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(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

34、*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;end% 輸出結(jié)果disp(' 時(shí)間 解析解 歐拉法 RK2法 RK4法')yt=t', y', y_euler', y_rk2', y_rk4'disp(yt)實(shí)驗(yàn)結(jié)果：2. 考慮如下二階系統(tǒng)： 在上的數(shù)字仿真解（已知：，），并將不同步長下的仿真結(jié)果與解析解進(jìn)行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法進(jìn)行計(jì)算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件采用RK4法進(jìn)行計(jì)算。程

35、序如下：clearh=input('請(qǐng)輸入步長h='); % 輸入步長M=round(10/h); % 置總計(jì)算步數(shù)t(1)=0; % 置自變量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分別對(duì)應(yīng)于為R=0和R=0.5)x1(1)=100; x2(1)=0; % 置狀態(tài)向量初值y_rk4_0(1)=x1(1); y_rk4_05(1)=x1(1); % 置數(shù)值解的初值 % 求解析解for k=1:M t(k+1)=t(k)+h; y_0(k+1)=100*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/

36、2).*cos(sqrt(3)/2*t(k+1)+100*sqrt(3)/3*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1);end% 利用RK4法求解% R=0for k=1:M k11=x2(k); k12=-x1(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_0(k+1)=x1(k+1);end% R=0.5for k=1:M k11=x2(k); k12=-x1(k)-x2(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2); k31=x2(k)+h*