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

1、實(shí)驗(yàn)一 MATLAB中矩陣與多項(xiàng)式的基本運(yùn)算實(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)式之間的乘除?；久钣?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 12、 ones(2)ans = 1 1 1 1 ones(4)ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2、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 13、 zeros(2)ans = 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 00 04、隨機(jī)陣 rand(2,3)ans = 0.2785 0.9575 0.15760.5469 0.9649 0.97

3、06 rand(2,3)ans =0.9572 0.8003 0.4218 0.4854 0.1419 0.91575、 diag(5)ans = 5 diag(5,5)ans = 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 06、（inv（A）為求A的逆矩陣） B=5 3 1;2 3 8;1 1 1,inv(B)B = 5 3 1 2 3 8 1 1 1ans = 0.6250 0.2500 -2.625

4、0 -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 4B = 5 3 3 8ans = -1.0000 0.7500 1.0000 -0.5000ans = -2.7500 3.0000 3.5000 -1.0000ans = 0.2258 0.2903 0.6452 0.2581ans = -2.7500 3.0000 3.5000 -1.0000ans = -2.0000 2.25005.0000 -1.75007

5、、 p =1,-6,-72,-27, roots(p)p = 1 -6 -72 -27ans = 12.1229 -5.7345 -0.3884 p=2,3,6,roots(p)p = 2 3 6ans = -0.7500 + 1.5612i -0.7500 - 1.5612i8、（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 4 0 -2 1 2ans = 1 -6 12 -8 B = 0, 1, 0 -4, 4, 0 -2, 1, 2 ans = x3-6*x2+12*x-89,、（conv是多項(xiàng)式

6、相乘，deconv是多項(xiàng)式相除） u=1 2 4 6 ,v=5 0 0 -6 7,conv(u,v)u = 1 2 4 6v = 5 0 0 -6 7ans = 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 6u = 5 0 0 -6 7ans = 5 -1010、（點(diǎn)乘是數(shù)組的運(yùn)算，沒有點(diǎn)的乘是矩陣運(yùn)算） a = 2 5;3 4, b =3 1;4 7,a.*b,a*ba = 2 5 3 4b = 3 1 4 7ans = 6 5 12 28ans = 26 3725 31 a = 2 3; b =

7、4 7;a.*b = 8 21;a*b %錯(cuò)誤a*b = 29;11、（who 可以看到你用過的一些變量，whos是把該變量及所存儲(chǔ)的大小等信息都顯示出來了） whoYour variables are:A B a ans b p u v whos Name Size Bytes Class Attributes A 2x2 32 double B 2x2 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;

8、6 5 4,disp(A),size(A),length(A)A = 2 5 3 6 5 4 2 5 3 6 5 4ans = 2 3ans = 3實(shí)驗(yàn)二 MATLAB繪圖命令實(shí)驗(yàn)任務(wù) 熟悉MATLAB基本繪圖命令，掌握如下繪圖方法： 1坐標(biāo)系的選擇、圖形的繪制； 2圖形注解（題目、標(biāo)號(hào)、說明、分格線）的加入； 3圖形線型、符號(hào)、顏色的選取。基本命令訓(xùn)練1、t=0:pi/360:2*pi;x=cos(t)+ cos(t*4);y=sin(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=

9、sin(2*t);plot3(x,y,z,g:)3、x = linspace(-2*pi,2*pi,40);y=sin(x);stairs(x,y) 4、t=0:pi/360:2*pi;x=cos(t)+ cos(t*4) + sin(t*4);y=sin(t)+ sin(t*4);plot(y,x,r:);xlabel(x軸);ylabel(y軸); 5、th=0:pi/1000:2*pi;r=cos(2*th);polar(th,r);title(四葉草圖)6、th=0:pi/20:2*pi;x=exp(j*th);plot(real(x),imag(x),r-.) ;grid; text

10、(0,0,中心) ; 7、 x=-2:0.01:2;8、 y=-2:0.01:2;9、 X,Y = meshgrid(x,y);Z = Y.*exp(-X.2-Y.2);C,h = contour(X,Y,Z);set(h,ShowText,on,TextStep,get(h,LevelStep)*2)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(對(duì)y取對(duì)數(shù));subplot(413) ;semilogx(x,y); grid;

11、title(對(duì)x取對(duì)數(shù));subplot(414) ;loglog(x,y);grid;title(對(duì)xy均取對(duì)數(shù)); 9、x = -3:0.3:3;bar(x,exp(-x.*x),g) 實(shí)驗(yàn)三 MATLAB程序設(shè)計(jì)實(shí)驗(yàn)任務(wù) 1熟悉MATLAB程序設(shè)計(jì)的方法和思路； 2掌握循環(huán)、分支語句的編寫，學(xué)會(huì)使用look for、help命令。程序舉例1、f=1,1;i=1;while f(i)+f(i+1)m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endendformat rataa = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1

12、/3 1/4 1/5 1/6 (分?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); endend format rataa = 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命令時(shí)，如果上次運(yùn)行結(jié)果沒有清除，輸出的結(jié)果就是上次運(yùn)行的結(jié)果！但是運(yùn)用clear命令清楚之前的運(yùn)行結(jié)果之后就會(huì)正常運(yùn)行。4、x=input(請(qǐng)輸入x的值:); if x=10 y=cos(x+1)+sqr

13、t(x*x+1); else y=x*sqrt(x+sqrt(x); end y請(qǐng)輸入x的值:2y = 2391/647x=input(請(qǐng)輸入x的值:); if x=10 y=fprintf(不在定義域內(nèi)，請(qǐng)重新輸入：);return else y=1/(x-10); end y請(qǐng)輸入x的值:2y = -1/85、p=0 0 0 1 3 0 2 0 0 9;for i=1:length(p),if p(1)=0,p=p(2:length(p); end;end;pp = Columns 1 through 5 1 3 0 2 0 Columns 6 through 7 0 9 p=0 0 0

14、1 3 0 2 0 0 9;p(p=0)=;pp = 1 3 2 96、 e2(500)ans = 1 1 2 3 5 8 13 21 34 55 89 144 233 377 lookfor ffibnoe2 - ffibno 計(jì)算斐波那契亞數(shù)列的函數(shù)文件 help e2 ffibno 計(jì)算斐波那契亞數(shù)列的函數(shù)文件 n可取任意自然數(shù) 程序如下（用法： lookfor 關(guān)鍵詞在所有M文件中找“關(guān)鍵詞”，比如：lookfor max（即尋找關(guān)鍵詞“max”）其實(shí)就和我們平時(shí)用CTRL+F來查找“關(guān)鍵詞”是一樣的而help是顯示matlab內(nèi)置的幫助信息 用法：help 命令，比如 help in

15、v ，作用就是調(diào)用inv這個(gè)命令的幫助）程序設(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ù)，除了給出其不是素?cái)?shù)的結(jié)論外，還應(yīng)給出至少一種其因數(shù)分解形式。例：輸入 6， 因?yàn)?不是素?cái)?shù)。則程序中除了有“It is not on

16、e prime”的結(jié)論外，還應(yīng)有：“6=2*3”的說明。function sushuwhile 1 x=input( 請(qǐng)輸入一個(gè)自然數(shù));if x a=sym(a11 a12;a21 a22);da=det(a)ea=eig(a)da = a11*a22-a12*a21 ea = 1/2*a11+1/2*a22+1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2) 1/2*a11+1/2*a22-1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2)a=sym(2 3;1 5);da=det(a)ea=eig(a)da =7ea = 7/2

17、+1/2*21(1/2) 7/2-1/2*21(1/2)2. 求方程的解（包括精確解和一定精度的解） r1=solve(x2+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/2rv =rv4 = .6180 -1.618rv20 = 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的行

18、列式det(B) %計(jì)算數(shù)值矩陣B的行列式A = a, b c, dB = 10 5 -8 11 ans = a*d-b*cans = 1504. syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %對(duì)s展開collect(s,x) %對(duì)s按變量x合并同類項(xiàng)(無同類項(xiàng))factor(ans) % 對(duì)ans分解因式ans =7*x4-13*x2*y2-24*y4ans = 7*x4-13*x2*y2-24*y4 ans = (8*y2+7*x2)*(x2-3*y2)5. 對(duì)方程 AX=b求解 A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=

19、linsolve(A,b) %調(diào)用linsolve函數(shù)求解Ab %用另一種方法求解X = 0.0675 0.1614 0.1037ans = 0.0675 0.16140.10376 對(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,a32,a33;b=b1;b2;b3;XX=Ab %用左除運(yùn)算求解(X=linsolve(A,b) %調(diào)用

20、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*a23) -(a11*a23*b3-a11*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*a

21、23) (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*a23)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(f

22、1) %按參數(shù)方程求導(dǎo)公式求y對(duì)x的導(dǎo)數(shù)ans = 1/2/(1+exp(x)(1/2)*exp(x) ans = -2*sin(x)-x*cos(x) ans = -3*cos(x)+x*sin(x) ans = -b*cos(t)/a/sin(t)三、SIMULINK的使用G1(s)G2(s)R(s)C(s)其中：R(s)為階躍輸入，C(s)為輸出 仿真圖：波形圖：實(shí)驗(yàn)五 MATLAB在控制系統(tǒng)分析中的應(yīng)用實(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)

23、穩(wěn)定性分析中的應(yīng)用基本命令 1. step 2. impulse 3. initial 4. lsim 5. rlocfind 6. bode 7. margin 8. nyquist 9. Nichols 10. cloop程序舉例1. 求下面系統(tǒng)的單位階躍響應(yīng) 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ì)算峰值tp = 1.6062ans =1.44412. 求如下系統(tǒng)的單位階躍響應(yīng) a=0,1;-6,-5;b=0;1;c=

24、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,x0)subplot(1 , 2 , 2) ; impulse(a , b , c , d)5：系統(tǒng)傳遞函數(shù)為：輸入正弦信號(hào)時(shí)，觀察輸出 信號(hào)的相位差。 num=1 ;

25、den=1 , 1 ;t=0 : 0.01 : 10 ;u=sin(2*t) ; hold onplot(t,u, r)lsim(num,den,u,t)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= rl

26、ocfind(num,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)Select a point in the graphic

27、s windowselected_point = -2.5924 - 0.0248ik = 0.7133p = -3.4160 -2.5918 -0.9961 + 0.4306i -0.9961 - 0.4306i8. 作如下系統(tǒng)的bode圖 n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d),grid on9. 系統(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

28、=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-);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);10. 已知系統(tǒng)模型為 求它

29、的幅值裕度和相角裕度 n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)Gm = 1.1433Pm = 7.1688Wcg = 1.7323Wcp =1.654111. 二階系統(tǒng)為：令wn=1,分別作出=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 onnyquist(n,d2) ; nyquist(n,d3) ; nyquist(n,d4) ;12. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為

30、 繪制系統(tǒng)的Nyqusit圖，并討論系統(tǒng)的穩(wěn)定性. G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis(square)13. 分別由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) ; %人工變量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=zp2

31、tf(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(num,den);impulse(num2,den2);p = -10.5969 +36.2148i -10.5969 -36.2148i -0.0562 15. 已知系統(tǒng)為：作該系統(tǒng)的nichols曲線。 n=1 ; d=1 , 1 , 0 ; nichols(n , d) ;16. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為：當(dāng)k=2時(shí)，分別作nichols曲線和波特圖。 num=1;de

32、n=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; %波特圖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)ans

33、 = -2.5214 -0.2393 + 0.8579i -0.2393 - 0.8579ians = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316i18. 系統(tǒng)的狀態(tài)方程為： 試確定系統(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)六 連續(xù)系統(tǒng)數(shù)字仿真的基本算法實(shí)驗(yàn)任務(wù) 1理解歐拉法和龍格-庫塔法的基本思想； 2理解數(shù)值積分算法的計(jì)算精度、速度、穩(wěn)定性與步長的關(guān)系；程序舉例1. 取h=0.2

34、，試分別用歐拉法、RK2法和RK4法求解微分方程的數(shù)值解，并比較計(jì)算精度。 注：解析解： cleart(1)=0; y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; h=0.001; % 步長修改為0.001for k=1:5 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1);endfor k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);endfor k=1:5 k1=y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k

35、1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;endfor 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*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;

36、end disp( 時(shí)間 解析解 歐拉法 RK2法 RK4法)yt=t, y, y_euler, y_rk2, y_rk4;disp(yt) 時(shí)間 解析解 歐拉法 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.0040 1.0040 1.0040 0.0050 1.0050 1.0050 1.0050 1.0050 2.

37、考慮如下二階系統(tǒng)： 在上的數(shù)字仿真解（已知：，），并將不同步長下的仿真結(jié)果與解析解進(jìn)行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法進(jìn)行計(jì)算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件 采用RK4法進(jìn)行計(jì)算。 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)/2).*cos(sqrt(3)/2*t(k+1)+100*sqr