1、中南大學系統仿真實驗報告 ans = 實驗一 matlab 中矩陣與多項式的基本運算 實驗任務 1. 了解 matlab 命令窗口和程序文件的調用。 2 熟悉如下 matlab 的基本運算： 矩陣的產生、數據的輸入、相關元素的顯示； 矩陣的加法、乘法、左除、右除； 特殊矩陣：單位矩陣、 1 矩陣、 0 矩陣、對角陣、隨機矩陣的產生和 運算； 多項式的運算：多項式求根、多項式之間的乘除。 基本命令訓練 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) 1 1 ans = 1 1 ones

2、(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 4、隨機陣 rand(2,3) a

3、ns = 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 的逆矩陣) b=5 3 1;2 3 8;1 1

4、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 0.2581 ans = -2.75

5、00 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 4 0 -2 1 2 ans =

6、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、（ 點乘是數組的運算，沒有點的乘是矩陣運算 ） a = 2

7、 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 2x2 32 double a 1x2

8、 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 實驗二 matlab 繪圖命令 實驗任務 熟悉 matlab 基本繪圖命令，掌握如下繪圖方法： 1 坐標系的選擇、圖形的繪制； 2 圖形注解（題目、標號、說明、分格線）的加入; 3 圖形線型、符號、顏色的選取 基本命令訓練 1、t=0

9、: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); y=si n( t)+ si

10、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); c,h = con tour(x,

11、 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 的原函數); subplot(412) ;semilogy(x,y); grid;title(對 y 取對數); 丫的原畫數 40 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - l| i p | il | i| i 九 _ 1- _ _ i _ l _ : ： _ j _ _ l _ u i| i |l

12、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 取對數 對弋觀對數 subplot(413) ;semilogx(x,y); 10 10 id 10 1lz 10 w 10 10 40 20 0 對好對數 2 掌握循環、分支語句的編寫，學會使用 look for 、 help 命令 grid;title(對 x 取對數); subplot(414) ;loglog(x,y);grid; title(對 xy 均取對數); 9、x = -3:0.3:3; bar(x,e

13、xp(-x.*x),g) 實驗三 matlab 程序設計 實驗任務 1 熟悉 matlab 程序設計的方法和思路； 程序舉例 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-1); end end forma

14、t rat 1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5 (分數格式形式。用有理數逼近顯示數據 ) 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 命令時，如果上次運行結果沒有清除，輸出的結果就是上次運行的 結果！但是運用 clear 命令清楚之前的運行結果之后就會正常運行。 4、 x

15、=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 ( 不在定義域內，請重新輸入： );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(2:length(p); end

16、; 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 計算斐波那契亞數列的函數文件 help e2 ffibno 計算斐波那契亞數列的函數文件 n 可取任意自然數 程序如下 （用法： lookfor 關鍵詞 在所有 m 文件中找 關鍵詞 ，比如： lookfor

17、 max （即尋找關鍵詞 max） 其實就和我們平時用 ctrl+f 來查找關鍵詞是一樣的 而 help 是顯示 matlab 內置的幫助信息 用法： help 命令，比如 help inv ，作用就是調用 inv 這個命令的幫助） 程序設計題 用一個 matlab 語言編寫一個程序：輸入一個自然數，判斷它是否是素數， 如果是，輸出 it is one prime ，如果不是，輸出 it is n ot o ne prime. 。要求通 過調用子函數實現。 最好能具有如下功能： 設計較好的人機對話界面， 程序中 含有提示 性的輸入輸出語句。 能實現循環操作， 由操輸入相關命令來控制 是否繼續進

18、行素數的判斷。 如果操希望停止這種判斷， 則可以退出程序。 如果所輸入的自然數是一個合數， 除了給出其不是素數的結論外， 還應給出至少 一種其因數分解形式。例：輸入 6 ，因為 6 不是素數。則程序中除了有 it is not one prime 的結論外，還應有： 6=2*3 的說明。 function sushu while 1 x=input（ 請輸入一個自然數 ）; if x2 disp（ 既不是質數又不是合數 ）; else if isprime(x)=1 disp（ 這是一個素數 ）; else disp( 這是一個合數，可以因式分解為： ) for n=2:sqrt(x) if

19、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( 是否繼續判斷？繼續請按 1 ,按任意鍵退出 ： ) ； if y=1 break end end 實驗四 matlab 的符號計算與 simulink 的使用 實驗任務 1 .掌握 matlab 符號計算的特點和常用基本命令; 2 .掌握 simulink 的使用。 程序舉例 1 . 求矩陣對應的行列式和特征根 a=sym(a11 a12;a21 a22); da=det

20、(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=solve(xa2+x-1) rv=

21、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 個符號變量 %定義 4 個數值變量 %建立符號矩陣

22、 a %建立數值矩陣 b %計算符號矩陣 a 的行列式 %計算數值矩陣 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 = (8*ya2+7*xa2)*(x

23、a2-3*ya2) 5. 對方程 ax=b 求解 a=34,8,4;3,34,3;3,6,8; b=4;6;2; x=linsolve(a,b) % 調用 linsolve 函數求解 ab % 用另一種方法求解 0.0675 0.1614 diff(f) %未指定求導變量和階數，按缺省規則處理 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 a21 a22 a23 a31 a32

24、a33 b1 b2 b3; a=a11,a12,a13;a21,a22,a23;a31,a32,a33; b=b1;b2;b3; xx=ab %用左除運算求解 (x=linsolve(a,b) %調用 linsolve 函數求的解 ) 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-a11*b2*a33-a21*a13*b

25、3-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(1+exp(x); %求 f 對

26、x 的二階導數 %求 f 對 x 的三階導數 %按參數方程求導公式求 y 對 x 的導數 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) %未指定求導變量和階數，按缺省規則處理 仿真圖: 波形圖: 其中: g i

27、( r 0.14 實驗五 matlab 在控制系統分析中的應用 實驗任務 1. 掌握 matlab 在控制系統時間響應分析中的應用; 2. 掌握 matlab 在系統根軌跡分析中的應用； 3. 掌握 matlab 控制系統頻率分析中的應用； 4. 掌握 matlab 在控制系統穩定性分析中的應用 基本命令 1. step 2. impulse 3. in itial 4. isim 5. rlocfi nd 6. bode 7. margin 8. nyquist 9. nichols 10. cloop 程序舉例 1. 求下面系統的單位階躍響應 1.5 g(s) t num=4 ; den=

28、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. 求如下系統的單位階躍響應 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 x 2 3. 求下面系統的單位脈沖

29、響應: g(s) 兒 num=4 ; den=1 , 1 ,4; impulse( nu m,de n) time (sec) response to initial conditions 4. 已知二階系統的狀態方程為: 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); impulse(a , b , c , d)

30、 5 :系統傳遞函數為: 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. 有一二階系統，求出周期為 4 秒的方波的輸出響應 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 函數功能 lsim( nu m,de n,u,t);

31、7. 已知開環系統傳遞函數，繪制系統的根軌跡，并分析其穩定性 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 results1.5 e 0.5 -0

32、.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 8 ! - impulse resp

33、 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.0248i 0.7133 -3.41

34、60 -2.5918 -0.9961 + 0.4306i -0.9961 - 0.4306i bode diagram 8. 作如下系統的 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. 系統傳函如下 g(s) s 1 0.5s e (s 2) 3 求有理傳函的頻率響應，然后在同一張圖上繪出以四階伯德近似表示的系統頻 率響應 num=1;de n=co nv(1 2,co nv(1 2,1 2); -20 w=logspace(-1,2);

35、t=0.5; ylabel(gai n); subplot(2,1,2); semilogx(w,p1,w,p2,g-);grid on; xlabel(freque ncy); ylabel(phase); 10. 已知系統模型為 求它的幅值裕度和相角裕度 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(n2,nu m); den t=(

36、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. 二階系統為: g(s) n s 2 2 n s s 令 wn=1, 分別作出 e =2,1 , 0.707 , c a 5 y a 時的 nyquist 曲線。 m n=1; d 仁 1

37、, 4,1 ; d2=1 , 2 , 1 ； d3=1 , 1.414,1; d4=1,1,1; nyquist diagram real axis 12. 已知系統的開環傳遞函數為 -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. 一多環系統，其結構圖如下，使用 nyquist 頻率曲線判斷系統的穩定性。 16.7s (0.85s 1)(0.25s 1)(0.0625s 1)繪制系統的 nyqusit 圖, 并討論系統的穩定性 . g=tf(1000,co nv(1,3,2,1,5)

38、; nyquist(g);axis(square) 13. 分別由 w 的自動變量和人工變量作下列系統的 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 an nt t9 9a a 卩 g

39、(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.5 1 0.5 s a 0 n g

40、 m -0.5 -1 -1.5 1.5 20 impulse resp onse -10.5969 +36.2148i -10.5969 -36.2148i -0.0562 15. 已知系統為: 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) 作該系統的 nichols 曲線。 n=1 ; d=1 , 1 , 0; ni chols( n , d); 16. 已知系

41、統的開環傳遞函數為: g(s) k s(s 1)(s 2) 當 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 bode diagram 90 - -

42、 % 波特圖 17. 系統的開環傳遞函數為: 分別確定 k=2 和 k=10 時閉環系統的穩定性 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. 系統的狀態方程為： x 1 4 3 0 x 1 1 x

43、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) 試確定系統的穩定性。 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 實驗六連續系統數字仿真的基本算法 實驗任務 1 .理解歐拉法和龍格 - 庫塔法的基本思想； 2 理解數值積分算法的計算精度、速度、穩定性與步長的關系; 程序舉例 1. 取 h=0.2 ,試分別用歐拉法、 rk2 法和 rk4

44、 法求解微分方程的數值解，并 比較計算精度。 注：解析解：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 k1= y_rk2(k)-2*t(

45、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)+h*k3); y_rk4(k+1)=

46、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.0040 1.0040 1.0040 0

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