計算機仿真技術基礎第3章-連續(xù)系統(tǒng)數(shù)值積分仿真方法學 課件

1、計算機仿真技術基礎 第三章連續(xù)系統(tǒng)數(shù)值積分仿真方法學第二章 連續(xù)系統(tǒng)數(shù)值積分仿真方法學第一節(jié) 數(shù)值積分法的基本原理第二節(jié) 數(shù)值積分法的單步算法第三節(jié) 數(shù)值積分法的多步算法如何把已建立起來的數(shù)學模型轉(zhuǎn)換成仿真運算模型（二次建模），以便為分析解決實際問題服務那是系統(tǒng)仿真學科的一個重要研究內(nèi)容。對于復雜的數(shù)學模型來說，求其解析解是很煩瑣和困難的，大多數(shù)情況下不求出解析解，或者根本不存在解析解，因此借助于數(shù)值解法對連續(xù)系統(tǒng)進行仿真研究。 用計算機不可能求出系統(tǒng)的解析解（連續(xù)）只能求出連續(xù)響應曲線上的有限個點，即數(shù)值解連續(xù)系統(tǒng)數(shù)值積分法：就是利用數(shù)值積分方法對常微分方程建立離散化形式的數(shù)學模型（差分方程

2、）并求出數(shù)值解。 最常用的數(shù)值解法有： 歐拉法、梯形法、Adams、RungeKutta法。 描述各類系統(tǒng)最基本的模型用微分方程或狀態(tài)空間表達式，二次建模就是要求出適合用數(shù)字計算機求解的模型，就需要把微分運算轉(zhuǎn)化成算術運算在用計算機求解。第一節(jié) 數(shù)值積分法的基本原理首先把需仿真研究的系統(tǒng)表示成一階微分方程組或狀態(tài)方程的形式。 以一階連續(xù)系統(tǒng)為例，微分方程及初值如下： 設方程在處的連續(xù)解為：希望找到一個近似公式來表示方程的近似解：為精確值的近似值為準確積分值的近似值所謂微分方程初值問題的數(shù)值解法就是尋求真解在一系列離散點上的近似解（數(shù)值解）相鄰兩個時間離散點的間隔稱為計算步距或步長，通常情況為定

3、值，也有變步長?？梢?，微分方程初值問題數(shù)值解法的主要問題歸結(jié)為對為此，要先把連續(xù)的微分方程用數(shù)值積分法轉(zhuǎn)化為離散的差分方程的初值問題，然后根據(jù)初始條件X0逐步遞推計算出后續(xù)時刻的數(shù)值解：如何求出定積分的近似解數(shù)值解法的共同特點是步進式的遞推算法，主要有：單步法、多步法和預估校正法，并有顯式和隱式之分。一、歐拉法（Euler Method)Euler法是最簡單的一種數(shù)值積分法的單步運算，雖然計算精度較差，但幾何意義明顯，便于理解，能說明構(gòu)造數(shù)值積分算法的基本思想。下面采用三種方法推導出Euler法的數(shù)值近似公式，以便對數(shù)值積分器的基本思想能透徹了解。第二節(jié) 數(shù)值積分法的單步算法 、Taylor級

4、數(shù)展開x(t)為解析解，將x(t)展開成Taylor級數(shù) 以一階連續(xù)系統(tǒng)為例，微分方程及初值如右：只取一次項，其余忽略寫成差分方程為這就是解微分方程初值問題的歐拉算法。、矩形近似把積分區(qū)間 h 取得足夠小，將 在 近似為常數(shù)用左矩形面積近似該區(qū)間的曲線面積對方程在上求積分也能得到：左矩形（也稱為前向歐拉法）近似及誤差將 在 近似為常數(shù)用右矩形面積近似該區(qū)間的曲線面積得到：這是右矩形歐拉公式，是一個隱式算法對積分右矩形（也稱為后向歐拉法）近似及誤差、切線近似（）取切線上處的值來近似在 的一個小鄰域內(nèi)，曲線x(t)可用 處的切線來表示， x(t) 在 處的斜率為：也能得到：前向歐拉法過點 以 fn

5、 為斜率的切線方程為：、切線近似（）過點 以 fn+1 為斜率的切線方程為：取切線上處的值來近似在 曲線 x(t) 可用 處的切線來表示， x(t) 在 處的斜率為：也能得到：后向歐拉法歐拉法（切線推導）的幾何意義歐拉法實際計算時的幾何意義例：設系統(tǒng)方程試用Euler法求其數(shù)值解，取步長 h=0.1,解:前向Euler法遞推式：有初始條件：可進行遞推：后向Euler法遞推式：隱式算法，需先解此非線性方程：由此公式可進行遞推。前向Euler法與精確值比較前向Euler法、后向Euler法與精確值比較前向Euler法在不同步長的結(jié)果比較二、梯形法 Euler法的計算精度較差，如果改用梯形面積代替每

6、個步距的曲線面積，就可提高精度。精確積分應為曲邊梯形的面積：現(xiàn)用直邊梯形的面積來近似：寫成差分方程為梯形近似及其誤差梯形法實質(zhì)是采用了兩點斜率平均值的結(jié)果，由于利用了兩點的信息，從而提高了計算精度。和 這一思想被廣泛地應用于許多算法中，實際計算時，如果在每個積分步矩中多取幾個點，分別求出其斜率，然后取不同的權(quán)值為： 后面RungeKutta法就是采用這樣的思想來進行計算的。梯形法的幾何意義也可按折線理解， 梯形法大大提高了精度，但為隱式算法，每次遞推計算時需解一次非線性方程，計算量較大由此考慮進行改進，先用Euler法計算出：的近似值代入導函數(shù)求出近似值再代入梯形公式求解。預估公式（Euler

7、法）校正公式（梯形法） 為預估校正法，也稱為改進的Euler法。例：設系統(tǒng)方程試用梯形法求其數(shù)值解，取步長h=0.1,解:梯形法遞推式：隱式算法，需先解此非線性方程：由此公式和初始條件可進行遞推，見仿真結(jié)果。例：設系統(tǒng)方程用改進歐拉法求數(shù)值解，取步長h=0.1,解:改進歐拉法：由此公式和初始條件可進行遞推，見仿真結(jié)果。前向Euler法、梯形法與精確值比較前向Euler法、改進歐拉法與精確值比較梯形法、改進Euler法與精確值比較三、Runge-Kutta法 、Taylor級數(shù)匹配原理由于輸入函數(shù)是 t 的函數(shù)，則將記做 得微分方程：如果對變量t、x具有各階導數(shù)，可推得x(t)的各階導數(shù)。 設已

8、知 進行Taylor級數(shù)展開： 若已知 的值，則當 h 較小時，可用級數(shù)展開的前 p+1項作為近似,令 則 即 以上公式(1)就稱為p階的Taylor展開法遞推公式 之間的誤差為： 局部截斷誤差與hp+1是同階無窮小量,記為O(hp+1) (1)歐拉法的Taylor級數(shù)展開只取一次項，其余忽略寫成差分方程為這就是解微分方程初值問題的歐拉算法。所以歐拉法稱為一階的Taylor展開法遞推公式 局部截斷誤差與h2是同階無窮小量,記為O(h2)梯形法的Taylor級數(shù)展開取一次項和二次項，寫成差分方程為所以梯形法稱為二階的Taylor展開法遞推公式O(h3) 可見Taylor展開法需用 在 的高階導數(shù)

9、計算 ，不便于數(shù)值計算。2、 Runge-Kutta法 于是Runge-Kutta法用 在一些點上的值表示 ，使局部截斷誤差的階數(shù)與Taylor展開法相等。避免了求高階導數(shù)，又保證高的精度。 對微分方程在區(qū)間的連續(xù)解為：在區(qū)間 取m 個點若已知則用它們的一次組合去近似 ，即？現(xiàn)在的問題是如何求設已知因為 是未知的，最簡單的用歐拉法構(gòu)造。由Euler法若按原有歐拉法計算下一步距的值得以二階Runge-Kutta法為例說明。精度較差，為改進精度，由歐拉法以步長 ah得到另一個點設已知某步從點開始，沿斜率為方向移動一個步長h，得到從圖中可見三點可得到比歐拉法精度較好的近似值。在一條直線上，可選取參數(shù)

10、 b1、 b構(gòu)造：二階R-K構(gòu)造法tX(t)tmtm+1xEm+1x(tm+1)tm+ahxEm+ax(tm)歐拉法計算以步長ah得到另一個點從點開始，沿斜率為方向移動一個步長h，得到可選取參數(shù) b1、 b構(gòu)造：以得到比歐拉法精度較好的近似值。以上計算歸結(jié)為：記則？如何選取參數(shù)a、b1、b2，可獲得最高的精度將k2在處展開成Taylor級數(shù)見代入式子見若滿足則與Taylor展開式前三項相同局部截斷誤差三個未知數(shù)，兩個方程，有多組解：局部截斷誤差是h的三階無窮小量,比歐拉法的精度高一階若取a=1時，得到：改進的Euler公式若取時，得：修正的Euler公式或中點公式顯式p階Runge-Kutta

11、法的一般形式為：3、常用的Runge-Kutta法：） Kutta三階法） Heun三階法3） 經(jīng)典顯式四階Runge-Kutta法四、微分方程數(shù)值積分的矩陣分析法 前述的各類數(shù)值積分公式都以一階系統(tǒng)（單個的微分方程）進行討論，而實際工程中大量的仿真對象是高階系統(tǒng)，可用一階微分方程組來描述。此時，數(shù)值積分公式有相應的矩陣形式。矩陣形式的數(shù)值積分公式：） 歐拉法公式前向歐拉法公式后向歐拉法公式） 梯形法公式）二階龍格庫塔法公式改進的歐拉法公式，是預估校正公式。） 四階龍格庫塔法公式（RK4）對于 n 階系統(tǒng)，狀態(tài)向量 x 為 n 維，計算中每前進一步 h ，要計算 4n 個 kij 值，對狀態(tài)空

12、間表達式：此時，RK4公式的個 k 值：例：系統(tǒng)方程取步長 h=0.1,試用RK4法求t=0.1，0.2時的解解:將原系統(tǒng)方程化為狀態(tài)方程形式：見仿真結(jié)果作業(yè)：P1493.2習題3-2：已知用前向歐拉法、梯形法求其數(shù)值解，取步長h=0.1解:前向歐拉法遞推式：習題3-2：已知用前向歐拉法、梯形法求其數(shù)值解，取步長h=0.1解:梯形法遞推式：隱式算法，需先解此非線性方程：單步法的特點：計算 n+1 時刻的值 yn+1 時，只用到第 n 時刻的 yn 和 fn 。如果能利用多步計算信息（歷史時刻值），則可能既加快仿真速度又獲得較高的仿真精度，這就是構(gòu)造多步法的出發(fā)點。第三節(jié) 數(shù)值積分法的多步算法實

13、際在逐步遞推過程中，計算 yn+1 時已經(jīng)獲得一系列的近似值： 以及 。多步法中以 Adams 法最具代表性，應用最為普遍。 對一階連續(xù)系統(tǒng)：連續(xù)解為：現(xiàn)過三點按插值原理構(gòu)造一個多項式 來逼近函數(shù)對函數(shù),再對多項式積分近似積分 一、Adams算法得：多項式 中的系數(shù)由下決定： 拉格朗日插值公式令：同時考慮：因為有：進行變量替換：顯然，對多項式的積分計算很容易。微分方程連續(xù)解為：寫成差分方程：這就是顯式兩步二階Adams遞推式。顯式 Adams 算法的系數(shù)值顯式 Adams 算法的遞推公式為： bibob1b2b3b40113/2-1/2223/12-16/125/12325/24-59/243

14、7/24-9/2441901/720-2774/7202616/720-1274/720-19/720k隱式 Adams 算法的系數(shù)值隱式 Adams 算法的遞推公式為： bib-1b0b1b2b30111/21/225/128/12-1/1239/2419/24-5/241/244251/720646/720-264/720106/720-19/720k無論用顯式或隱式 k 階Adams 法求解微分方程初值問題數(shù)值，需要先知道 k+1 個初始值。例如三步Adams 法：于是初始值只能從初始條件得到 ，還需知道才能求出：需用單步法求出，才能使多步法的遞推計算能夠進行。為保證多步法的精度，注意選

15、擇相應精度的單步法計算初始值。例：設系統(tǒng)方程用顯式二階Adams法求解，取步長h=0.1,解:顯式二階Adams法：起步初始值由梯形公式求出：下面就可以用Adams公式進行遞推：有初始條件：可進行初值計算： 為計算 的值，用到 時刻以前的值來推導，可獲得性能更好的算法。由 Adams 法得出更為一般的形式：即： 再令 二、線性多步法解出： 為計算 的值，用到和相應的導數(shù)值而且公式關于 x 、 f 是線性的，稱為線性 k 步法來推導即： 當： 稱為顯式線性 k 步法當： 稱為隱式線性 k 步法 時的情形。是線性 k 步法再看前面的 當： 顯式 Adams 算法的遞推公式為：隱式 Adams 算法

16、的遞推公式為：？接下來的問題就是線性 k 步法如何確定其中的常數(shù) 應用 Taylor 級數(shù)匹配原理，使局部截斷誤差盡可能的小。將處展開成Taylor 級數(shù)，和再帶入誤差公式得：其中：對于若能選取 使：則此線性多步法為 p 階 k 步算法。3、構(gòu)造線性多步法 確定步數(shù) k ，由Taylor級數(shù)匹配原理的 得到方程求出待定系數(shù) ，得出盡可能高階的算法。以二步法為例說明構(gòu)造步驟：5個未知數(shù)，4個方程， 得： 得到一般的線性二步法形式：是三階二步法； 又可得：是四階二步法：稱為Milne算法隱式算法為三階隱式算法就是三階隱式Adams法?如何得到二步顯式算法得到二步顯式算法但這是一個數(shù)值不穩(wěn)定算法，計

17、算中出現(xiàn)的微小誤差會迅速增長。因此，二步顯式算法無法達到三階，最多達到二階。只用三個方程二階二步顯式算法的一般形式：二階二步Adams顯式算法為Adams法為向后微分公式顯然是隱式公式若要求向后微分公式為K階算法，則當K=1為向后Eular法K=1，2，.6的向后微分公式系數(shù)值 11-1122/31/3-4/3136/11-2/119/11-18/111412/253/25-16/2536/25-48/2515161k 計算 時，僅已知和相應的導數(shù)值 稱為線性多步法的起步值，需用單步法來求出。 若線性多步法是P階算法，計算起步值的算法不應低于P階，否則影響計算精度。仿真步長 h 的選取，是否影

18、響仿真結(jié)果？先看一個例子：第四節(jié) 數(shù)值積分法穩(wěn)定性分析系統(tǒng)穩(wěn)定，精確解：用前向歐拉法求解： ch3_6.m由此可見，仿真步長 h 的選取，會影響仿真結(jié)果。 用前向歐拉法求解，當 h 0.2時不穩(wěn)定，是由于步長太大，從而截斷誤差太大造成的。一、數(shù)值解法穩(wěn)定性的含義 數(shù)值解的穩(wěn)定性：在擾動（初始誤差、舍入誤差、截斷誤差）的影響下，計算過程中的累積誤差不會隨計算步數(shù)的增加而無限增長。 微分方程的數(shù)值積分方法，實質(zhì)是微分方程的差分化，然后從初始條件遞推迭代。不同數(shù)值解法對應著不同的差分方程，是否穩(wěn)定取決于該差分方程的特征根是否滿足穩(wěn)定性要求。（處于Z平面上以原點為圓心的單位圓內(nèi)）著重研究單步法的穩(wěn)定性對步長的限制。二、數(shù)值解法穩(wěn)定性分析這樣做的根據(jù)是：1）試驗模型簡單，對其數(shù)值不穩(wěn)定的方法，不可用；2）一般的初始問題在其解的存在