計(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法及四階龍格庫塔法構(gòu)建系統(tǒng)仿 真模型的方法，并對(duì)仿真結(jié)果進(jìn)行分析。二. 實(shí)驗(yàn)容系統(tǒng)電路如圖2.1所示。電路元件參數(shù)：直流電壓源E = W,電阻R = 10G,電感L = 0.01H,電容C = l/zF。電路元件初始值：電感電流z(O) = OA,電容電壓c(0) = 0V 0 系統(tǒng)輸出量為電容電壓i/c(r)o連續(xù)系統(tǒng)輸出響應(yīng)的解析解為：(2

2、-1)uc(t) = Ux(cosof + sin 必xa/勁)其中,Ra =2L1R LC69 =9RLM)dcO圖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)定性問題；3) 分別編寫歐拉法.梯形法、二階顯式Adams法及顯式四階Runge-Kutta法的m函數(shù)文 件，并存入磁盤中。ni函數(shù)文件要求輸入?yún)?shù)為系統(tǒng)狀態(tài)方程的系數(shù)矩陣、仿真時(shí)間及仿 真

3、步長(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)/1(0 ：(0兒是系統(tǒng)的廠維輸出向量A為nxn階參數(shù)矩陣，又稱動(dòng)態(tài)矩陣，為n x m階輸入矩陣，6*為/公階輸出矩陣，D 為rxm階交聯(lián)矩陣。根據(jù)圖所示電路，系統(tǒng)狀態(tài)方程模型：x(t) = Ax(t) + BEy(t) = Cx(t)式中，狀態(tài)變量x = xhx2r =iL,ucr.輸出變量y(t) = iic,系數(shù)矩陣為:-/?/L-1/L 1/LA =,B =1/C00,

4、C = o lo(1) 歐拉法利用前向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為:卜”+i =兀+ % = ( + Ah)?！?+ hB%1九+】=Cxt+i +式中，力為積分步長(zhǎng)，2為單位矩陣。利用后向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為: 1 )h/n+1+ 12(4) 顯式四階Runge-Kutta法利用顯式四階Runge-Kutta法構(gòu)建線性系統(tǒng)的仿真模型為:K =九,心)=化+ bek2 =+ ，?！?+ 人)=A (x)n + kji / 2)+ BEk3 = fg + ，心 + &)= A (xn + k2h / 2)+ BE % = f(tm + h, xm + hk3) = A (xm + k

5、3h) + BE=札+ 匕+汰2 + 2他+心)五.實(shí)驗(yàn)過程1. 實(shí)驗(yàn)程序(1) 前向歐拉法function =RLC(R,L,C,U,t,h)R=10；L=0.01；C=l. 0e6；U=l；t=0.01；h = 2. Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1;E = 1 0；0 1;%前向歐拉法%for i=l:1:nxl (l：n,1) = 0；endfor k=l：mxl (l：n,k+l) = xl (l：n,k) + (A* xl (1：n,k)+B)*h；endfor k=l:l：myl(k) =

6、D*xl(l：n,k)；end%解析解%P = R/(2*L)；w=sqrt(1/(L*C)-(R/(2*L)*2);for k=l:l：my (k) = U*(l-exp(-p*(k-l) *h) * ( cos (w*(kl) *h) + sin(w*(kl)*h) *p/w);end%輸出曲線%for k=l:l：mt (k) = (k-1) *h；endsubplot (2,3,1) ,plot (t,y,legend(y解析解/ , ryl前向歐拉)titleC前向歐拉法)(2) 后向歐拉法function =RLC(R,L,C,U,t,h)R=10；L二0.01；C=l. Oe-6

7、；U=l；t 二0 01；h = 2.Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E = 1 0；0 1;%后向歐拉法%for i=l：l：nx2(l：n, 1) = 0；endAl = inv(EA*h)；for k=l：mx2(l：n,k+l) = Al*(x2(l：n,k) + B*h);endfor k=l：1：my2(k) = D*x2(l：n,k);end%解析解%P = R/ (2*L);w=sqrt(1/(L*C)-(R/(2*L)2)；for k=l：1：my (k)= LI* (1-exp (-

8、p* (k-1) *h)*( cos (w*(kl) *h)+sin(w*(k-1)*h)*p/w)；end%輸出曲線%for k=l:1：mt (k) = (k-1) *h；endsubplot (2,3,2),plot(t,y,* g*,t,y2,r*)legend( y解析解，y2后向歐拉) titleC后向歐拉法)(3) 梯形法function =RLC(R,L,C,U,t,h)R=10；L=0.01；C=l. 0e6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h);n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E

9、= 1 0；0 1;%梯形法%for i=l:l：nx3(l：n, 1) = 0；endA2 = inv(E-A*h/2)；for k=l：mx3(l：n,k+l) = A2*( x3(l：n,k) + B*h + A*x3(l：n,k)*h/2)；endfor k=l：1：my3(k) = D*x3(1：n,k)；end%解析解%P = R/ (2*L);w=sqrt (1/ (L*C)-(R/(2*L) 2)；for k=l：l：my (k)= LI* (1-exp (-p* (k-1) *h)*( cos (w*(k-l) *h)+sin(w*(kl)*h)*p/w)；end%輸出曲線%

10、for k=l：1：mt (k) = (k-1) *h；endsubplot (2,3,3),plot (t,y,1gf,t,y3,r) legend( y解析解，r, r y3梯形法) title(梯形法)(4) 二階顯式Adams法function =RLC(R,L,C,U,t,h)R=10；0.01；C=l. Oe-6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E = 1 0；0 1;%二階顯示Adams法%for i=l：1:nx4(l：n, 1) = 0；endf

11、or k=l：mx4(l：n,k+l) = A2*(x4(l：n,k) + B*h + A*x4(l:n,k)*h/2):endfor k=3：mfml = 23*(A*x4(l：n,k) + B)；fm2 = -16*(A*x4(l：n,k-l)+ B)；fm3 = 5*(A*x4(l：n,k-2)+ B)；x4(l：n,k+l) = x4(1：n,k)+(fml+fm2+fm3)*h/12；endfor k=l：1：my4 (k) = D*x4(l：n,k);end%解析解%P = R/ (2*L);w=sqrt (1/ (L*C)-(R/(2*L) 2)；for k=l：l：my (k)

12、= LI* (1-exp (-p* (k-1) *h)*( cos (w*(k-l) *h)+sin(w*(kl)*h)*p/w)；end%輸出曲線%for k=l：1：mt (k) = (k-1) *h；endsubplot(2,3,4),plot(t,y,g，t,y4,i)legend ( y 解析解，r,r y4Adams 法 J titleC 二階顯式 Adams 法)(5) 四階 Runge-Kutta 法function =RLC(R,L,C,U,t,h)R=10；0.01；C=l. Oe-6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h)；n = 2；

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

14、+h. *(kl+2*k2+2*k3+k4). /6；endfor k=l：1:my5(k) = D*x5(l：n,k)；end%解析解%P = R/ (2*L);w=sqrt(1/(L*C)-(R/(2*L)2)；for k=l：l：my(k)=U*(1-exp(-p*(k-1)*h)*( cos(w*(k-l)*h)sin (w*(k-1)*h)*p/w)；end%輸出曲線%for k=l：1：mt (k) = (k-1) *h；endsubplot (2,3,5),plot(t,y,g*,t,y5,r) legend ( y 解析解，/ y5RungeKutta 法) title(顯式四

15、階 Runge-Kutta 法)2. 仿真圖形取積分步長(zhǎng)h=2*10 4s,可以得到以下幾個(gè)仿真圖形：(1)前向歐拉法11 1y解儕昶.”仍同吹垃1 1 1 1 1 1r11rrrr11OO.OO10.0020.003O.OCM0.0050.0060.0070.0080.0090.01(2)后向歐拉法y解析解y馮向歐拉0.0010.0020.0030.0040.0050.0060.0070.0080.0090.01(3)梯形法梯形法(4) 二階顯式Adams法x 10二陸 SB re AdamsiAl11ll11lyAdams；/*riirriirrOO.OOl0.0020.0030.0040.0050.0060.0070.0080.009O.Ol(5) 四階 Runge-Kutta 法六.實(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-Kutta法建模最為復(fù)雜，仿真時(shí)間也較長(zhǎng)。3