1、仿真教學PPT課件 3-1 基于數值積分法的連續系統仿真基于數值積分法的連續系統仿真 3-2 基于離散相似法的連續系統仿真基于離散相似法的連續系統仿真 3-3 系系統非線性環節的仿真統非線性環節的仿真 仿真教學PPT課件3.1 基于數值積分法的基于數值積分法的連續系統仿真連續系統仿真 3.1.1 數值積分基本原理數值積分基本原理 3.1.2 數值積分方法的選擇數值積分方法的選擇 3.1.3 基于數值積分法的連續系統仿真基于數值積分法的連續系統仿真 3.1.4 數值積分法的數值積分法的MATLAB函數函數仿真教學PPT課件 連續系統的數學模型一般都能以微分方程的連續系統的數學模型一般都能以微分方

2、程的形式給出，所以連續系統仿真算法問題通常可歸形式給出，所以連續系統仿真算法問題通常可歸結為如何用計算機來求解微分方程的問題，也就結為如何用計算機來求解微分方程的問題，也就是對一階微分方程如何進行求解。是對一階微分方程如何進行求解。是是解決該問題的重要方法之一。解決該問題的重要方法之一。 例如，假設有一系統，它的數學模型可用例如，假設有一系統，它的數學模型可用（1 1）式所示微分方程來描述）式所示微分方程來描述3.1.1 數值積分基本原理數值積分基本原理3.1 基于數值積分的基于數值積分的連續系統仿真連續系統仿真仿真教學PPT課件 數數 值值積積 分分方方 法法是是解解 決決在在 已已知知 初

3、初值值 的的情情 況況下下 ，對對),( tyf進進行行 近近 似似積積分分 ，對對 )(ty進進 行行數數值值求求解解的的方方法法。在在數數學學上上稱稱為為微微分分方方程程初初值值問問題題的的數數值值方方法法。 仿真教學PPT課件仿真教學PPT課件仿真教學PPT課件仿真教學PPT課件仿真教學PPT課件三、龍格三、龍格庫塔法（庫塔法（RungeKutta法）法） 龍格龍格庫塔法的基本思想是：庫塔法的基本思想是：用幾個點上的函用幾個點上的函數值的線性組合來代替函數的各階導數，然后再按數值的線性組合來代替函數的各階導數，然后再按泰勒級數展開確定其中的系數。泰勒級數展開確定其中的系數。 下面以二階下

4、面以二階RK法為例介紹其基本原理。法為例介紹其基本原理。)()(,00tyttytyfdtdy( 5 )仿真教學PPT課件將將K2用二元函數泰勒級數展開式展開，并只取前三用二元函數泰勒級數展開式展開，并只取前三項，則有：項，則有：mmmmyfhKbtfhbtyfK121,2)(首先假設式（首先假設式（1 1）的解具有下面的形式：）的解具有下面的形式：),(),()(1122122111hbthKbyfKtyfKKaKahyymmmmmm式中，式中，a a1 1、a a2 2、b b1 1、b b2 2為待定系數。為待定系數。仿真教學PPT課件將將 K1、K2 代入式代入式（5），），得：得：m

5、mmmyfhKbtfhbtyfK121,2)(),(),()(1122122111hbthKbyfKtyfKKaKahyymmmmmmmmmmmmmmmmmmmmmmmyfhtyfbatfhbatyhfaayyfhtyfbtfhbtyfhatyhfayy2222122121211),(),()(),(),(),( 6 )仿真教學PPT課件比較式比較式（5）與與（6）可以得到：可以得到：另一方面，將另一方面，將ym+1在在ym附近進行泰勒展開，并只取附近進行泰勒展開，并只取前三項，則有：前三項，則有：mmmmyhyhyy 21212,)(21)(hyftyftftyhfymmmmmmm21211

6、221221babaaa仿真教學PPT課件 上述上述3個方程中有個方程中有4個未知數，因而解不是唯一個未知數，因而解不是唯一的。若限定的。若限定 ，則可得其中一組解：，則可得其中一組解：21aa 1,212121bbaa將它們代入式（將它們代入式（5），可得一組計算公式：），可得一組計算公式：),(),()(2121211hthKyfKtyfKKKhyymmmmmm仿真教學PPT課件截斷誤差為截斷誤差為 O(hO(h3 3) )截斷誤差為截斷誤差為 O(hO(h5 5) )仿真教學PPT課件（6 6） 梯形公式可看作二階龍格梯形公式可看作二階龍格庫塔公式，截斷誤差正比庫塔公式，截斷誤差正比h3

7、（2 2）步長）步長 h可變。可變。RK的特點：的特點：仿真教學PPT課件四、亞當姆斯（四、亞當姆斯（AdamsAdams）法）法 在計算在計算ym+1時，只利用前一步的時，只利用前一步的ym的值，經過若干的值，經過若干步的計算以后，可以求出一系列的值步的計算以后，可以求出一系列的值y1，y2，.，yn。如果充。如果充分利用前面多步的值來計算分利用前面多步的值來計算ym+1，則可以達到既提高計算速度，則可以達到既提高計算速度又能獲得較高精度的目的，這就是又能獲得較高精度的目的，這就是的基本思想。在多步的基本思想。在多步法中，應用較廣的是法中，應用較廣的是Adams法。法。 Adams法有法有和

8、和兩種。兩種。 Adams顯式一般形式：顯式一般形式： （5） 式中，式中， 為顯式為顯式Adams公式系數，部分數據如下表公式系數，部分數據如下表rmrmmmmfffhyy1101i仿真教學PPT課件仿真教學PPT課件如二次（如二次（Adams）亞當姆斯公式：）亞當姆斯公式：三次（三次（Adams）亞當姆斯公式：）亞當姆斯公式：即是即是多步法計算公式多步法計算公式。 多步法與單步法相比，欲達相同精度，計算工多步法與單步法相比，欲達相同精度，計算工作量較少，在相同條件下多步法比單步法要快。作量較少，在相同條件下多步法比單步法要快。2115162312mmmmmfffhyy32119375955

9、24mmmmmmffffhyy仿真教學PPT課件計算穩定性計算穩定性 仿真計算時，是否仍然穩定呢？先看下面的例子：仿真計算時，是否仍然穩定呢？先看下面的例子： 從穩定性理論，我們知道如何去從系統的微分從穩定性理論，我們知道如何去從系統的微分方程或傳遞函數去判斷該系統的穩定性。那么，對方程或傳遞函數去判斷該系統的穩定性。那么，對于一個穩定的連續系統，當用某數值積分方法進行于一個穩定的連續系統，當用某數值積分方法進行五、數字仿真中的幾個問題五、數字仿真中的幾個問題仿真教學PPT課件仿真教學PPT課件（4）仿真教學PPT課件仿真教學PPT課件仿真教學PPT課件小結：小結： 一般穩定性與步長關系密切（

10、除恒穩公式之外），若一般穩定性與步長關系密切（除恒穩公式之外），若用兩種顯著不同的步長所得到的數值解有明顯差別，則可用兩種顯著不同的步長所得到的數值解有明顯差別，則可能是這種數值方法不穩定；反之如果基本相同，則一般視能是這種數值方法不穩定；反之如果基本相同，則一般視為是穩定的為是穩定的 。仿真教學PPT課件仿真教學PPT課件仿真教學PPT課件 3.1.2 數值積分方法的選擇數值積分方法的選擇一、積分方法的選擇一、積分方法的選擇在步長相同的條件下，積分方法的在步長相同的條件下，積分方法的階數越高，精度越高；另外，多步法的精度比單步法階數越高，精度越高；另外，多步法的精度比單步法高，隱式算法的精度

11、高于顯式算法。因此，當需要高高，隱式算法的精度高于顯式算法。因此，當需要高精度時，可采用高階的多步隱式算法和較小的步長。精度時，可采用高階的多步隱式算法和較小的步長。若精度要求不高，一般可選擇低階算法。若精度要求不高，一般可選擇低階算法。為加快計算速度，可在精度要求不高為加快計算速度，可在精度要求不高時，盡量選擇低階的計算工作量少的方法。時，盡量選擇低階的計算工作量少的方法。仿真教學PPT課件數值解的穩定性必須保證，否則計算結果將數值解的穩定性必須保證，否則計算結果將失去真實意義。不同的方法有不同的穩定性，要失去真實意義。不同的方法有不同的穩定性，要通過合適選擇步長來保證穩定性。通過合適選擇步

12、長來保證穩定性。單步法有自啟動能力，多步法沒有自啟動能單步法有自啟動能力，多步法沒有自啟動能力，必須借助于單步法啟動運算之后，才能開始力，必須借助于單步法啟動運算之后，才能開始工作。一般簡單的仿真程序多用單步法。工作。一般簡單的仿真程序多用單步法。仿真教學PPT課件單步法在整個計算中，步長可在一定范圍內單步法在整個計算中，步長可在一定范圍內變化；而多步法則對步長的變化有嚴格的要求。變化；而多步法則對步長的變化有嚴格的要求。若要求仿真時步長可變，最好用單步法。若要求仿真時步長可變，最好用單步法。 綜上所述，積分方法的選擇與多種因素有關，綜上所述，積分方法的選擇與多種因素有關，各因素之間又相互影響

13、。究竟選哪一種方法，要各因素之間又相互影響。究竟選哪一種方法，要由具體系統及具體要求而定。由具體系統及具體要求而定。仿真教學PPT課件 二、積分步長的選擇二、積分步長的選擇 步長的選擇很重要。步長過大會增大截斷誤差，步長的選擇很重要。步長過大會增大截斷誤差，甚至出現數值不穩定現象，甚至出現數值不穩定現象， 過小了又因增加了步過小了又因增加了步數，而使舍入誤差增大。所以存在一個最佳步長使數，而使舍入誤差增大。所以存在一個最佳步長使得總誤差最小（如圖所示）。得總誤差最小（如圖所示）。 在用經驗方法選取步長時，一種方法是根據系在用經驗方法選取步長時，一種方法是根據系統方程中統方程中最小時間常數最小時

14、間常數Tmin來決定步長，一般取：來決定步長，一般取： h=（0.20.05）Tmin仿真教學PPT課件數值積分計算時，積分步長有數值積分計算時，積分步長有固定步長固定步長和和變步長變步長兩種工作方式。兩種工作方式。 就是在整個仿真計算過程中，積分步長就是在整個仿真計算過程中，積分步長 h 始終保持不變。始終保持不變。 就是在仿真積分計算的每一步，根據計算誤差的大小改變就是在仿真積分計算的每一步，根據計算誤差的大小改變步長步長h。目的目的：在保持一定的計算精度的前提下，盡可能地選取較大：在保持一定的計算精度的前提下，盡可能地選取較大的步長。的步長。方法方法：首先估計計算誤差：首先估計計算誤差

15、；判斷誤差是否在允許的誤差；判斷誤差是否在允許的誤差 0范范圍內；若在允許誤差范圍內，則該步計算有效，否則計算無圍內；若在允許誤差范圍內，則該步計算有效，否則計算無效，改變步長，重新計算。效，改變步長，重新計算。 因此，在變步長的積分計算中，必須解決因此，在變步長的積分計算中，必須解決誤差估計方法誤差估計方法和和步長調整策略步長調整策略兩個問題。兩個問題。仿真教學PPT課件三、直接采用三、直接采用MATLAB語言根據算法原理編程語言根據算法原理編程 31030yydtdycleart0=0;y0=1/3;h=0.1;N=1.5/h;t(1)=0;y(1)=1/3;for i=2:N+1 k1=

16、 -30*y0; y1=y0+h*k1; y0=y1; y(i)=y1; t(i)=t(i-1)+h;endy1cleart0=0;y0=1/3;h=0.1;N=1.5/h;t(1)=0;y(1)=1/3;for i=2:N+1 k1=-30*y0; k2=-30*(y0+h*k1/2); k3=-30*(y0+h*k2/2); k4=-30*(y0+h*k3); y1=y0+h*(k1+2*k2+2*k3+k4)/6; y0=y1; y(i)=y1; t(i)=t(i-1)+h;endy1stableEuler.m stableRK.m仿真教學PPT課件START讀入讀入 、 、 、 0t0

17、ynh m=0),()2,2()2,2(),(3423121hthkyfkhtkhyfkhtkhyfktyfkmmmmmmmm )22(643211kkkkhyymmmn-1 打印打印1my m=m+1RK程序框圖程序框圖仿真教學PPT課件習習 題題3-1 用用Euler法求初值問題法求初值問題1,01)0(2tuutudtdu的數值解。取步長的數值解。取步長 h=0.1。Euler法：法： u（1） =1.7848RK法：法： u（1） =1.7321 仿真教學PPT課件RK法：法：86466447. 0)2(yEuler法法：86806219. 0)2(y精確解：精確解：86466472.

18、 0)2(y習習 題題3-2 用用Euler法（法（h=0.025）及經典）及經典Runge-Kutta法法（h=0.1）計算初值問題：）計算初值問題：0)0(1yyy在在 t=2 時刻的值，并和精確解比較。時刻的值，并和精確解比較。仿真教學PPT課件仿真教學PPT課件對一個對一個n維向量維向量X，每前進一個步距，每前進一個步距h，至少要求，至少要求n4個個ijk的值。的值。仿真教學PPT課件3.1.4 數值積分法的數值積分法的MATLAB函數函數 MATLAB工具箱提供了各種數值積分方法的常用函數，這些函數的基工具箱提供了各種數值積分方法的常用函數，這些函數的基本功能是用數值計算方法求解常系

19、數微分方程（本功能是用數值計算方法求解常系數微分方程（ODE）或微分方程組。）或微分方程組。ODE45一種顯式一種顯式RK（4-5）公式，單步法，變步長，適用于求解非）公式，單步法，變步長，適用于求解非剛性微分方程。對于大多數問題能獲得滿意的解。剛性微分方程。對于大多數問題能獲得滿意的解。ODE23一種顯式一種顯式RK（2-3）公式，單步法，變步長，適用于求解非）公式，單步法，變步長，適用于求解非剛性微分方程。在允許計算誤差較大或方程具有輕微剛性時，剛性微分方程。在允許計算誤差較大或方程具有輕微剛性時，效果比效果比ODE45好。好。ODE113一種變階一種變階Adams-Bashforth-M

20、oulon算法，多步法，適用于求算法，多步法，適用于求解非剛性微分方程。在允許計算誤差較小時，比解非剛性微分方程。在允許計算誤差較小時，比ODE45更有更有效。效。ODE15s多步法，用于求解多步法，用于求解剛性微分方程剛性微分方程。當。當ODE45、ODE113無法無法解決問題時，可嘗試用解決問題時，可嘗試用ODE15s求解。求解。ODE23s單步法，用于求解單步法，用于求解剛性微分方程剛性微分方程。若在使用。若在使用ODE15s無效情無效情況下，可嘗試使用況下，可嘗試使用ODE23s來求解某些剛性方程。在允許計來求解某些剛性方程。在允許計算誤差較大時，比算誤差較大時，比ODE15s更有效。

21、更有效。仿真教學PPT課件函數調用格式：函數調用格式：說明說明：（1）solver為數值積分函數名為數值積分函數名，即，即ODE45、ODE23、ODE113、ODE15s、ODE23s。（2）F是描述常微分方程的函數是描述常微分方程的函數m文件名文件名。（3）tspan為一個行向量，描述運算的起止時間為一個行向量，描述運算的起止時間。例如。例如0 20。（4）y0為初始值為初始值。)0, (,ytspanFsolveryt仿真教學PPT課件解決問題的一般步驟：解決問題的一般步驟： 首先把描述系統的高階微分方程轉變為一階微首先把描述系統的高階微分方程轉變為一階微分方程組形式。分方程組形式。 建

22、立一個建立一個函數函數m m文件文件用于描述該一階微分方程用于描述該一階微分方程 組。組。1.1. 再建立一個再建立一個腳本腳本m m文件文件，調用某個，調用某個solversolver函數求函數求解微分方程。解微分方程。仿真教學PPT課件【例例3.4-1】 已知二階微分方程已知二階微分方程0)1 (2yyyy 初始條件：初始條件：1)0(,0)0(yy求時間區間求時間區間t=0 20微分方程的解微分方程的解。解：分三個步驟求解解：分三個步驟求解（1）將微分方程表示為一階微分方程組）將微分方程表示為一階微分方程組1221221)1(yyyyyy（2）建立描述微分方程組的）建立描述微分方程組的函

23、數函數m文件文件vdp.m function dy=vdp(t,y) dy=y(2);(1-y(1)*y(1)*y(2)-y(1);（3）建立）建立腳本腳本m文件文件example23.m，調用解題器指令求解調用解題器指令求解yt,y=ode45(vdp,0 20,0,1);plot(t,y(:,1),r-,t,y(:,2),b:);xlabel(t);ylabel(y);legend(y1,y2);仿真教學PPT課件【例例3.4-2】 已知一個剛體在無外力作用下運動方程及初始條件為已知一個剛體在無外力作用下運動方程及初始條件為2133123212yyyyyyyyy初始條件：初始條件：1) 0(, 1) 0(,0) 0(321yyy求時間區間求時間區間t=0 12微分方程的解。微分方程的解。解：分兩個步驟求解解：分兩個步驟求解（1）建立描述微分方程組的函數）建立描述微