實驗三哈工大計算機仿真.docx

實驗3 利用數(shù)值積分算法的仿真實驗(一 實驗?zāi)康模?）熟悉MATLAB的工作環(huán)境；（2）掌握MATLAB的 .M文件編寫規(guī)則，并在命令窗口調(diào)試和運行程序；（3）掌握利用歐拉法、梯形法、二階顯式Adams法及四階龍格庫塔法構(gòu)建系統(tǒng)仿真模型的方法,并對仿真結(jié)果進行分析。二、實驗內(nèi)容系統(tǒng)電路如圖3所示。電路元件參數(shù)：直流電壓源，電阻，電感，電容。電路元件初始值：電感電流，電容電壓。系統(tǒng)輸出量為電容電壓。試?yán)脷W拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法構(gòu)建系統(tǒng)仿真模型，并求出離散系統(tǒng)的輸出量響應(yīng)曲線。連續(xù)系統(tǒng)輸出響應(yīng)的解析解為： （2-1）其中， ， 。 三、實驗要求1）利用歐拉法、梯形法、二階顯式Adams法及顯式四階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)的標(biāo)題。四、實驗原理在連續(xù)系統(tǒng)的數(shù)字仿真算法中,較常用的有歐拉法、梯形法、二階顯式Adams法及顯式四階Runge-Kutta法等。歐拉法、梯形法和二階顯式Adams法是利用離散相似原理構(gòu)造的仿真算法，而顯式四階Runge-Kutta法是利用Taylor級數(shù)匹配原理構(gòu)造的仿真算法。對于線性系統(tǒng)，其狀態(tài)方程表達式為： (4-1)式（4-1）中，是系統(tǒng)的n維狀態(tài)向量, 是系統(tǒng)的m維輸入向量，是系統(tǒng)的r維輸出向量。A為階參數(shù)矩陣，又稱動態(tài)矩陣，B為階輸入矩陣，C為階輸出矩陣，D為階交聯(lián)矩陣。利用前向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為： （4-2）式中，為積分步長，為單位矩陣。利用后向歐拉法構(gòu)建線性系統(tǒng)的仿真模型為： （4-3）利用梯形法構(gòu)建線性系統(tǒng)的仿真模型為： （4-4）利用二階顯式Adams法構(gòu)建線性系統(tǒng)的仿真模型為： （4-5）式中： （4-6）二階顯式Adams法為多步計算方法，利用多步計算方法對系統(tǒng)進行仿真時，需要與之具有相同計算精度的單步計算方法輔助計算。二階顯式Adams法的計算精度為二階，可以采用梯形法或改進的Euler法等輔助計算。利用改進的Euler法構(gòu)建線性系統(tǒng)的仿真模型為： （4-7）利用顯式四階Runge-Kutta法構(gòu)建線性系統(tǒng)的仿真模型為： （4-8）五、實驗結(jié)果取不同的積分步長h，仿真結(jié)果如下：(1) 仿真時間t=0.01s，積分步長h=10-6s（2） 仿真時間t=0.01s，積分步長h=2*10-6s（3）仿真時間t=0.01s，積分步長h=10-5s（4）仿真時間t=0.01s，積分步長h=5*10-5s（5）仿真時間t=0.01s，積分步長h=10-4s（6）仿真時間t=0.01s，積分步長h=2*10-4s（7）仿真時間t=0.01s，積分步長h=4*10-4s六、實驗結(jié)果分析在h=5e-5,h=5e-6,h=5e-7時得到的圖像，可以看出，在h=5e-5時，前向歐拉法和后向歐拉法圖像明顯誰真，而梯形法和二階顯示Adams法圖像有輕微失真，步距仍然較大。在步距為時，前向歐拉法和后向歐拉法圖像有部分失真，前向歐拉法失真較嚴(yán)重。而梯形法和二階顯示Adams法圖像與連續(xù)型函數(shù)曲線相似度極高，仿真效果非常好。當(dāng)步距為時，所有方法仿真效果都非常好，與連續(xù)型函數(shù)曲線相似度都極高?？梢姡骄鄷r，所有仿真方法都可以應(yīng)用。從仿真模型實現(xiàn)的難易性、模型的穩(wěn)定性、模型的精度及離散時間間隔等方面，對比分析上述方法構(gòu)造的離散系統(tǒng)模型的優(yōu)缺點。 難易性：通過單個模型的分別仿真，可以得出顯式四階Runge-Kutta法建模最為復(fù)雜，仿真時間也較長，對步距要求較低。其次復(fù)雜的是二階顯式Adams法和梯形法，仿真時間稍短，梯形法取梯形面積，誤差也較小。前向歐拉法和后向歐拉法模型的復(fù)雜程度差不多，仿真時間也差不多。模型的穩(wěn)定性：當(dāng)步距h=5.0e-5時，前向歐拉法和后向歐拉法明顯失真，隨著步距的減小，二階顯式Adams法，梯形法和顯式四階Runge-Kutta法的波形變化不大，而前向歐拉法和后向歐拉法的波形得到明顯改善。所以顯式四階Runge-Kutta法，二階顯式Adams法和梯形法的穩(wěn)定性較好，前向歐拉法和后向歐拉法的穩(wěn)定性較差。模型的精度和離散時間間隔：步距為h=5.0e-6時，顯式四階Runge-Kutta法精度最高，其次是二階顯式Adams法和梯形法。步距為h=5.0e-7時，前向歐拉法和后向歐拉法仿真精度才達到要求。所以，顯式四階Runge-Kutta法，二階顯式Adams法和梯形法模型的精度較高，離散時間間隔要求低，其中，顯式四階Runge-Kutta法模型的精度最高，其次是二階顯式Adams法，由于是二次函數(shù)較復(fù)雜，函數(shù)曲線與真實曲線較為接近；再次精確的是梯形法，取梯形面積，誤差也較??；前向歐拉法和后向歐拉法模型的精度較低，由于取的是矩形面積，離散時間間隔要求高。 七、實驗結(jié)論二階顯示Adams法精度最高，其次是梯形法。在h = 5.0e-6已然如此，當(dāng)步長在小的時候反而使運行時間延長，效果不一定好。這就說明方法的選取與精度和復(fù)雜性，以及可行性等密切相關(guān)。八、附錄m文件源程序：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;% 前向歐拉法 % for i=1:1:n x1(1:n,1) = 0; end for k=1:m x1(1:n,k+1) = x1(1:n,k) + (A* x1(1:n,k)+B)*h; end for k=1:1:m y1(k) = D*x1(1:n,k); end % 后向歐拉法 % 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 % 梯形法 % 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 % 二階顯示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)*h/12; end for k=1:1:m y4(k) = D*x4(1:n,k); end % 四階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*( 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 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,1),plot(t,y,g,t,y1,r) legend(y解析解,y1前向歐拉) title(前向歐拉法) subplot(2,3,2),plot(t,y,g,t,y2,r) legend(y解析解,y2后向歐拉) title(后向歐拉法) subplot(2,3,