中南大學(xué)系統(tǒng)仿真實(shí)驗(yàn)報告

1、本文格式為word版，下載可任意編輯中南大學(xué)系統(tǒng)仿真實(shí)驗(yàn)報告 ans = 試驗(yàn)一 matlab 中矩陣與多項式的基本運(yùn)算 試驗(yàn)任務(wù) 1. 了解 matlab 命令窗口和程序文件的調(diào)用。 2 熟識如下 matlab 的基本運(yùn)算： 矩陣的產(chǎn)生、數(shù)據(jù)的輸入、相關(guān)元素的顯示； 矩陣的加法、乘法、左除、右除； 特別矩陣：單位矩陣、 1 矩陣、 0 矩陣、對角陣、隨機(jī)矩陣的產(chǎn)生和 運(yùn)算； 多項式的運(yùn)算：多項式求根、多項式之間的乘除。 基本命令訓(xùn)練 1、 eye(2) ans = 1 0 0 1 eye(4) ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 2、 ones(2)

2、1 1 ans = 1 1 ones(4) ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ones(2,2) ans = 1 1 1 1 ones(2,3) ans = 1 1 1 1 1 1 ones(4,3) ans = 1 1 1 1 1 1 1 1 1 1 1 1 3、 zeros(2) 0 0 0 0 zeros(4) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 zeros(2,2) ans = 0 0 0 0 zeros(2,3) ans = 0 0 0 0 0 0 zeros(3,2) ans = 0 0 0 0 00

3、 4、隨機(jī)陣 rand(2,3) ans = 0.2785 0.9575 0.1576 0.5469 0.9649 0.9706 rand(2,3) ans = 0.9572 0.8003 0.4218 0.4854 0.1419 0.9157 5、 diag(5) ans = 5 diag(5,5) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diag(2,3) ans = 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 6、 (inv (a)為求 a 的逆矩陣)

4、b=5 3 1;2 3 8;1 1 1,inv(b) 5 3 1 2 3 8 1 1 1 ans = 0.6250 0.2500 -2.6250 -0.7500 -0.5000 4.7500 0.1250 0.2500 -1.1250 a=2 3;4 4,b=5 3;3 8,inv(a),inv(b);ab,a/b,inv(a)*b,b*inv(a) a = 2 3 4 4 b = 5 3 3 8 ans = -1.0000 0.7500 1.0000 -0.5000 ans = -2.7500 3.0000 3.5000 -1.0000 ans = 0.2258 0.2903 0.6452

5、0.2581 ans = -2.7500 3.0000 3.5000 -1.0000 ans = -2.0000 2.2500 5.0000 -1.7500 7、 p =1,-6,-72,-27, roots(p) p = 1 -6 -72 -27 ans = 12.1229 -5.7345 -0.3884 p=2,3,6,roots(p) p = 2 3 6 ans = -0.7500 + 1.5612i -0.7500 - 1.5612i 8、( a 為 n*n 的方陣) a=0 1 0;-4 4 0;-2 1 2,poly(a),b=sym(a),poly(b) a = 0 1 0 -4

6、 4 0 -2 1 2 ans = 1 -6 12 -8 b = 0, 1, 0 -4, 4, 0 -2, 1, 2 ans = x a 3-6*x a 2+12*x-8 9, 、( conv 是多項式相乘， deconv 是多項式相除) u=1 2 4 6 ,v=5 0 0 -6 7,conv(u,v) u = 1 2 4 6 v = 5 0 0 -6 7 ans = 5 10 20 24 -5 -10 -8 42 v=1 2 4 6 ,u=5 0 0 -6 7,deconv(u,v) v = 1 2 4 6 u = 5 0 0 -6 7 ans = 5 -10 10、（ 點(diǎn)乘是數(shù)組的運(yùn)算，

7、沒有點(diǎn)的乘是矩陣運(yùn)算 ） a = 2 5;3 4, b =3 1;4 7,a.*b,a*b a = 2 5 3 4 b = 3 1 4 7 ans = 6 5 12 28 ans = 26 37 25 31 a = 2 3; b = 4 7; a.*b = 8 21; a*b %錯誤 a*b = 29; 11、（ who 可以看到你用過的一些變量， 來了） who your variables are: a b a ans b p u whos name size bytes whos 是把該變量及所存儲的大小等信息都顯示出 class attributes 2x2 32 double b 2

8、x2 32 double a 1x2 16 double ans 1x2 16 double b 1x2 16 double p 1x3 24 double u 1x5 40 double v 1x4 32 double 12、 a=2 5 3;6 5 4,disp(a),size(a),length(a) a = 2 5 3 6 5 4 2 5 3 6 5 4 ans = 2 3 ans = 3 試驗(yàn)二 matlab 繪圖命令 試驗(yàn)任務(wù) 熟識 matlab 基本繪圖命令，把握如下繪圖方法： 1 坐標(biāo)系的選擇、圖形的繪制； 2 圖形注解（題目、標(biāo)號、說明、分格線）的加入; 3 圖形線型、符號、

9、顏色的選取 基本命令訓(xùn)練 1、t=0:pi/360:2*pi; x=cos(t)+ cos(t*4); y=si n( t)+ sin (t*4); xlabel(x 軸);ylabel(y 軸); plot(y,x),grid; 2、 t=0:0.1:100; x=3*t;y=4*t;z=si n( 2*t); plot3(x,y,z, g:) 15 i 0 5 0 05 1 1 5 2 3、x = linspace(-2*pi,2*pi,40); y=si n( x); stairs(x,y)4、 t=0:pi/360:2*pi; x=cos(t)+ cos(t*4) + sin (t*4

10、); y=si n( t)+ si n( t*4); plot(y,x, r:); xlabel(x 軸);ylabel(y 軸); 6、th=0:pi/20:2*pi; x=exp(j*th); plot(real(x),imag(x),r-.); grid; text(0,0,中心); 5、 th=0:pi/1000:2*pi; r=cos(2*th); polar(th,r); title( 四葉草圖 ) 270 四葉草圖 10 7、x=-2:0.01:2; 8、y=-2:0.01:2; 9、x,y = meshgrid(x,y); z = y.*exp(-x. a 2-y. a 2);

11、 c,h = con tour(x, yz); set(h,showtext,o n,textstep,get(h,levelstep)*2)_ 1os i ,5 2 n.s o o.s 8、x = 0:0.2:10; y = 2*x+3; subplot(411);plot(x,y); grid;title(y 的原函數(shù)); subplot(412) ;semilogy(x,y); grid;title(對 y 取對數(shù)); 丫的原畫數(shù) 40 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - l| i p | il | i| i 九 _ 1- _ _ i _ l _ : ： _ j

12、_ _ l _ u i| i |l i , il _ - j i i _ h- _ i i q 【 i 1 f i i i i ii i i it 10 1 1 2 3 4 5 6 r 6 9 10 , 對 y 取對數(shù) 對弋觀對數(shù) subplot(413) ;semilogx(x,y); 10 10 id 10 1lz 10 w 10 10 40 20 0 對好對數(shù) 2 把握循環(huán)、分支語句的編寫，學(xué)會使用 look for 、 help 命令 grid;title(對 x 取對數(shù)); subplot(414) ;loglog(x,y);grid; title(對 xy 均取對數(shù)); 9、x =

13、 -3:0.3:3; bar(x,exp(-x.*x),g) 試驗(yàn)三 matlab 程序設(shè)計 試驗(yàn)任務(wù) 1 熟識 matlab 程序設(shè)計的方法和思路； 程序舉例 1、 f=1,1; i=1; while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1; end f,i f = columns 1 through 14 1 1 2 3 5 8 13 21 34 55 89 144 233 377 columns 15 through 16 610 987 i = 15 2、 m=3; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j

14、-1); end end format rat 1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5 (分?jǐn)?shù)格式形式。用有理數(shù)靠近顯示數(shù)據(jù) ) m=5; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j-1); end end format rat a a = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7 1/5 1/6 1/7 1/8 3、程序中沒有 format rat 命令時，假如上次運(yùn)行結(jié)果沒有清除，輸出的結(jié)果就是上次運(yùn)行的 結(jié)果！但是運(yùn)用 clear 命令清晰之前的

15、運(yùn)行結(jié)果之后就會正常運(yùn)行。 4、 x=input（請輸入 x 的值:）; if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); 1/4 1/5 1/6 end y 請輸入 x 的值 :2 y = 2391/647 x=input( 請輸入 x 的值:); if x=10 y= fprintf ( 不在定義域內(nèi)，請重新輸入： );return else y=1/(x-10); end y 請輸入 x 的值 :2 -1/8 5、 p=0 0 0 1 3 0 2 0 0 9; for i=1:length(p), if p(1)=0,p=p

16、(2:length(p); end ; end ; p p = columns 1 through 5 1 3 0 2 0 columns 6 through 7 0 9 p=0 0 0 1 3 0 2 0 0 9; p(p=0)=;p p = 1 3 2 9 6、 e2(500) ans = 1 1 2 3 5 8 13 21 34 55 89 144 233 377 lookfor ffibno e2 - ffibno 計算斐波那契亞數(shù)列的函數(shù)文件 help e2 ffibno 計算斐波那契亞數(shù)列的函數(shù)文件 n 可取任意自然數(shù) 程序如下 （用法： lookfor 關(guān)鍵詞 在全部 m 文件中

17、找 關(guān)鍵詞 ，比如： lookfor max （即查找關(guān)鍵詞 max） 其實(shí)就和我們平常用 ctrl+f 來查找關(guān)鍵詞是一樣的 而 help 是顯示 matlab 內(nèi)置的關(guān)心信息 用法： help 命令，比如 help inv ，作用就是調(diào)用 inv 這個命令的關(guān)心） 程序設(shè)計題 用一個 matlab 語言編寫一個程序：輸入一個自然數(shù)，推斷它是否是素數(shù)， 假如是，輸出 it is one prime ，假如不是，輸出 it is n ot o ne prime. 。要求通 過調(diào)用子函數(shù)實(shí)現(xiàn)。 最好能具有如下功能： 設(shè)計較好的人機(jī)對話界面， 程序中 含有提示 性的輸入輸出語句。 能實(shí)現(xiàn)循環(huán)操作，

18、 由操輸入相關(guān)命令來掌握 是否連續(xù)進(jìn)行素數(shù)的推斷。 假如操盼望停止這種推斷， 則可以退出程序。 假如所輸入的自然數(shù)是一個合數(shù)， 除了給出其不是素數(shù)的結(jié)論外， 還應(yīng)給出至少 一種其因數(shù)分解形式。例：輸入 6 ，由于 6 不是素數(shù)。則程序中除了有 it is not one prime 的結(jié)論外，還應(yīng)有： 6=2*3 的說明。 function sushu while 1 x=input（ 請輸入一個自然數(shù) ）; if x2 disp（ 既不是質(zhì)數(shù)又不是合數(shù) ）; else if isprime(x)=1 disp（ 這是一個素數(shù) ）; else disp( 這是一個合數(shù)，可以因式分解為： ) f

19、or n=2:sqrt(x) if rem(x,n)=0 num3=x; num1=n; num2=x/n; disp(num2str(num3),=,num2str(num1),x,num2str(num2) end end end end y=input( 是否連續(xù)推斷？連續(xù)請按 1 ,按任意鍵退出 ： ) ； if y=1 break end end 試驗(yàn)四 matlab 的符號計算與 simulink 的使用 試驗(yàn)任務(wù) 1 .把握 matlab 符號計算的特點(diǎn)和常用基本命令; 2 .把握 simulink 的使用。 程序舉例 1 . 求矩陣對應(yīng)的行列式和特征根 a=sym(a11 a1

20、2;a21 a22); da=det(a) ea=eig(a) da = a11*a22-a12*a21 ea = 1/2*a11+1/2*a22+1/2*(a11 a 2-2*a11*a22+a22 a 2+4*a12*a21) a (1 /2) 1/2*a11+1/2*a22-1/2*(a11 a 2-2*a11*a22+a22 a 2+4*a12*a21) a (1 /2) a=sym(2 3;1 5); da=det(a) ea=eig(a) da = 7 ea = 7/2+1/2*21a(1/2) 7/2-1/2*21a(1/2) 2. 求方程的解(包括精確解和肯定精度的解) r1=

21、solve(xa2+x-1) rv=vpa(r1) rv4=vpa(r1,4) rv20=vpa(r1,20) r1 = 1/2*5 八( 1/2)-1/2 -1/2*5 八( 1/2)-1/2 rv = .6894848260 -1.689484826 rv4 = .6180 -1.618 rv20 = .689484820 -1.68948482 3 a=sym( a );b=sym( b );c=sym( w=10;x=5;y=-8;z=11; a=a,b;c,d b=w,x;y,z det(a) det(b) a = a, b c );d=sym( d ); %定義 4 個符號變量 %

22、定義 4 個數(shù)值變量 %建立符號矩陣 a %建立數(shù)值矩陣 b %計算符號矩陣 a 的行列式 %計算數(shù)值矩陣 b 的行列式 c, d x 10 -8 11 ans = a*d-b*c ans = 150 4. syms x y;s=(-7*x a 2-8*y a 2)*(-x a 2+3*y a 2); expand(s) % 對 s 綻開 collect(s,x) % 對 s 按變量 x 合并同類項 （ 無同類項 ） factor(ans) % 對 ans 分解因式 ans = 7*xa4-13*xa2*ya2-24*ya4 ans = 7*xa4-13*xa2*ya2-24*ya4 ans

23、= (8*ya2+7*xa2)*(xa2-3*ya2) 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ù)求解 ab % 用另一種方法求解 0.0675 0.1614 diff(f) %未指定求導(dǎo)變量和階數(shù)，按缺省規(guī)章處理 0.1037 ans = 0.0675 0.1614 0.1037 6 對方程組求解 a11*x1+a12*x2+a13*x3=b1 a21*x1+a22*x2+a23*x3=b2 a31*x1+a32*x2+a33*x3=b3 syms a11 a12 a13 a2

24、1 a22 a23 a31 a32 a33 b1 b2 b3; a=a11,a12,a13;a21,a22,a23;a31,a32,a33; b=b1;b2;b3; xx=ab %用左除運(yùn)算求解 (x=linsolve(a,b) %調(diào)用 linsolve 函數(shù)求的解 ) xx = (a12*a23*b3-a12*b2*a33+a13*a32*b2-a13*a22*b3+b1*a22*a33-b1*a32*a23)/( a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a 23) -(a11*a23*b3-a1

25、1*b2*a33-a21*a13*b3-a23*a31*b1+b2*a31*a13+a21*b1*a33)/ (a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12* a23) (a32*a21*b1-a11*a32*b2+a11*a22*b3-a22*a31*b1-a12*a21*b3+a31*a12*b2)/( a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a 23) 7 syms a b t x y z; f=sqrt

26、(1+exp(x); %求 f 對 x 的二階導(dǎo)數(shù) %求 f 對 x 的三階導(dǎo)數(shù) %按參數(shù)方程求導(dǎo)公式求 y 對 x 的導(dǎo)數(shù) ans = 1/2/(1+exp(x) a (1/2)*exp(x) ans = -2*si n(x)-x*cos(x) ans = -3*cos(x)+x*si n(x) ans = -b*cos(t)/a/si n(t) 三、 simulink 的使用 f=x*cos(x); diff(f,x,2) diff(f,x,3) f1=a*cos(t);f2=b*si n(t); diff(f2)/diff(f1) diff(f) %未指定求導(dǎo)變量和階數(shù)，按缺省規(guī)章處理

27、仿真圖: 波形圖: 其中: g i ( r 0.14 試驗(yàn)五 matlab 在掌握系統(tǒng)分析中的應(yīng)用 試驗(yàn)任務(wù) 1. 把握 matlab 在掌握系統(tǒng)時間響應(yīng)分析中的應(yīng)用; 2. 把握 matlab 在系統(tǒng)根軌跡分析中的應(yīng)用； 3. 把握 matlab 掌握系統(tǒng)頻率分析中的應(yīng)用； 4. 把握 matlab 在掌握系統(tǒng)穩(wěn)定性分析中的應(yīng)用 基本命令 1. step 2. impulse 3. in itial 4. isim 5. rlocfi nd 6. bode 7. margin 8. nyquist 9. nichols 10. cloop 程序舉例 1. 求下面系統(tǒng)的單位階躍響應(yīng) 1.5 g

28、(s) t num=4 ; den=1 , 1 , 4; step( num , den) y , x , t=step( num , den); tp=spli ne(y , t , max(y) % 計算峰值時間 max(y) % 計算峰值 tp = step response 0.5 time (sec) p m 1.6062 ans = 1.4441 0.18 0.16 2. 求如下系統(tǒng)的單位階躍響應(yīng) x i 0 1 x i 0 u x 2 6 5 x 2 1 a=0,1;-6,-5;b=0;1;c=1,0;d=0; y,x=step(a,b,c,d); plot(y) 1,0 x i

29、 x 2 3. 求下面系統(tǒng)的單位脈沖響應(yīng): g(s) 兒 num=4 ; den=1 , 1 ,4; impulse( nu m,de n) time (sec) response to initial conditions 4. 已知二階系統(tǒng)的狀態(tài)方程為: e u p m a 0.25 0.2 0.15 0.1 0.05 0 -0.05 l i 1 - - r - 1 廠、 / impulse response c=1 , 0 ； d=0; xo=1 ,0 ； subplot(1 , 2,1); in itial(a , b , c ,d,x0) subplot(1 , 2,2); impu

30、lse(a , b , c , d) 5 :系統(tǒng)傳遞函數(shù)為: g(s) 占 輸入正弦信號時，觀看輸出 信號的相位差。 num=1 ; den=1 ,1; t=0 : 0.01 : 10 ; u=s in( 2*t) ; hold on plot(t,u, r) lsim( nu m,de n,u,t) real axis 6. 有一二階系統(tǒng)，求出周期為 4 秒的方波的輸出響應(yīng) 2s 2 5s 1 2 s 2s 3 num=2 5 1; den=1 2 3; t=(0:.1:10); period=4; u=(rem(t,period)=period./2); % 看 rem 函數(shù)功能 lsim

31、( nu m,de n,u,t); 7. 已知開環(huán)系統(tǒng)傳遞函數(shù)，繪制系統(tǒng)的根軌跡，并分析其穩(wěn)定性 g(s) k(s 2) (s 2 4s 3) 2 num=1 2; den 仁 1 4 3; den=conv (de n1,de n1); figure(1) -6 rlocus(num,den) k,p= rlocfi nd(nu m,de n) l l l 匚一 l c - - - r t r r -10 -8 -6 -4 -2 4 6 root locus 0 2 2 2 0 0 2 2 - - cpxa vyanma n -8 2.5 g(s) lin ear simulati on r

32、esults1.5 e 0.5 -0.5 -1 -1.5 -2 time (sec) p m impulse response (k=55) 5 5 figure(2) k=55; num 仁 k*1 2; den=1 4 3; den 1=c onv (de n,den); nu m,de n=cloop( nu m1,de n1,-1); impulse( nu m,de n) title(impulse resp onse (k=55) -1.5 -1 0 200 400 600 time (sec) 800 1000 1200 5 5 0 0 5 5 0 0 o - - 6 x 10

33、8 ! - impulse resp on se(k=56) figure(3) k=56; num 仁 k*1 2; den=1 4 3; den 1=c onv (de n,den); nu m,de n=cloop( nu m1,de n1,-1); impulse( nu m,de n) -2 -4 -6 -8 l 0 500 1000 time (sec) title(impulse resp on se(k=56) 1500 2021 2500 select a point in the graphics win dow selected_po int = -2.5924 - 0.

34、0248i 0.7133 -3.4160 -2.5918 -0.9961 + 0.4306i -0.9961 - 0.4306i bode diagram 8. 作如下系統(tǒng)的 bode 圖 g(s) n=1 , 1 ; d=1 , 4,11 , 7; bode( n , d),grid on frequency (rad/sec) s 1 s 3 4s 2 11s 7 9. 系統(tǒng)傳函如下 g(s) s 1 0.5s e (s 2) 3 求有理傳函的頻率響應(yīng)，然后在同一張圖上繪出以四階伯德近似表示的系統(tǒng)頻 率響應(yīng) num=1;de n=co nv(1 2,co nv(1 2,1 2); -20

35、w=logspace(-1,2); t=0.5; ylabel(gai n); subplot(2,1,2); semilogx(w,p1,w,p2,g-);grid on; xlabel(freque ncy); ylabel(phase); 10. 已知系統(tǒng)模型為 求它的幅值裕度和相角裕度 n=3.5; d=1 2 3 2; gm,pm,wcg,wcp=margi n(n,d) g(s) 3.5 s 3 2s 2 3s 2 bode plot m1,p1=bode( nu m,de n, 2); p1=p1-t*w*180/pi; n 2,d2=pade(t,4); nu mt=c onv

36、(n2,nu m); den t=(c onv (de n, d2); 2 freque xlabel(freque gm = 1.1433 pm = 7.1688 wcg = 1.7323 wcp = 1.6541 nyq uist (n, d1); hold on nyq uist (n, d2) ; nyq uist (n, d3) ; nyq uist (n, d4); nyquist diagram 80 | - - - - 11. 二階系統(tǒng)為: g(s) n s 2 2 n s s 令 wn=1, 分別作出 e =2,1 , 0.707 , c a 5 y a 時的 nyquist

37、 曲線。 m n=1; d 仁 1 , 4,1 ; d2=1 , 2 , 1 ； d3=1 , 1.414,1; d4=1,1,1; nyquist diagram real axis 12. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為 -20 _一 - 一 一 -一_一_ - 60 00 5 5 s s 2) s s 3 3 2 2 s s v vr ra an nk ky ya a 14. 一多環(huán)系統(tǒng)，其結(jié)構(gòu)圖如下，使用 nyquist 頻率曲線推斷系統(tǒng)的穩(wěn)定性。 16.7s (0.85s 1)(0.25s 1)(0.0625s 1)繪制系統(tǒng)的 nyqusit 圖, 并爭論系統(tǒng)的穩(wěn)定性 . g=tf(100

38、0,co nv(1,3,2,1,5); nyquist(g);axis(square) 13. 分別由 w 的自動變量和人工變量作下列系統(tǒng)的 nyquisf 曲線： m g(s) 1 s(s 1) n=1 ; d=1 , 1 ,0; nyquist(n ,d) ; % 自動變量 n=1 ; d=1 , 1 ,0; w=0.5 : 0.1 : 3; nyq uist (n , d , w) ; % 人工變量 -2 -1 -1.5 -0.9 -0.8 -0.7 -0.6 -0.5 real axis -0.4 -0.3 -0.2 -0.1 5 5 - -0 0 s si ix xa a u u a

39、 an nt t9 9a a 卩 g(s) 40 30 20 10 0 figure(2) nu m2,de n2=cloop( nu m,de n); impulse( nu m2,de n2); k1=16.70.0125;z1=0; p1=-1.25 -4 -16; nu m1,de n1=zp2tf(z1,p1,k1); nu m,de n=cloop( nu m1,de n1); z,p,k=tf2zp( nu m,de n);p figure(1) nyq uist (nu m,de n) -2 -1 nyquist diagram -0.5 0 0.5 real axis 1 1

40、.5 1 0.5 s a 0 n g m -0.5 -1 -1.5 1.5 20 impulse resp onse -10.5969 +36.2148i -10.5969 -36.2148i -0.0562 15. 已知系統(tǒng)為: eanmlpm -5 -10 -15 0 0.1 0.2 0.3 time (sec) nichols chart 0.4 0.5 0.6 10 -2 10 0 10 2 open-loop phase (deg) frequency (rad/sec) g(s) s(s 1) 作該系統(tǒng)的 nichols 曲線。 n=1 ; d=1 , 1 , 0; ni chol

41、s( n , d); 16. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為: g(s) k s(s 1)(s 2) 當(dāng) k=2 時，分別作 nichols 曲線和波特圖 num=1; den=conv (co nv(1 0,1 1),0.5 1); subplot(1,2,1); ni chols( nu m,de n);grid; b % n ichols 曲線 g subplot(1,2,2); 2 g=tf(nu m,de n); bode(feedback(g,1,-1);grid; nichols chart koflr ea mhnaam msuvesab 50 50 - - 00 80 -270 b

42、ode diagram 90 - - % 波特圖 17. 系統(tǒng)的開環(huán)傳遞函數(shù)為: 分別確定 k=2 和 k=10 時閉環(huán)系統(tǒng)的穩(wěn)定性 d 仁 1 , 3,2,0 ; n 仁 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.2393 + 0.8579i -0.2393 - 0.8579i ans = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316i 18. 系統(tǒng)的狀態(tài)方程為：

43、x 1 4 3 0 x 1 1 x 2 1 0 0 x 2 0 u x 3 0 1 0 x 3 0 x 1 y 0 1 2 x 2 x 3 g(s) k s(s 1)(s 2) 試確定系統(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) % 求特征根 ran k(ctrb(a,b) ans = 0 -1 -3 ans = 3 試驗(yàn)六連續(xù)系統(tǒng)數(shù)字仿真的基本算法 試驗(yàn)任務(wù) 1 .理解歐拉法和龍格 - 庫塔法的基本思想； 2 理解數(shù)值積分算法的計算精度、速度、穩(wěn)定性與步長的關(guān)系; 程序舉例 1. 取 h=0.2 ,試

44、分別用歐拉法、 rk2 法和 rk4 法求解微分方程的數(shù)值解，并 比較計算精度。 注：解析解：y(t) .1 2t clear t(1)=0 ； y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; h=0.001; % 步長修改為 0.001 for k=1:5 y(t) y(t) y(0) 2t y(t) 1 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 for k=1:5

45、 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 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*k22)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+

46、h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6; end disp( 時間 解析解 yt=t, y, y_euler, y_rk2, y_rk4; disp(yt) 歐拉法 rk2 法 rk4 法 ) 時間 解析解 歐拉法 rk2 法 rk4法0 1.0000 1.0000 1.0000 1.0000 0.0010 1.0010 1.0010 1.0010 1.0010 0.0020 1.0020 1.0020 1.0020 1.0020 0.0030 1.0030 1.0030 1.0030 1.0030 0.0040 1.0040 1.00

47、40 1.0040 1.0040 0.0050 1.0050 1.0050 1.0050 1.0050 y(t) 2ry(t) y(t) 0 在 0 t 10 上的數(shù)字仿真解(已知： y(0) y(0) 0)， 并將不同步長下的仿真結(jié)果與解析解進(jìn)行精度比較。 說明： 已知該微分方程的解析解分別為： 100， yt y t 100cost ( 當(dāng) r 0) 100e 2t cs 仝 t 10 e*s in 仝 t 2 3 2 ( 當(dāng) r 0.5) 采納 rk4 法進(jìn)行計算，選擇狀態(tài)變量: 2. 考慮如下二階系統(tǒng)： x 1 y x 2 y 則有如下狀態(tài)空間模型及初值條件 x 1 x 2 x 1 (0) 100 x 2 x 1 2rx 2 x 2 (0) 0 y x 1 采納 rk4 法進(jìn)行計算。 clear h=input( 請輸入步長 h=); m=round(10/h); t(1)=0; y_0(1)=100; y_05(1)=100; 和 y_05 分別對應(yīng)于為 r=0 和 r=0.5) x1(1)=100; x2(1)=0; y_rk4_0(1)=x