計算機仿真技術(shù)實驗報告-實驗三

1、仿真技術(shù)與應(yīng)用實驗報告計算機仿真技術(shù)實驗報告實驗三 利用數(shù)值積分算法的仿真實驗 - 13 -實驗三 利用數(shù)值積分算法的仿真實驗一 實驗目的1) 熟悉MATLAB的工作環(huán)境；2) 掌握MATLAB的 .M文件編寫規(guī)則，并在命令窗口調(diào)試和運行程序；3) 掌握利用歐拉法、梯形法、二階顯式Adams法及四階龍格庫塔法構(gòu)建系統(tǒng)仿真模型的方法,并對仿真結(jié)果進行分析。二 實驗內(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）對比分析利用歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法構(gòu)建系統(tǒng)仿真模型的仿真精度與模型運行的穩(wěn)定性問題；3）分別編寫歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法的.m函數(shù)文件，并存入磁盤中。.m函數(shù)文件要求輸入?yún)?shù)為系統(tǒng)狀態(tài)方程的系數(shù)矩陣、仿真時間及仿真步長。編寫.m命令文件，在該命令文件中調(diào)用已經(jīng)編寫完成的上述.m函數(shù)文件，完成仿真實驗；4） subplot和plot函數(shù)將輸出結(jié)果畫在同一個窗口中，每個子圖加上對應(yīng)的標題。四.實驗原理（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級數(shù)匹配原理構(gòu)造的仿真算法。對于線性系統(tǒng)，其狀態(tài)方程表達式為: 其中：是系統(tǒng)的n維狀態(tài)向量是系統(tǒng)的m維輸入向量是系統(tǒng)的r維輸出向量

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

5、步計算方法，利用多步計算方法對系統(tǒng)進行仿真時，需要與之具有相同計算精度的單步計算方法輔助計算。二階顯式Adams法的計算精度為二階，可以采用梯形法或改進的Euler法等輔助計算。利用改進的Euler法構(gòu)建線性系統(tǒng)的仿真模型為： 其中，。由式計算出和后，便可以轉(zhuǎn)入由二階顯式Adams法構(gòu)造的離散系統(tǒng)模型計算，即系統(tǒng)差分方程。其計算方程為： （）（4）顯式四階Runge-Kutta法 利用顯式四階Runge-Kutta法構(gòu)建線性系統(tǒng)的仿真模型為： 五實驗過程1.實驗程序（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.仿真圖形取積分步長h=2*10-4s，可以得到以下幾個仿真圖形：(1)前向歐拉法（2）后向歐拉法（3）梯形法（4）二階顯式Adams法（5）四階Runge-Kutta法6 實驗結(jié)論 1.從仿真的穩(wěn)定性看，當選取不同的積分步長時，歐拉法穩(wěn)定性最低，梯形法穩(wěn)定性其次，而顯式四階Runge-Kutta法、二階顯示Adams法穩(wěn)定性較好。 2.從仿真的難易性看，歐拉法為單步計算法，用到一個過去的值，計算起來比較簡單。而梯形法則是用兩條折線所謂面積來近似，與歐拉法相比較為困難。二階顯示Adams法需要知道k個初始值，不能自起步，二次函數(shù)很復雜，因此此方法較復雜。而顯式四階Runge-Kutt