計(jì)算機(jī)仿真技術(shù)實(shí)驗(yàn)報(bào)告-實(shí)驗(yàn)三

1、仿真技術(shù)與應(yīng)用實(shí)驗(yàn)報(bào)告計(jì)算機(jī)仿真技術(shù)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)三 利用數(shù)值積分算法的仿真實(shí)驗(yàn) - 13 -實(shí)驗(yàn)三 利用數(shù)值積分算法的仿真實(shí)驗(yàn)一 實(shí)驗(yàn)?zāi)康?) 熟悉MATLAB的工作環(huán)境；2) 掌握MATLAB的 .M文件編寫規(guī)則，并在命令窗口調(diào)試和運(yùn)行程序；3) 掌握利用歐拉法、梯形法、二階顯式Adams法及四階龍格庫塔法構(gòu)建系統(tǒng)仿真模型的方法,并對(duì)仿真結(jié)果進(jìn)行分析。二 實(shí)驗(yàn)內(nèi)容系統(tǒng)電路如圖2.1所示。電路元件參數(shù)：直流電壓源，電阻，電感，電容。電路元件初始值：電感電流，電容電壓。系統(tǒng)輸出量為電容電壓。連續(xù)系統(tǒng)輸出響應(yīng)的解析解為： （2-1）其中，。三、要求1）利用歐拉法、梯形法、二階顯式Adams法及顯式

2、四階Runge-Kutta法構(gòu)建系統(tǒng)仿真模型，并求出離散系統(tǒng)的輸出量響應(yīng)曲線；2）對(duì)比分析利用歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法構(gòu)建系統(tǒng)仿真模型的仿真精度與模型運(yùn)行的穩(wěn)定性問題；3）分別編寫歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法的.m函數(shù)文件，并存入磁盤中。.m函數(shù)文件要求輸入?yún)?shù)為系統(tǒng)狀態(tài)方程的系數(shù)矩陣、仿真時(shí)間及仿真步長(zhǎng)。編寫.m命令文件，在該命令文件中調(diào)用已經(jīng)編寫完成的上述.m函數(shù)文件，完成仿真實(shí)驗(yàn)；4） subplot和plot函數(shù)將輸出結(jié)果畫在同一個(gè)窗口中，每個(gè)子圖加上對(duì)應(yīng)的標(biāo)題。四.實(shí)驗(yàn)原理（1）連續(xù)系統(tǒng)解析解連續(xù)

3、系統(tǒng)輸出響應(yīng)的解析解為： 其中，（2）原系統(tǒng)的傳遞函數(shù) 根據(jù)所示電路圖，我們利用電路原理建立系統(tǒng)的傳遞函數(shù)模型，根據(jù)系統(tǒng)的傳遞函數(shù)是在零初始條件下輸出量的拉普拉斯變換與輸入量的拉普拉斯變換之比，可得該系統(tǒng)的傳遞函數(shù)：(3)系統(tǒng)的仿真模型在連續(xù)系統(tǒng)的數(shù)字仿真算法中,較常用的有歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法等。歐拉法、梯形法和二階顯式Adams法是利用離散相似原理構(gòu)造的仿真算法，而顯式四階Runge-Kutta法是利用Taylor級(jí)數(shù)匹配原理構(gòu)造的仿真算法。對(duì)于線性系統(tǒng)，其狀態(tài)方程表達(dá)式為: 其中：是系統(tǒng)的n維狀態(tài)向量是系統(tǒng)的m維輸入向量是系統(tǒng)的r維輸出向量

4、A為階參數(shù)矩陣，又稱動(dòng)態(tài)矩陣，B為階輸入矩陣，C為階輸出矩陣，D為階交聯(lián)矩陣。根據(jù)圖所示電路，系統(tǒng)狀態(tài)方程模型: 式中，狀態(tài)變量，輸出變量，系數(shù)矩陣為：，。(1) 歐拉法利用前向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為： 式中，為積分步長(zhǎng)，為單位矩陣。利用后向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為： 對(duì)于前向歐拉法，系數(shù)矩陣為：，D=0。對(duì)于后向歐拉法，系數(shù)矩陣為：，。(2) 梯形法利用梯形法構(gòu)建線性系統(tǒng)的仿真模型為： 對(duì)圖所示的系統(tǒng)，利用梯形法構(gòu)造的系統(tǒng)差分方程具有形式： 其系數(shù)矩陣為：，D = 0。（3）二階顯式Adams法利用二階顯式Adams法構(gòu)建線性系統(tǒng)的仿真模型為： 式中： 二階顯式Adams法為多

5、步計(jì)算方法，利用多步計(jì)算方法對(duì)系統(tǒng)進(jìn)行仿真時(shí)，需要與之具有相同計(jì)算精度的單步計(jì)算方法輔助計(jì)算。二階顯式Adams法的計(jì)算精度為二階，可以采用梯形法或改進(jìn)的Euler法等輔助計(jì)算。利用改進(jìn)的Euler法構(gòu)建線性系統(tǒng)的仿真模型為： 其中，。由式計(jì)算出和后，便可以轉(zhuǎn)入由二階顯式Adams法構(gòu)造的離散系統(tǒng)模型計(jì)算，即系統(tǒng)差分方程。其計(jì)算方程為： （）（4）顯式四階Runge-Kutta法 利用顯式四階Runge-Kutta法構(gòu)建線性系統(tǒng)的仿真模型為： 五實(shí)驗(yàn)過程1.實(shí)驗(yàn)程序（1）前向歐拉法function =RLC(R,L,C,U,t,h)R=10;L=0.01;C=1.0e-6;U=1;t=0.01

6、;h = 2.0e-4;m = fix(t/h);n = 2;A = -R/L -1/L;1/C 0;B = 1/L;0;D = 0 1;E = 1 0;0 1;% 前向歐拉法 %for i=1:1:nx1(1:n,1) = 0;endfor k=1:mx1(1:n,k+1) = x1(1:n,k) + (A* x1(1:n,k)+B)*h;endfor k=1:1:my1(k) = D*x1(1:n,k);end% 解析解 %p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)2);for k=1:1:my(k) = U*(1-exp(-p*(k-1)*h) * ( cos

7、(w*(k-1)*h) + sin(w*(k-1)*h)*p/w);end%輸出曲線 %for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,1),plot(t,y,'g',t,y1,'r')legend('y解析解,','y1前向歐拉')title('前向歐拉法')（2）后向歐拉法function =RLC(R,L,C,U,t,h) R=10; L=0.01; C=1.0e-6; U=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2; A = -R/

8、L -1/L;1/C 0; B = 1/L;0; D = 0 1; E = 1 0;0 1;% 后向歐拉法 % for i=1:1:n x2(1:n,1) = 0; end A1 = inv(E-A*h); for k=1:m x2(1:n,k+1) = A1*(x2(1:n,k) + B*h); end for k=1:1:m y2(k) = D*x2(1:n,k); end % 解析解 % p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)2); for k=1:1:m y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) +

9、sin(w*(k-1)*h)*p/w); end %輸出曲線 % for k=1:1:m t(k) = (k-1)*h; end subplot(2,3,2),plot(t,y,'g',t,y2,'r') legend('y解析解,','y2后向歐拉') title('后向歐拉法')（3）梯形法function =RLC(R,L,C,U,t,h) R=10; L=0.01; C=1.0e-6; U=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2; A = -R/L -1/L;1

10、/C 0; B = 1/L;0; D = 0 1; E = 1 0;0 1; % 梯形法 % for i=1:1:n x3(1:n,1) = 0; end A2 = inv(E-A*h/2); for k=1:m x3(1:n,k+1) = A2*( x3(1:n,k) + B*h + A*x3(1:n,k)*h/2); end for k=1:1:m y3(k) = D*x3(1:n,k); end % 解析解 % p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)2); for k=1:1:m y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w

11、*(k-1)*h) + sin(w*(k-1)*h)*p/w); end %輸出曲線 % for k=1:1:m t(k) = (k-1)*h; end subplot(2,3,3),plot(t,y,'g',t,y3,'r') legend('y解析解,','y3梯形法') title('梯形法')（4）二階顯式Adams法function =RLC(R,L,C,U,t,h) R=10; L=0.01; C=1.0e-6; U=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2;

12、 A = -R/L -1/L;1/C 0; B = 1/L;0; D = 0 1; E = 1 0;0 1;% 二階顯示Adams法 % for i=1:1:n x4(1:n,1) = 0; end for k=1:m x4(1:n,k+1) = A2*(x4(1:n,k) + B*h + A*x4(1:n,k)*h/2); end for k=3:m fm1 = 23*(A*x4(1:n,k)+ B); fm2 = -16*(A*x4(1:n,k-1)+ B); fm3 = 5*(A*x4(1:n,k-2)+ B); x4(1:n,k+1) = x4(1:n,k)+(fm1+fm2+fm3)

13、*h/12; end for k=1:1:m y4(k) = D*x4(1:n,k); end % 解析解 % p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)2); for k=1:1:m y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w); end %輸出曲線 % for k=1:1:m t(k) = (k-1)*h; end subplot(2,3,4),plot(t,y,'g',t,y4,'r') legend('y解析解,',&

14、#39;y4Adams法') title('二階顯式Adams法')（5）四階Runge-Kutta法function =RLC(R,L,C,U,t,h) R=10; L=0.01; C=1.0e-6; U=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2; A = -R/L -1/L;1/C 0; B = 1/L;0; D = 0 1; E = 1 0;0 1;% 四階Runge-Kutta法 % for i=1:1:n % 狀態(tài)變量初值 x5(1:n,1) = 0; end for k=1:m x5(1:n,k+1) = A2*(

15、 x5(1:n,k) + B*h + A*x5(1:n,k)*h/2); end for k=1:1:m k1=A*x5(1:n,k+1); k2=A*(x5(1:n,k+1)+h*k1/2); k3=A*(x5(1:n,k+1)+h*k2/2); k4=A*(x5(1:n,k+1)+h*k3); x5(1:n,k+1)=x5(1:n,k+1)+h.*(k1+2*k2+2*k3+k4)./6; end for k=1:1:m y5(k) = D*x5(1:n,k); end % 解析解 % p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)2); for k=1:1:m

16、y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w); end %輸出曲線 % for k=1:1:m t(k) = (k-1)*h; end subplot(2,3,5),plot(t,y,'g',t,y5,'r') legend('y解析解,','y5Runge-Kutta法 ') title('顯式四階Runge-Kutta法')2.仿真圖形取積分步長(zhǎng)h=2*10-4s，可以得到以下幾個(gè)仿真圖形：(1)前向歐拉法（2）后向歐拉法（3）梯形法（4）二階顯式Adams法（5）四階Runge-Kutta法6 實(shí)驗(yàn)結(jié)論 1.從仿真的穩(wěn)定性看，當(dāng)選取不同的積分步長(zhǎng)時(shí)，歐拉法穩(wěn)定性最低，梯形法穩(wěn)定性其次，而顯式四階Runge-Kutta法、二階顯示Adams法穩(wěn)定性較好。 2.從仿真的難易性看，歐拉法為單步計(jì)算法，用到一個(gè)過去的值，計(jì)算起來比較簡(jiǎn)單。而梯形法則是用兩條折線所謂面積來近似，與歐拉法相比較為困難。二階顯示Adams法需要知道k個(gè)初始值，不能自起步，二次函數(shù)很復(fù)雜，因此此方法較復(fù)雜。而顯式四階Runge-Kutt