系統(tǒng)仿真實(shí)驗(yàn)報告

1、系統(tǒng)仿真技術(shù)實(shí)驗(yàn)報告 信息科學(xué)與工程學(xué)院 目錄實(shí)驗(yàn)一 MATLAB中矩陣與多項(xiàng)式的基本運(yùn)算1實(shí)驗(yàn)二 MATLAB繪圖命令2實(shí)驗(yàn)三 MATLAB程序設(shè)計3實(shí)驗(yàn)四 MATLAB的符號計算與SIMULINK的使用5實(shí)驗(yàn)五 MATLAB在控制系統(tǒng)分析中的應(yīng)用8實(shí)驗(yàn)六 連續(xù)系統(tǒng)數(shù)字仿真的基本算法15實(shí)驗(yàn)一對不同命令操作訓(xùn)練1.eye(3)ans = 1 0 0 0 1 0 0 0 12.ones(2)ans = 1 1 1 1ones(3,2)ans = 1 1 1 1 1 1 3.zeros(3,2)ans = 0 0 0 0 0 04.rand(3,2)ans = 0.8147 0.9134 0.9

2、058 0.6324 0.1270 0.09755.diag(3) ans = 36. A=ones(3,3);B=rand(3,3); ans=Bans = 0.0357 0.6787 0.3922 0.8491 0.7577 0.6555 0.9340 0.7431 0.17127. A/Bans = 0.1195 1.5417 -0.3355 0.1195 1.5417 -0.33550.1195 1.5417 -0.3355 AB警告: 矩陣為奇異工作精度。 ans = NaN NaN NaN NaN NaN NaN Inf Inf -Inf inv(A)*B、B*inv(A)同上8.

3、roots(1)ans = Empty matrix: 0-by-1A*Bans = 1.8188 2.1796 1.2189 1.8188 2.1796 1.2189 1.8188 2.1796 1.2189>> B*Aans = 1.1067 1.1067 1.1067 2.2623 2.2623 2.2623 1.8483 1.8483 1.84839.who您的變量為:A B ans th x whos Name Size Bytes Class Attributes A 3x3 72 double B 3x3 72 double ans 3x3 72 double th

4、1x21 168 double x 1x21 336 double complex 10.size(A)ans = 3 3 length(A)ans = 3實(shí)驗(yàn)二1、二維直角坐標(biāo) t=-pi:0.3:pi;y=1./(1+exp(-t);subplot(221),plot(t,y,'b');xlabel('t');title('plot(t,y)');subplot(222),stem(t,y,'g');xlabel('t');title('plot(t,y)');subplot(223),semi

5、logy(t,y,'m');xlabel('t');title('plot(t,y)');subplot(224),stairs(t,y,'k');xlabel('t');title('plot(t,y)')實(shí)驗(yàn)三 function f=SY3()while 1 b=input('輸入一個自然數(shù):'); if b<2 fprintf('%d既不是質(zhì)數(shù)也不是合數(shù)',b); elseif isprime(b)=1 fprintf('%d是一個素數(shù)',

6、b); else fprintf('%d是一個合數(shù)，可以因式分解為:',b) fprintf('n') for n=2:sqrt(b) if rem(b,n)=0 num1=n; num2=b/n; fprintf('%d = %d×%d ',b,num1,num2) fprintf('n') end end end fprintf('n') y=input('繼續(xù)判斷？繼續(xù)按1，退出按0:'); if y=0 break; endend 命令窗口輸入 SY3輸入一個自然數(shù):7878是一個

7、合數(shù)，可以因式分解為:78 = 2×39 78 = 3×26 78 = 6×13 繼續(xù)判斷？繼續(xù)按1，退出按0:實(shí)驗(yàn)四1求矩陣對應(yīng)的行列式和特征根a=sym('a11 a12;a21 a22');da=det(a)ea=eig(a)da = a11*a22 - a12*a21 ea = a11/2 + a22/2 - (a112 - 2*a11*a22 + a222 + 4*a12*a21)(1/2)/2 a11/2 + a22/2 + (a112 - 2*a11*a22 + a222 + 4*a12*a21)(1/2)/2 2 求方程的解（包括精

8、確解和一定精度的解）r1=solve('x2-x-1')rv=vpa(r1)rv4=vpa(r1,4)rv20=vpa(r1,20)r1 = 5(1/2)/2 + 1/2 1/2 - 5(1/2)/2 rv = 1.6180339887498948482045868343656 -0.61803398874989484820458683436564 rv4 = 1.618 -0.618 rv20 = 1.6180339887498948482 -0.61803398874989484823 a=sym('a');b=sym('b');c=sym(

9、'c');d=sym('d'); %定義4個符號變量w=10;x=5;y=-8;z=11; %定義4個數(shù)值變量A=a,b;c,d %建立符號矩陣AB=w,x;y,z %建立數(shù)值矩陣Bdet(A) %計算符號矩陣A的行列式det(B) %計算數(shù)值矩陣B的行列式A = a, b c, d B = 10 5 -8 11 ans = a*d - b*c ans = 150 4 syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %對s展開collect(s,x) %對s按變量x合并同類項(xiàng)(無同類項(xiàng))factor(ans) % 對an

10、s分解因式5 對方程 AX=b求解A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=linsolve(A,b) %調(diào)用linsolve函數(shù)求解X = 138/2045 66/409 212/2045 Ab %用另一種方法求解ans = 138/2045 66/409 212/2045 6syms 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對x的二階導(dǎo)數(shù)diff(f,x,3) %求f對x的三階導(dǎo)數(shù)f1=a*cos(t);f2=b*sin(t);diff(f2)/d

11、iff(f1) %按參數(shù)方程求導(dǎo)公式求y對x的導(dǎo)數(shù) ans = exp(x)/(2*(exp(x) + 1)(1/2) ans = - 2*sin(x) - x*cos(x) ans = x*sin(x) - 3*cos(x) ans = -(b*cos(t)/(a*sin(t) 結(jié)果： 實(shí)驗(yàn)五1、求下面系統(tǒng)的單位階躍響應(yīng)%程序如下：num=5 ; den=1 , 1 , 4 ;(分子4改為5)step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %計算峰值時間max(y) %計算峰值 tp = 1.6579>

12、;> max(y) %計算峰值ans =1.80402. 求如下系統(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) 3. 求下面系統(tǒng)的單位脈沖響應(yīng)：%程序如下：num=4 ; den=1 , 1 ,4 ;impulse(num,den) 4. 已知二階系統(tǒng)的狀態(tài)方程為：求系統(tǒng)的零輸入響應(yīng)和脈沖響應(yīng)。 %程序如下：a=0 , 1 ; -10 , -2 ; b=0 ; 1 ; c=1 , 0 ; d=0 ; x0=1 ,0; subplot(1 , 2 , 1) ; initial(a , b , c ,d

13、,x0) subplot(1 , 2 , 2) ; impulse(a , b , c , d) 6. 有一二階系統(tǒng)，求出周期為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);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

14、,den) figure(2)k=55;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=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)') 8. 作如下系統(tǒng)的bode圖 %程

15、序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d) 9. 系統(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);semil

16、ogx(w,20*log10(m1),w,20*log10(m2),'g-');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'); 11. 二階系統(tǒng)為：令wn=1,分別作出=2 , 1 , 0.707 , 0.5時的nyquist曲線。 %程

17、序如下：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) ; 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') 13. 分別由w的自動變量和人工變量作下列系統(tǒng)的nyquist曲線： %程序如下： n=1 ; d=1

18、, 1 ,0 ;nyquist(n ,d) ; %自動變量n=1 ; d=1 , 1 ,0 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工變量 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 figure(1)nyquist(num,den)figure(2)num2,den2=cloop(n

19、um,den);impulse(num2,den2); 15. 已知系統(tǒng)為：作該系統(tǒng)的nichols曲線。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ; 16. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為：當(dāng)k=2時，分別作nichols曲線和波特圖。 %程序如下：num=1;den=conv(conv(1 0,1 1),0.5 1);subplot(1,2,1);nichols(num,den);grid; % nichols曲線subplot(1,2,2);g=tf(num,den);bode(feedback(g,1,-1);grid;

20、%波特圖17. 系統(tǒng)的開環(huán)傳遞函數(shù)為： 分別確定k=2和k=10時閉環(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)ans = -2.5214 + 0.0000i -0.2393 + 0.8579i -0.2393 - 0.8579ians = -3.3089 + 0.0000i 0.1545 + 1.7316i 0.1545 - 1.7316i18. 系統(tǒng)的狀態(tài)方程為： 試確

21、定系統(tǒng)的穩(wěn)定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根rank(ctrb(a,b) ans = 0 -1 -3ans = 3實(shí)驗(yàn)六1. 取h=0.2，試分別用歐拉法、RK2法和RK4法求解微分方程的數(shù)值解，并比較計算精度。 注：解析解： %程序如下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;

22、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% 利用RK2法求解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*k

23、1/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+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;end% 輸出結(jié)果disp(' 時間 解析解 歐拉法 RK2法 RK4法')yt=t', y', y_euler', y_rk2', y_rk4'disp(yt)程序運(yùn)行結(jié)果如下： 時間

24、 解析解 歐拉法 RK2法 RK4法 0 1.0000 1.0000 1.0000 1.0000 0.2000 1.1832 1.2000 1.1867 1.1832 0.4000 1.3416 1.3733 1.3483 1.3417 0.6000 1.4832 1.5315 1.4937 1.4833 0.8000 1.6125 1.6811 1.6279 1.6125 1.0000 1.7321 1.8269 1.7542 1.73212. 考慮如下二階系統(tǒng)： 在上的數(shù)字仿真解（已知：，），并將不同步長下的仿真結(jié)果與解析解進(jìn)行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法

25、進(jìn)行計算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件 采用RK4法進(jìn)行計算。程序如下：clearh=input('請輸入步長h='); % 輸入步長M=round(10/h); % 置總計算步數(shù)t(1)=0; % 置自變量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分別對應(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)=1

26、00*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/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*