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

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

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

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

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

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

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

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

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

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

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

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

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

13、際在逐步遞推過程中，計(jì)算 yn+1 時(shí)已經(jīng)獲得一系列的近似值： 以及 。多步法中以 Adams 法最具代表性，應(yīng)用最為普遍。 對(duì)一階連續(xù)系統(tǒng)：連續(xù)解為：現(xiàn)過三點(diǎn)按插值原理構(gòu)造一個(gè)多項(xiàng)式 來逼近函數(shù)對(duì)函數(shù),再對(duì)多項(xiàng)式積分近似積分 一、Adams算法得：多項(xiàng)式 中的系數(shù)由下決定： 拉格朗日插值公式令：同時(shí)考慮：因?yàn)橛校哼M(jìn)行變量替換：顯然，對(duì)多項(xiàng)式的積分計(jì)算很容易。微分方程連續(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 個(gè)初始值。例如三步Adams 法：于是初始值只能從初始條件得到 ，還需知道才能求出：需用單步法求出，才能使多步法的遞推計(jì)算能夠進(jìn)行。為保證多步法的精度，注意選

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

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

17、算中出現(xiàn)的微小誤差會(huì)迅速增長(zhǎng)。因此，二步顯式算法無法達(dá)到三階，最多達(dá)到二階。只用三個(gè)方程二階二步顯式算法的一般形式：二階二步Adams顯式算法為Adams法為向后微分公式顯然是隱式公式若要求向后微分公式為K階算法，則當(dāng)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 計(jì)算 時(shí)，僅已知和相應(yīng)的導(dǎo)數(shù)值 稱為線性多步法的起步值，需用單步法來求出。 若線性多步法是P階算法，計(jì)算起步值的算法不應(yīng)低于P階，否則影響計(jì)算精度。仿真步長(zhǎng) h 的選取，是否影

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