連續(xù)系統(tǒng)數值積分仿真方法學

1、第三章 連續(xù)系統(tǒng)數值積分仿真方法學 3.1 連續(xù)系統(tǒng)數值積分法基本原理 3.2 Runge-Kutta積分法 3.4 數值積分法穩(wěn)定性分析 3.5 數值積分法的選擇與計算步長的確定 3.6 通用仿真程序8/29/20221概述連續(xù)系統(tǒng)仿真：借助于計算機，求解數學模型關鍵問題：如何選擇一種算法，將系統(tǒng)模型轉化為能在計算機上運行的離散模型兩類仿真算法：數值積分法：歐拉法、龍格-庫塔法離散相似法：離散狀態(tài)法、屠斯丁法8/29/20222 連續(xù)系統(tǒng)的數學模型一、微分方程：常見：單變量、線性定常微分方程。 y+5y+6=3x它是根據數學力學原理推導而出二、傳遞函數G（S）=Y（S）/X（S）三、狀態(tài)方程

2、中間狀態(tài)輸入之關系輸出方程：中間狀態(tài)輸出關系8/29/20223連續(xù)系統(tǒng)仿真步驟1 建立數學模型：微分方程、狀態(tài)方程、傳遞函數2 建立仿真模型：選擇仿真算法，將數學模型轉化為仿真模型3 編寫運行程序：選擇某種語言，寫出計算機仿真程序，進行運行，獲得仿真結果 8/29/20224 數值積分法系統(tǒng)仿真數值積分法：面向方程的系統(tǒng)仿真、面向結構圖的系統(tǒng)仿真。數值積分法：8/29/20225數值積分法的任務：尋求真解在一系列離散點上數值解數值積分法要解決的問題：如何對函數F（t，y）進行數值積法基本方法：單步法、多步法、預估校正法（精度、速度、穩(wěn)定性不一樣）基本思路：以歐拉法為例進行說明。8/29/20

3、226在若干離散點處，求出若干的y值來代替連續(xù)變量y（t），數值積分法就是一種把連續(xù)變量離散成離散變量的一種方法。一、 歐拉法8/29/202278/29/20228舉例由于是近似的，肯定有誤差，可減小步長h=tK+1-tK，但步長過小，計算量太大，通過此法提高精度是有限的。8/29/20229二、 改進歐拉法（梯形法）用梯形面積代替小區(qū)間曲線積分，可以比歐拉法提高精度。8/29/202210一次迭代求近似解先用歐拉法進行預估，再代入梯形公式它是一種多步法，不能自啟動，必須用其它方法求以前的多步解。以上只適于一階微分方程組，高階要降階處理。8/29/202211三、 數值積分法的幾個基本概念1

4、 算法自啟動只要知道初值，能遞推一系列離散解，無需其它算法的輔助（如歐拉）2 單步法與多步法采用的遞推公式，都是步進的，如果由Yn能求得Yn+1，不需其它信息（單步）如果，為了求Yn+1，需要Yn、Yn-1等信息（多步）8/29/2022123 顯式與隱式顯式：Yn+1遞推公式中，用不到Yn+1的數值隱式：Yn+1的遞推公式中，用到了Yn+1的數值（需要預估）4 截斷誤差（泰勒展開）歐拉法是一階精度的，改進的歐拉法是二階精度的 8/29/2022135 舍入誤差：由于計算機字長引起的誤差6 初始誤差：給定的初始值與真值之間的差異7 數值解的穩(wěn)定性：由于誤差的影響，計算步長選擇不合適，有可能使數

5、字仿真結果出現不穩(wěn)定的現象給定步長h，如果計算結果對初始誤差或計算誤差不敏感，算法是穩(wěn)定的。不穩(wěn)定的算法，誤差會惡性膨脹，計算結果發(fā)散，仿真失敗。8/29/2022143.2龍格-庫塔法（RUNGE-KUTTA）間接利用泰勒展開式，用幾個點上的斜率線性組合代替Y（t）的各階導數然后用TAYLOR級數展開式確定線性組合中各加權系數既提高了精度，又省去了高階導數計算8/29/202215一、龍格-庫塔數值積分公式將y（t）在t0，y0處用臺勞級數展開，保留h2及以前的項。8/29/2022168/29/2022178/29/202218另一套遞推公式8/29/202219保留到h4，此時龍格庫塔公

6、式8/29/2022208/29/202221龍格-庫塔法的特點A 單步法可自啟動，即知道初值后，直接從微分方程初值開始計算下去。B 不同步驟的步長可以不一樣，但每一步需同一步長來計算系數。C 計算YK+1分兩部分：求YK；步長乘各系數的加權平均值8/29/202222D 精度取決于h的大小及解決方法當精度相同時，四階的計算量是二階的兩倍，但步長可以是二階的10倍。當然，也不是階數越高越好，高于四階，計算量猛增，一般也不用。E 歐拉公式實質就是一階龍格-庫塔公式總之，四階，精度高，計算量適中，編程容易，步長易改變，穩(wěn)定性較好。8/29/202223二、四階龍格-庫塔法的向量形式用于微分方程組（

7、多維）向量形式，n階系統(tǒng)微分方程：8/29/2022243.4 數值積分法穩(wěn)定性分析情況：系統(tǒng)本來是穩(wěn)定的，可仿真結果卻是發(fā)散的原因：積分步長選得不合適造成的穩(wěn)定性含義：在擾動影響下，其計算過程中的累積誤差不會隨計算步長的增加而無限增長。一個數值法是否穩(wěn)定取決于該差分方程的特征根是否滿足穩(wěn)定性要求8/29/202225二、穩(wěn)定性分析1 歐拉法的前差公式：積分步長H必須小于系統(tǒng)時間常數的兩倍2 后差公式：恒穩(wěn)定3 梯形公式：恒穩(wěn)定4 中心差公式：不存在穩(wěn)定區(qū)域8/29/202226三、 數值積分法的選擇如何選擇，尚未有具體方法。幾點原則：A 精度問題精度：截斷誤差、舍入誤差、累積誤差一般，步長小

8、，截斷誤差??；計算機字長長，舍入誤差小；累積誤差與積分時間長短有關（步長小，累積誤差大）在一定積分方法條件下，綜合考慮總誤差，使之最小化。8/29/202227B 計算速度問題取決于步長，以及每一步所需時間，每一步的時間取決積分方法，它主要取決于計算導數的次數。C 數值解穩(wěn)定性問題數值積分法實際是將微分方程化為差分方程來求解。一個穩(wěn)定的微分方程，相應的差分方程不一定穩(wěn)定。不同的積分方法，具有不同的穩(wěn)定域8/29/202228D 自啟動問題自啟動：直接從微分方程的初值開始。多步法：需要初始值加上以前時刻的值，無法自啟動?？傊?，數值積分法的選擇同許多因素有關，究竟選擇哪一種要具體分析。一般采用四階

9、龍格庫塔法當對仿真問題了解不多時，也可先采用變步長。8/29/2022293.6 通用仿真程序一、仿真程序的主要功能和結構1 它能將輸入的系統(tǒng)模型及參數2 一般應能提供多種仿真算法供選擇3 將仿真結果進行存儲、打印、繪圖4 具有重復運行的控制功能，以便研究系統(tǒng)參數變化對系統(tǒng)性能的影響結構：主模塊初始化、運行程序（數值積分子程序、顯示子程序）、輸出程序（數據顯示、繪圖程序）8/29/202230二、面向方程的系統(tǒng)仿真1 狀態(tài)方程的仿真程序2 傳遞函數3 微分方程8/29/2022311 狀態(tài)方程的仿真程序設系統(tǒng)的狀態(tài)方程和輸出方程：8/29/202232X(t)n維列向量;An*n系數矩陣Bn*

10、1控制矩陣;C1*n輸出矩陣U(t)輸入變量;y(t)輸出變量采用四階龍格庫塔法進行仿真計算8/29/202233龍格庫塔公式8/29/202234每個微分方程有四個系數一共有4n個系數若已知初值x0和輸入函數u（I），即可求各時刻的xK和yK。一個完整的仿真程序：初始化程序、執(zhí)行運算的計算程序8/29/202235初始化程序在仿真開始時，對模型初始化，改變仿真時鐘；設置狀態(tài)變量初值和參數值，步長，打印間隔，仿真終止時間；輸出程序按指定格式輸出仿真結果如圖形、表格、曲線等。一般對一個復雜的仿真程序要分成若干子程序，分別進行調試和修改。8/29/202236子程序清單C subroutine o

11、f subroutine RK(N,K,A,B,C,X,D,Z,P) real K(n,0:4),A(N,N),B(N),C(N),X(N) real H,T,T1 integer Q 打印間隔步數 輸入H,T,T1,Q,A(I,J),B(I),X(I)8/29/202237 do 5 I=1,n do 5 j=1,n k(I,j)=.05 continue L(1)=.0 L(2)=H/2 L(3)=h/2 L(4)=h goto 258/29/2022386 do 20 m1=1,Q do 200 j=1,4 T2=T+L(J) do 150 I=1,N P(I,j)=L(j)*K(I,J

12、-1) Z(I)=P(I,J)*X(I)150 continue8/29/202239 do 200 I=1,N D(I,J)=.0 do 160 m=1,N D(I,J)=D(I,J)+A(I,M)*T(M)160continue K(I,J)=D(I,J)+B(I)*U(T2)200continue8/29/202240 do 10 I=1,N X(I)=X(I)+H*(K(I,1)+2*K(I,2) +2*K(I,3)+K(I,4)/610 continue T=T+H IF(T.GT.T1)goto 2520 continue8/29/20224125 y=.0 do 30 I=1,N

13、 y=y+C(I)*X(I)30 continue write (*,40)T,Y40 format ( ) if(T.LT.T1)goto 6 return end8/29/202242例：連續(xù)系統(tǒng)微分方程步長為0.001秒，打印間隔為0.01秒終止時間為1秒。試用四階龍格庫塔法編制其仿真程序8/29/202243狀態(tài)方程和輸出方程8/29/202244主程序C simulation dess in fortran 77C main frogram real K(2,0:4),A(2,2),B(2),C(2),X(2), D(2,4),z(2),P(2,4) N=2 call RK(N,K,

14、A,B,C,X,D,Z,P) end8/29/2022452 傳遞函數的仿真程序先將傳遞函數轉換成狀態(tài)方程，再用前面的方法仿真一般采用相變量法獲得狀態(tài)方程和輸出方程：X=AX+BU Y=CX帶有負反饋系數v的閉環(huán)系統(tǒng)：u=r-vy3 微分方程的仿真程序容量推廣到含有非線性系統(tǒng)的仿真中組成:主程序、微分方程子程序、運行子程序、積分子程序、輸出子程序8/29/202246三、面向混合模型的仿真程序分析GSPGeneral Simulation Program北京理工大學仿真實驗室開發(fā)GSP 在 windows下運行，用VB語言編寫，菜單驅動 界面標準 使用方便下面簡單地看一下數值積法算法程序該程序

15、是仿真軟件中的核心程序8/29/202247RK4程序Sub rk4()Dim hi(4) as single, z(30)as sigle, T1 as single Hi(1)=0 : Hi(2)=0.5*h:Hi(3)=0.5*h : Hi(4)=hFor i =1 to nxZ(i)=xx(i)Next i8/29/202248T1=tFor j = 1 to 4Istep = j t=T1+Hi( j )For i = 1 to nxXx( i )=z( i )+Hi( j )*dxx( i )Next iDeriv右函數子程序8/29/202249For i = 1 to nxXk( j , i )=dxx( i )Next iNext j For i = 1 to nxXx( i )=z( i )+h*xk(1,i)+2*xk(2,i) +2*xk(3,i)+xk(4,i)/6Next iEnd sub8/29/202250輔助函數及公用子程序模型文件（如右函數子程序DERIV）作業(yè)3.13.23.33.48/29/202251RK4 PROGRAMFunction f1 (x1,x2,x3,x4,t as single)F1=End function(以上為公用部分)Dim h as singleDim t as singleDim I as integerH=