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

1、.計(jì)算機(jī)仿真技術(shù)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)三利用數(shù)值積分算法的仿真實(shí)驗(yàn)精品文檔.實(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法及四階龍格庫(kù)塔法構(gòu)建系統(tǒng)仿真模型的方法 , 并對(duì)仿真結(jié)果進(jìn)行分析。二 實(shí)驗(yàn)內(nèi)容系統(tǒng)電路如圖2.1所示。電路元件參數(shù)：直流電壓源 E1V ，電阻 R 10，電感L 0.01H ，電容 C1F 。電路元件初始值： 電感電流 iL (0)0 A ，電容電壓 uc (0)0V 。系統(tǒng)輸出量為電容電壓 uc (t ) 。連續(xù)系統(tǒng)輸出

2、響應(yīng) uc (t ) 的解析解為：uc (t )U s (1e at(cos tsinta /)（2-1 ）R1R2其中， a。，LC2L2LRLiL (t)DCEuc (t )C圖 2.1 RLC 串聯(lián)電路三、要求1） 利用歐拉法、梯形法、二階顯式 Adams法及顯式四階 Runge-Kutta 法構(gòu)建系統(tǒng)仿真模型，并求出離散系統(tǒng)的輸出量響應(yīng)曲線；精品文檔.2） 對(duì)比分析利用歐拉法、梯形法、二階顯式 Adams法及顯式四階 Runge-Kutta 法構(gòu)建系統(tǒng)仿真模型的仿真精度與模型運(yùn)行的穩(wěn)定性問(wèn)題；3） 分別編寫歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta 法的 .m

3、 函數(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ù)系統(tǒng)輸出響應(yīng) uc (t) 的解析解為：uc (t)U s (1e at(costsintx /)其中， xR ，1R 22LLC2L（ 2）原系統(tǒng)的傳遞函數(shù)根據(jù)所示電路圖，我們利用電路原理建立系統(tǒng)的傳遞函數(shù)模型，根據(jù)系統(tǒng)的傳遞函數(shù)是在零初始條件下輸出量的拉普拉斯變換與輸入量的拉普

4、拉斯變換之比，可得該系統(tǒng)的傳遞函數(shù)：U C (s)1/ LCG(s)s2R / Ls 1/ LCE(s)(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)(t)x tAxBux (t0 ) x 0()()()y tCx tDu t其中：xx1 (t)x2 (t)xn (t) T 是系統(tǒng)的

5、 n 維狀態(tài)向量精品文檔.u(t)u1 (t )u2 (t)um (t ) T 是系統(tǒng)的 m維輸入向量y(t)y1 (t )y2 (t )yr (t ) T 是系統(tǒng)的 r 維輸出向量A 為 nn 階參數(shù)矩陣，又稱動(dòng)態(tài)矩陣， B 為 nm 階輸入矩陣， C 為 rn 階輸出矩陣， D為 rm 階交聯(lián)矩陣。根據(jù)圖所示電路，系統(tǒng)狀態(tài)方程模型 :x (t)Ax (t) BEy(t)Cx (t )式中，狀態(tài)變量 x x , x TiL,uCT，輸出變量 y(t)uC，系數(shù)矩陣為：1 2AR / L1/ L， B1/ L， C01 。1/ C00(1) 歐拉法利用前向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為：xm

6、1xm&1 Ah xmhBu mxmhym 1Cx m 1Du m 1式中， h 為積分步長(zhǎng)， 1 為單位矩陣。利用后向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為：&1xm1xm1h1 Ahxm hBu m 1xmym1Cxm 1Dum1對(duì)于前向歐拉法，系數(shù)矩陣為： A z 1hA ， B zhB ，C zC ，D=0。對(duì)于后向歐拉法，系數(shù)矩陣為： Az111 hA ， Bz1 hAhB ， C z C 。(2) 梯形法利用梯形法構(gòu)建線性系統(tǒng)的仿真模型為：精品文檔.xm 1xmh&&2xmxm 11xm 11h A1h A x mh B um um 1222ym 1Cxm1

7、Dum1對(duì)圖所示的系統(tǒng)，利用梯形法構(gòu)造的系統(tǒng)差分方程具有形式：xm 1At xm2Bt Eym 1Ct xm1其系數(shù)矩陣為： A1hA /121 hA / 2 ，t1hB / 2 ， CtBt1 hA / 2C，D=0。（3）二階顯式 Adams法利用二階顯式 Adams法構(gòu)建線性系統(tǒng)的仿真模型為：xm 1xmh23F m 16Fm 1 5Fm 212ym1Cx m 1Du m 1FmAx mBu m式中：Fm1Ax m 1Bu m1Fm2Ax m2Bu m2二階顯式 Adams法為多步計(jì)算方法，利用多步計(jì)算方法對(duì)系統(tǒng)進(jìn)行仿真時(shí)，需要與之具有相同計(jì)算精度的單步計(jì)算方法輔助計(jì)算。二階顯式 Ada

8、ms法的計(jì)算精度為二階，可以采用梯形法或改進(jìn)的 Euler 法等輔助計(jì)算。利用改進(jìn)的 Euler 法構(gòu)建線性系統(tǒng)的仿真模型為：k1f ( tm, xm )Ax mBEk2f ( tmh,xhk)A xmk hBEm11xm 1xmh ( k1k 2 )2其中， m 0,1 。由式計(jì)算出 x1和 x2 后，便可以轉(zhuǎn)入由二階顯式Adams法構(gòu)造的離散系統(tǒng)模型計(jì)算，即系統(tǒng)差分方程。其計(jì)算方程為：精品文檔.F mAx m BEF m 1Ax m 1B EF m2Ax m 2BExm 1xmh23F m 16Fm 15Fm212（ m1 ）ym 1Cx m 1（4）顯式四階 Runge-Kutta 法利

9、用顯式四階 Runge-Kutta 法構(gòu)建線性系統(tǒng)的仿真模型為：k1f ( tm, x m)Axm BEk 2f( tmh , xmh k1 )A xmk1h / 2BE22k3f ( tmh , xmh k 2)A xmk 2h / 2BE22k 4f( tmh, xmhk3)A xmk3 hBExm 1x mh (k122k4)6k2k 3ym 1Cxm1五實(shí)驗(yàn)過(guò)程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;h = 2.0e-4;m = fix(t/h);n = 2;A = -R/L -

10、1/L;1/C 0;B = 1/L;0;D=01;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(w*(k-1)*h) + sin(w*(k-1)*h)*p/w); end%輸

11、出曲線 %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/L -1/L;1/C 0;B = 1/L;0;D=01;E = 1 0;0 1;% 后向

12、歐拉法 %for i=1:1:nx2(1:n,1) = 0;endA1 = inv(E-A*h);for k=1:mx2(1:n,k+1) = A1*(x2(1:n,k) + B*h);精品文檔.endfor k=1:1:my2(k) = D*x2(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(w*(k-1)*h) + sin(w*(k-1)*h)*p/w);end%輸出曲線 %for k=1:1:mt(k) = (k-1)*h;endsu

13、bplot(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/C 0;B = 1/L;0;D=01;E = 1 0;0 1;%梯形法 %for i=1:1:nx3(1:n,1) = 0;endA2 = inv

14、(E-A*h/2);for k=1:mx3(1:n,k+1) = A2*( x3(1:n,k) + B*h + A*x3(1:n,k)*h/2); endfor k=1:1:my3(k) = D*x3(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(w*(k-1)*h)+sin(w*(k-1)*h)*p/w);end%輸出曲線 %for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,3),plot(t,y,'

15、;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;A = -R/L -1/L;1/C 0;B = 1/L;0;D=01;E = 1 0;0 1;% 二階顯示 Adams法 %for i=1:1:nx4(1:n,1) = 0;endfor k=1:mx4(1:n,k+1)

16、 = A2*(x4(1:n,k) + B*h + A*x4(1:n,k)*h/2); endfor k=3:mfm1 = 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)*h/12;endfor k=1:1:my4(k) = D*x4(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)*

17、h)*(cos(w*(k-1)*h)+sin(w*(k-1)*h)*p/w);end%輸出曲線 %for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,4),plot(t,y,'g',t,y4,'r')legend('y解析解 ,','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

18、/h);n = 2;A = -R/L -1/L;1/C 0;B = 1/L;0;D=01;E = 1 0;0 1;% 四階 Runge-Kutta 法 %for i=1:1:n %狀態(tài)變量初值x5(1:n,1) = 0;endfor k=1:mx5(1:n,k+1) = A2*( x5(1:n,k) + B*h + A*x5(1:n,k)*h/2); endfor k=1:1:mk1=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

19、:n,k+1)+h.*(k1+2*k2+2*k3+k4)./6;endfor 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:my(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:mt(k) = (k-1)*h;endsubplot(2,3,5),plot(t,y,'g',t,y5,'r')legend(

20、9;y解析解 ,','y5Runge-Kutta法 ')title('顯式四階 Runge-Kutta 法')2. 仿真圖形取積分步長(zhǎng) h=2*10-4 s，可以得到以下幾個(gè)仿真圖形：(1) 前向歐拉法x 101y 解析解,y1前向歐拉0-1-2-3-4-500.0010.0020.0030.0040.0050.0060.0070.0080.0090.01（ 2）后向歐拉法精品文檔.1.6后向歐拉法y解 析 解 ,1.4y2后 向 歐 拉1.210.80.60.40.200.0010.0020.0030.0040.0050.0060.0070.0080.

21、0090.010（ 3）梯形法梯 形 法21.8y解 析 解 ,y3 梯 形 法1.61.41.210.80.60.40.2000.0010.0020.0030.0040.0050.0060.0070.0080.0090.01（ 4）二階顯式 Adams法x 1026二階顯式Adams法0.5y 解析解 ,y4Adams法016前向歐拉法x 101-0.50y 解析解,后向歐拉法y1前向歐拉2y 解,梯形法析解2-1y2 后向歐拉y 解,1.5析解-1-2y3 梯形法1.5-31-1.5-41-50.5-200.0050.010.500.0050.010000.0050.01-2.5-300.0010.0020.0030.0040.0050.0060.0070.0080.0090.01顯式2（ 5）四階 Runge-Kutta 法精品文檔.2顯 式四階 Runge-Kutta法y 解 析 解 ,y5 龍 格庫(kù) 塔1.516前 向 歐拉 法x 10y 解 析 解 ,10y1 前 向歐 拉后 向歐 拉 法2梯 形