計算機模擬技術

1、計算機模擬技術課程名：計算機模擬技術計算機模擬是在科學研究中常采用的一種技術，特別是在科學試驗環(huán)節(jié)，利用計算機模擬非常有效。所謂計算機模擬就是用計算機來模仿真實的事物，用一個模型（物理的-實物模擬；數學的-計算機模擬）來模擬真實的系統(tǒng)，對系統(tǒng)的內部結構、外界影響、功能、行為等進行實驗，通過實驗使系統(tǒng)達到優(yōu)良的性能，從而獲得良好的經濟效益和社會效益。計算機模擬方面的研究始于六十年代，早期的研究主要用于國防和軍事領域（如航空航天、武器研制、核試驗等），以及自動控制等方面。 隨著計算機應用的普及， 應用范圍也在擴大， 現在已遍及自然科學和社 會科學的各個領域。在農業(yè)方面，我國從80年代開始進行作物生

2、長發(fā)育模擬模型和生產管理系統(tǒng)的研究，目前有一定基礎的：在小麥方面有北農大、中科院；棉花方面有中國農業(yè)大學、中國棉花所；水稻方面有 江西農科院；在土壤水份、水資源及灌溉方面西北農業(yè)科技大學。目前影響較大的有比較成形的有江蘇省 農科院。目前的主要成果有：我國主要農作物栽培模擬優(yōu)化決策系統(tǒng)RCSODS（水稻）和WCSODS（小麥-江蘇省農科院）、MCSODS（玉米-河南省農科院）、CCSODS（棉花-中國農業(yè)大學）等。計算機模擬特別適合于實驗條件苛刻、環(huán)境惡劣（如真空、高溫、高壓、有毒有害的場所）、試驗周期長，花費大的場合。農作物的生產系統(tǒng)就很適合于計算機模擬：農作物的生產受各種條件的影響，不同作物

3、、不同品種也 有差異。比如，要想提高一種作物的產量，就先要作試驗，通過試驗了解這種作物的特性：抗旱性、耐寒 性、對氮、磷、鉀哪種肥更有效等。但農業(yè)的田間實驗不能保證精度（除人為可控條件外，還有許多隨機因素）、周期長（周期一年），耗費大?？赏ㄟ^計算機模擬來實現：先建立這種作物生產系統(tǒng)的數學模型 （依靠專業(yè)知識或試驗數據。一般來說，諸如作物產量和農業(yè)環(huán)境的關系可用微分方程或其它方程來描述），通過計算機模擬來找出這種作物的生長與農業(yè)環(huán)境相互作用的關系，以及各種條件之間的協(xié)迫情況。不僅可大大 節(jié)省實驗經費、加快研究進度 （周期一年的實驗結果幾秒鐘內即可得到），這種模擬軟件的開發(fā)還可與農業(yè)生產管理系統(tǒng)，

4、決策系統(tǒng)相聯(lián)系，實現對農作物生產的預測、分析、調控、設計的數字化和科學化。作為一門課程，不是研究某個特定系統(tǒng)的模擬問題，而是了解計握基礎知識，掌握建模及動態(tài)模擬的一般方法。第一章計算機模擬概述1.1計算機模擬技術研究對象在一個計算機模擬問題中，我們研究的對象是一個系統(tǒng)。系統(tǒng):一些具有特定功能的、相互之間按一定規(guī)律聯(lián)系著的實體的集合。如作物的生產系統(tǒng)可看作由作物、環(huán)境、技術、經濟等要素構成的。各要素之間相互影響、相互聯(lián)系，稱為系統(tǒng)的相關性；一個系統(tǒng)是一個整體，整體內的各個部分不能分割，各因素之間必須相互協(xié)調，不能在任何一個環(huán)節(jié)出問題，才能 使系統(tǒng)達到優(yōu)良的狀態(tài)，稱為系統(tǒng)的完整性。目標計算機模擬的

5、目標是了解系統(tǒng)的各個實體之間的相互制約關系，從而使系統(tǒng)在預定的目標下達到最優(yōu)和完善。 如在作物生產系統(tǒng)中，怎樣控制、實施各水、肥、栽培技術等，從而使產量最高，以 獲得最優(yōu)的經濟效益。方法模擬的方法是先建立系統(tǒng)與環(huán)境相互作用的數學模型，用數學模型來類比、模仿現實系統(tǒng)（一個數學模型就是從數學上表達系統(tǒng)各因素之間的數量關系，或各因素之間協(xié)調的規(guī)則；從整個模擬過程來看就是一個算法， 或一系列數據，這些數據綜合描述一個系統(tǒng)過程或現象的重要行為 ），然后在數學模型 和對系統(tǒng)深刻了解的基礎上，開發(fā)模擬軟件，用影響系統(tǒng)目標的因素作為輸入，通過計算機技術來表達系 統(tǒng)各因素作用的狀態(tài)。 從數學的角度來看，模擬的過

6、程就是對數學模型求解的過程，并把系統(tǒng)過程演示出 來。 基礎知識 可見，對一個系統(tǒng)進行計算機模擬，(1) 要對模擬對象有深刻的了解。 如：一個公交車調試系統(tǒng)，要編制一個好的調度程序，必須先對現 行系統(tǒng)作周密的調查，搞清哪些是影響調度程序的成分及實體，如現有車輛數、每車載客數、每趟車花費 的時間、沿途乘客的密集程度、乘客的一般去向、乘客高峰期的人數等。只有經過周密的調查研究，才能 形成一個完整的模擬系統(tǒng)； 作物生長模擬系統(tǒng)中， 也要搞清影響作物生長的各因素， 具備豐富的專業(yè)知識， 否則不能建立起精確的模型 。在對一個系統(tǒng)的動態(tài)特性不完全清楚的情況下，有必要通過實驗獲取數據， 以用于數學模型的建立。

7、(2) 要有數學知識。一般研究的系統(tǒng)較復雜，不能用簡單的函數或方程來描述，要綜合使用各種數學 方法，才能使模型準確、可靠。(3) 計算機知識和編程技巧。軟件應完整地實現系統(tǒng)的數學描述，輸出應直觀、形象，如三維可視化 輸出等。 軟件的開發(fā)是項目的重要工作， 也直接影響模擬的效果。 軟件的編寫可使用任何編程語言， 如 C、 VB、Java 等。專門的語言如 GPSS，GPSS 是面向對象問題的離散事件的專用模擬語言，優(yōu)其適用于排隊 系統(tǒng)。1961年，IBM公司發(fā)表GPSS的第一個版本，后來又有其它公司的各種版本。標準的版本有52個模塊， 每個模塊用特定的名稱和圖形來表示其功能 ?，F在一般常用的語言

8、都有模擬庫 (已編好的用于實現模 擬功能的函數 )專門用于模擬軟件的開發(fā)。1.2 系統(tǒng)的分類可模擬的系統(tǒng)各種各樣，不同類型的系統(tǒng)用不同的模型來描述。系統(tǒng)的分類方法很多，重要的分類方法是按系統(tǒng)的狀態(tài)是否隨時間變化 來分：一個行為與時間有關的系統(tǒng)稱為動態(tài)系統(tǒng)待研系統(tǒng)： 靜態(tài)系統(tǒng)： 系統(tǒng)的行為與時間無關； 用靜態(tài)模型來描述， 一般為數學方程、 邏輯表達式等。 如，電路的布爾表達式；電路中電壓與電流的關系；系統(tǒng)的穩(wěn)態(tài)解公式等 動態(tài)系統(tǒng)：連續(xù)系統(tǒng)： 系統(tǒng)的狀態(tài)隨時間連續(xù)變化；常用微分方程來描述；方程對所有 時間點有效。如，衛(wèi)星運行軌道，作物生長量等 確定性系統(tǒng)： 系統(tǒng)的輸出完全由其輸入來描述，即系統(tǒng)輸入

9、與 輸出按某種規(guī)則一一對應 集中參數模型： 用常微分方程來描述，即方程中的 導數不是偏導數分布參數模型： 用偏微分方程來描述，但一般用 集中參數模型近擬表示隨機系統(tǒng)： 系統(tǒng)的輸出是隨機的，有規(guī)律的存在一族隨機變量， 且隨機變量序列與時間有關 (隨機過程 ) 。 在確定性系統(tǒng)的模擬中使用隨機變量的研究方法 稱為蒙特卡羅方法離散系統(tǒng)： 系統(tǒng)狀態(tài)的變化只在離散的時間發(fā)生；動態(tài)方程只在離散點上有效； 一般為隨機系統(tǒng)。如，庫存問題；企業(yè)的管理系統(tǒng)等 離散時間系統(tǒng)： 時間步長固定；常用差分方程來描述 離散事件系統(tǒng)： 用事件來表示系統(tǒng)在時間間隔內的變化；常用 概率模型來描述1.3 建立數學模型對一個系統(tǒng)，確

10、定其類型后有助于選擇合適的方法建模，建模的一般步驟可分兩個大的階段：(1) 實質內容模型階段：首先對模擬對象進行調查(了解系統(tǒng)，搜索模擬所需的信息 )、實驗 (參數估計等統(tǒng)計推斷方法，確定參數及參數的敏感性 )、分析(將信息分類、量化，確定描述系統(tǒng)的規(guī)則 )，盡可能全 面地掌握系統(tǒng)的基本特性、運動規(guī)律以及中間狀態(tài) (最好是有對系統(tǒng)有深入了解的專業(yè)人員參與 )，通過分 析和邏輯推理，揭示系統(tǒng)內的規(guī)律。(2) 形式數量階段：在調查、實驗、分析的基礎上，進一步揭示系統(tǒng)內部的數量關系，并對其進行數 學處理，即對系統(tǒng)用數學形式來描述：用變量描述系統(tǒng)狀態(tài)，用各種數學方程定量表示各變量之間的相互 聯(lián)系，用遞

11、歸方程描述系統(tǒng)狀態(tài)的發(fā)展趨勢。多數情況下，建立一個好的(能真實的描述系統(tǒng)， 有代表性，能準確的模擬系統(tǒng)的數量信息，與實際系統(tǒng)較吻合)數學模型不易，優(yōu)其是復雜多變的系統(tǒng)或系統(tǒng)本身的特性尚不完全清楚的情況(如對農作物開發(fā)新品種)。所以在數學模型建立后，必須進行模擬驗證，如與真實系統(tǒng)相差較大， 則要重新建?；蛐薷哪Ｐ汀Ｋ砸粋€好的數學模型必須經過多次模擬，不斷修改、完善，才能得到。一個數學模型是描述系統(tǒng)行為的一個算法或一系列方程，工程中對系統(tǒng)建模應先建立系統(tǒng)的需求規(guī)格 說明，在制定模擬規(guī)劃之前予以充分討論，過程中需考慮以下因素：(1) 模型中需考慮哪些有效的因素：一個系統(tǒng)可能有不同的行為(如一個作物

12、生產系統(tǒng)中有作物的長勢、產量、質量、病蟲害等)，但不一定所有這些因素都要建立在模型中。模型中需要考慮的因素是：真正 能簡化系統(tǒng)的建模；對系統(tǒng)易于建模、測試和維護；使用較少的計算資源；對研究的系統(tǒng)有直接作用的。(2) 建模細化到什么程度：根據系統(tǒng)的需求確定細化的程度，是只須建立一個簡單的模型，還是要對 系統(tǒng)行為精確描述。在模型的準確性和花費之間求得平衡。(3) 與系統(tǒng)有相互作用的哪些外部環(huán)境 考慮在建模中：如在作物生長模型中，氣候、水、肥等。(4) 在建模中采取什么技術：首先，是基于物理的方程，還是基于測試數據。如果基于物理的方程， 是用微分方程，還是差分方程，考慮不考慮隨機因素等。這個問題常由

13、專業(yè)人員根據系統(tǒng)的專業(yè)知識來決 定；如果基于試驗數據，則建立經驗方程。(5) 建模時必須獲得什么樣的數據：如林木生長模型中，需要胸徑、樹高、材積等數據。另外還需考 慮這些數據的方便的輸出形式，以及模型需要的計算資源，如模型需要占用的內存、磁盤空間、消耗的 CPU時間等。(6) 建模和測試模型需多長時間，多少人力、物力、財力：隨著模型復雜性的增加，成本也會增加。 一般可從簡單開始，隨著應用逐步完善。(7) 如何驗證和確認模型：必須確認建立的模型被正確實現，模型所描述的行為與真實系統(tǒng)匹配到可 接收的程度，才能有價值。那么以什么標準來衡量?，F在已有驗證、確認與簽定(W&A)技術來證明和驗證模

14、擬的精確性。以上問題應列在系統(tǒng)需求說明書中，并應用于最高層次：在實際的模擬問題中，可能將系統(tǒng)分解為子 系統(tǒng)、組件，在對各組件、子系統(tǒng)建模時仍須遵從以上規(guī)則。模擬步驟：圖1-1迭代開發(fā)模型(1) 明確系統(tǒng)(2) 建立模型(3) 模型變換(4) 軟件設計開發(fā)(5) 測試檢驗(6) 評估、對結果的評價和分析一個模擬項目中各項工作的過程應該是一 個迭代過程，如圖1-1所示。下面通過實例來說明模擬過程。1.4應用舉例假設兩機相距10公里以內可實施攻擊,1.4.1兩物體追逐問題。設有一架殲擊機追蹤一架敵方轟炸機， 且須在12分鐘內完成追擊任務， 否則認為追擊失敗。 設兩機初始位置如圖 1-2所示。問題：對

15、轟炸機的任 一條特定航線，模擬殲擊機的追擊過程。圖1-2兩機追蹤模擬分析：在這個系統(tǒng)中，是對指定的轟炸機的一條航線而言，模擬殲擊機的追擊情況：殲擊機按什么航 線飛行，何時完成追擊任務。能否殲滅敵機由在規(guī)定的時間內兩機隨時間變動的距離而定，所以實施功擊 的時間、距離是要輸出的數據。在時刻t兩機距離由殲擊機在t時的位置決定，而t時的位置依賴于其速度和航向。為使模型簡單，作如下簡化：（1）設兩機在同一平面飛行。于是三維問題轉化為二維問題。（2）設殲擊機的速度（必須考慮的因素）Vf是常數（20km/mim）。變速須用微分方程描述，而常速即可用 一般方程表達，求解更簡單。（3）設殲擊機的航向（必須考慮的

16、因素）在厶t （設1分鐘）改變一次，而在1分鐘以內操作不變。這樣曲 線運動即變?yōu)檎劬€運動。（4）轟炸機航向（航線）可任意，比如現給定一條航線，如表1-1所示。表1-1轟炸機航線t012345678901112Xb (t)809099108116125133141151160169179180Yb (t)50484541353227212225293033 殲擊機初始位置：Xf（O）=O, Yf （0）=0（初始條件）建立數學模型：設在t時刻兩機位置為（XF（t）,YF（t）、（XB（t）,YB（t），兩機連線與水平線夾角為則1分鐘后殲擊機 的位置為：XF(t+1) = X F(t) +Vf C

17、osB(1.1)YF(t+1) = Y F(t) +VFSinB(1.2)由(XF(t),YF(t)計算(XF(t+1),YF(t+1)時涉及角0，從圖中看出：Sin 0 = Y b (t) - Yf (t) / D(t)(1.3)Cos0 = X b (t) - X f (t) / D(t)(1.4)而2 2D(t) = SQRY B(t) - Y f (t) + X b (t) - X f (t) (1.5)式（1.5） - （1.1）即是這個問題的數學模型。先算出兩機之距離， 不斷判斷是否在12分鐘內到達可攻擊的距離之內。程序流程圖如圖1-3所示。此例是連續(xù)系統(tǒng)，確定性模型，用一組方程來

18、描述。下例是一個離散事件系統(tǒng)，是隨機系統(tǒng)，用概率 模型來描述。開始VB程序見實例。142庫存問題。商業(yè)部門為了合理利用有限的流動資金，每項商品都要在庫存與銷售之間求得平衡: 庫存量太大會增加管理費用、積壓資金；庫存量太小又可能造成缺貨，也會造成銷售損失、信譽損失。所 以，當庫存量不滿足某一時段的顧客需求時， 就要到廠家訂貨。 這就需要采取一種策略： 當庫存量（比如布 匹）降到P匹布時（稱為定貨點）則向廠家訂貨，定貨量為 Q匹（稱為定貨量）。假設現在有5種策略（方案）。從中選擇一種使花費最少。策略PQ11251502125250315025041752505175300表1-2庫存策略已知條件：

19、（1）從發(fā)出訂貨單到收到貨物需3天，即第i天訂貨，第i+3天收到。（2）每匹布的保管費為0.75元，缺貨損失為1.8元/匹，訂貨費（包括手續(xù)費、采購差旅費及其他費用） 為750元。 需求量為一個0-99之間的均勻分布隨機數。（4）原始庫存量為115匹，并設第一天沒發(fā)出訂貨。0w x< 99分析：庫存問題是商業(yè)上的一個重要問題在數學上有專門的模型研究，存儲論也是運籌學的一個重要 分支。庫存問題用計算機模擬最合適，若通過銷售來找最優(yōu)方案，勢必造成經濟損失。這是典型的離散事 件系統(tǒng)，用概率模型來描述。這里已給定概率模型：X U : 0,99:即密度函數為：p(x) = 1/99分布函數為F(x

20、) = : 0 / 99以下框圖對每種策略模擬180天,x<00 < x w 99x>99選出效益最好的方案，如圖1-5所示。圖1-5庫存問題模擬程序框圖C程序見實例。以上兩個例子一個連續(xù)系統(tǒng)，狀態(tài)隨時間連續(xù)變化，用方程來描述；一個離散系統(tǒng)，到貨和銷售都按 一些離散的步驟進行，存在隨機性，只能用概率模型來描述。但解決問題的過程、分析方法類似：先分析 系統(tǒng)，使對系統(tǒng)有充分的了解；建立數學模型，設定一些初始條件（如t = 0時的狀態(tài)，開始時的庫存量）；系統(tǒng)狀態(tài)的變化對應一組方程或一組規(guī)則，隨著時是的變化，系統(tǒng)狀態(tài)改變，當一個周期結束時，收集模 擬過程的統(tǒng)計數據 （即問題的解）。如

21、果程序設計的好，就會使模擬非常逼真，就象一個真正的系統(tǒng)在演示 一樣。從理論上講，任何問題都可用計算機模擬。計算機模擬具有經濟、安全可靠、周期短等優(yōu)點。對任何 問題，只要建立起數學模型，改變參數值及變量值，就可模擬各種情況下的系統(tǒng)運行情況，從模擬輸出的結果，可分析系統(tǒng)內各因素的權重及其制約關系，幫助決策者作出合理的決策，克服盲目性，使系統(tǒng)在實 際運行過程中取得最好的效益。143傳染病傳播問題。假設某一地區(qū)有一種傳染病正在流行，那么政府、醫(yī)務部門就要采取措施來 控制這種疾病的傳播，要使采取的措施能夠有效的控制傳染速度.，就必須先搞清被傳染的人數跟哪些因素 有關，被傳染的人數是一個什么樣的發(fā)展趨勢，

22、從而來有效的預測和控制，下面建立描述被傳染人數的數 學模型。傳染病的傳播涉及因素很多，不可能通過一兩次簡單的假設就能建好完善的數學模型。這里的作法是:先作出最簡單的假設，看會得到什么樣的結果，然后對不合理的地方再行修改，逐步得到較滿意的模型。先討論一個粗略的模型。模型I:假設：(1) 每個病人在單位時間內傳染其他人的人數為一個常數ko。(2) 人得病后，經久不愈，該人在傳染期內不會死亡。記時刻t的得病人數為y(t)，開始模擬時有yo個傳染病人，則在 t時間內增加的病人人數為y(t+ t) - y(t) = k o y(t) t由導數定義得(在假設(1)、(2)下的數學模型)：-dy/dt =

23、ko y(t)(1.6)y(0) = yo初值問題(1.6)的解為：y(t) = y o eko t(分析：)這個結果表明，得病人數將按指數形式無限增加，當t is時，y(t) is,顯然與實際不符，說明上面的假設條件不合理。事實上，一個地區(qū)的總人數大致可視為常數(不考慮疾病傳播期間出生的、遷出的、死亡的)，所以一個病人在單位時間內能傳播的人數ko是在改變的：在初期，ko較大，隨著病人的增多，健康人數的減少，被傳染的機會也將減少，所以ko將變小。對上述模型進行改進：記 t時刻的健康人數人 S(t)，當總人數不變時，ko隨S(t)減少而減少。模型II :假設：(1) 總人數為常數n, t時刻的健

24、康人數為S(t)，得病人數為Y(t)，則Y(t) + S(t) = n(2) 單位時間內一個病人傳染的人數與當時的健康人數成正比，比例系數為(3) 同模型I的假設。由假設，(1.6)式中的ko應為k S (t)，即：/ dy(t) / dt = k S(t) Y(t)、y(o)= y o將(1.7)式代入(1.8)式，得(上述假設下的數學模型)：/ dy(t) / dt = k Y(t) (n-Y(t) y(o) = y o初值問題(1.9)的解為：-k n tY(t) = n / 1+(n / y o-1) e (分析：)由(1.1o)式得：2-k n t-k n t 2dy/dt = kn

25、 (n / yo-1) e / 1 + ( n / y o-1 ) e (1.7)k(醫(yī)學上稱為傳染強度.)(1.8)(1.9)(1.10)(1.11)(1.1o)式是被傳染人數隨時間變化的關系；(1.11)式是被傳染病人的變化率與時間的關系，如圖1-6和圖1-7所示。圖1-6病人人數變化曲線圖1-7病人變化曲線這個模型可預報疾病高峰到來的時間：令d2y / dt2 = 0，得極大值點：ti = ln ( n / y 0 - 1 ) / k n(1.12)由(12)式可知，當傳染強度k或人數n較大時，t1變化較小，表明傳染高峰到來的快，這與實際情況相吻合。但由(1.10)式知當tis時，y(t

26、) tn,這意味著最終人人都被傳染，這與實際不符，原因在于假設(3)不合理。再改進：可將人員分為三類：第一類為傳染者(y);第二類為易受傳染者(s),即這一類是非傳染者，但能得病而成為傳染者；第三類為除前兩者之外的人(r),包括患病死去的、痊愈后具有長期免疫力的、被隔離的。用y(t)、s(t)、r(t)分別表示這三類人數。模型III :假設：(1) 總人數為常數n,則：y(t) + s(t) + r(t) = n(1.13)(2) 同模型II的假設(2)(3) 單位時間內病愈免疫的人數與當時病人的人數成正比，比例系數為m (恢復系數)由假設，有dr / dt = m y(t)(1.14)由于引

27、入了 r(t)，則模型II的方程(1.8)應改為：dy / dt = k S (t) y (t) dr / dt(1.15)(1.15)式表示單位時間內病人人數的增加應等于被傳染的人數減去病愈的人數。從(1.13) - (1.15)中消去dy,并設 S (0) = So, r (0) = ro 得dS/dt = - k S(t) y(t) ，dr/dt = m y(t)y(t) + s(t) + r(t) = n (1.16)S(0) = S0r (0) = r 0'模型(1.16)較好地描述了傳染病傳播問題。通過以上實例，對數學模型的建立就有了一個基本的思路，對計算機模擬技術、模擬的

28、過程、問題的 方法步驟也有了一個概括的了解。后面兩章下分別對連續(xù)系統(tǒng)和離散系統(tǒng)討論基本的模擬方法。習題：(1)對例1.4.1中的飛機追擊問題采用極坐標(r, 0 )，相應的方程為：dr/dt = Vb Cos 0 -Vfr d 0 /dt = -Vb Sin 0其中r為兩機距離，0為兩機連線的角度。編程模擬兩機追擊情況。(2)球擺：如圖1-8所示，設繩長為L,夾角為0，球質量為m,初始速 度0 (0) = 0 ,初始偏角0 0，確定擺動周期。提示：影響擺動周期的因素： 擺球重力； 擺球質量； 擺球尺寸相對于繩很小，故可視為質點建模； 繩子質量忽略； 磨擦力忽略； 空氣阻力不考慮。建模：影響擺球

29、運動的重力 (使擺球回到中心位置的力)F的方向與繩子垂直；對繩子產生拉緊力的重 力與繩子平行，它不影響擺球運動，忽略。F = -mg Sin 0，將牛頓定律 F = ma代入上式得：a = -g Sin 0 ,其中a是擺球的切線加速度，由于L = 0 "，故得運動方程：0 " = -L/g Sin 00 (0) = 0 o ,0 '(0) = 0第二章連續(xù)系統(tǒng)的計算機模擬技術連續(xù)系統(tǒng)的狀態(tài)隨時間連續(xù)的動狀變化，這種變化依賴相關的因素，且遵從一定的規(guī)律，所以一般來 說可用數學方程描述(代數方程、微分方程，此外還有傳遞函數、結構圖及狀態(tài)方程等)。問題：對數學模型怎樣求解

30、,求解過程也就是模擬的過程，即程序的算法。對前面的飛機追擊的問題，其模型是一組代數方程，而且是顯式.的，所以，作法是：假設 1分鐘改變一次航向，每隔 1分鐘計算一次系統(tǒng)的狀態(tài)，而且 由模型知道了殲擊機前一分鐘的狀態(tài)就可以算出后一分鐘的位置。這樣，通過迭代，就可求出問題的解。對一般的模擬問題，模型用微分方程表達(方程中除含有自變量外， 還有導數或偏導數 一分布參數型可 通過一些變換轉換為集中參數型)，而對一個常微分方程(只有線性的或幾種特殊的能求出通解及特解，即 解析解)一般來說求解析解是不可能的，(即使能求出解析解)在計算機上求解要用 數值積分法:將連續(xù)的系統(tǒng)離散化，把微分方程轉化為差分方程，

31、通過迭代運算，求出問題的數值解。通過正確的控制步長、選擇 合理的算法即能達到要求的精度，這就是連續(xù)系統(tǒng)模擬的主要技術。數值積分法種類很多，本章介紹幾種簡單常用的方法。2.1 歐拉(Euler)法設連續(xù)系統(tǒng)的數學模型為dy/dt = f (t,y)(2.1)耳L y(0) = y o為了模擬系統(tǒng)狀態(tài)y隨時間t的變化，需求解微分方程(2.1)的數值解(不是解析解)，為此，把模擬周期分為 若干小區(qū)間，比如分為N個相等的小區(qū)間，如圖 2-1所示。只須在每個離散的點上計算系統(tǒng)的狀態(tài)。在區(qū)間(tn, tn+1)上求積分，得tn "1積分的幾何意義是小曲邊梯形的面積。如果小區(qū)間取的足夠小，則在區(qū)間

32、 面積代替小曲邊梯形的面積，于是得在y(t n+1) - y(t n) =f (t, y)dt(tn, tn+1)之間的f (t,y)可近似的看成常數f(tn,yn)，這樣可用小矩形的 tn+1時的積分值為y(tn+1 )yn + f (tn,yn) h = y n+1(2.2)其中h = T / N，將(2.2)式寫為差分方程形式，得yn+1 = yn + f (tn,yn) hn = 0,1,2,(2.3)這就是歐拉公式：它是一個遞推的差分方程，任何一個新的數值解 yn+1均基于前一個數值解 yn以及導數值f(tn,yn)求得，只要給出初始條件 y0及步長h,就可算出y1,由y1可算出 北

33、,如此迭代計算 % , 直到滿足所需計算的范圍。歐拉法也叫折線法，特點是方法簡單，計算量小，但計算精度也底。擴展I:改進的歐拉法(梯形法)為了提高計算的精度，可對歐拉法進行改進：從幾何意義上看到，用小梯形的面積來代替小矩形的面 積，必能提高精度。這樣(2.3)式即可寫成：yn+1 = yn + h / 2 f (tn,yn) + f (tn+1,yn+1) =yn + h / 2 f n + f n+1(2.4)但(2.4)式是隱式公式，即公式右端含有yn+1，而這是未知的待求量，故梯形法不能自行啟動運算，要借助于其它算法：如用歐拉法算出yPn+1作為梯形法中yCn+1的預估值，再進行計算，這

34、樣，公式可寫為：預估：yPn+1=yn +f (tn,yn) h校正：yCn+1=yn +h / 2 f (t n,y n)+ f (t n+1,yPn+1)卜(2.5)=yn + h / 2 f n + f Pn+l(2.5)式也稱為預估校正法，計算量稍大(需要附加的計算)，但精度要比歐拉法高，穩(wěn)定性好。擴展II :多個輸入的多階系統(tǒng)在(2.1)式所表示的數學模型中，y是系統(tǒng)的狀態(tài)，稱之為狀態(tài)變量，它隨時間t而連續(xù)變化。(2.1)式是 最簡單的情況，只有一個狀態(tài)變量，而它直接依賴時間而變化，在實際模擬中，狀態(tài)變量可能不只一個， 而影響系統(tǒng)狀態(tài)的因素更多，這些因素是隨時間變化的同，而系統(tǒng)狀態(tài)的

35、改變也正是由于這些因素隨時間 變化而變化。如在林木的生長系統(tǒng)中，要考察的狀態(tài)可能有樹高、胸徑、材積等；在作物的生長系統(tǒng)中， 考察的狀態(tài)可能有株高、葉子的片數、葉子的寬度等，而影響系統(tǒng)狀態(tài)的因素如水、肥、氣候、土質、病 蟲害等，而這些因素都隨時間而變化，如水份會隨時間而被蒸發(fā)，肥會隨時間被作物吸收，氣候等自然環(huán) 境更是隨時間變化無常。所以與(2.1)式相應的數學模型應為：dyk / dt = f (xi(t), X2(t),xm(t), y i, y2,yq)k = 1,2,q(2.6)y k (0) = y k, o其中y k為系統(tǒng)狀態(tài)變量，Xi為系統(tǒng)輸入信號，這樣的系統(tǒng)稱為具有 m個輸入、q

36、階系統(tǒng)。(2.6)式是q個微 分方程的方程組，q個初始條件?？梢园?2.6)式表示成向量的形式，就和(2.1)式在形式上一致了：記Y為q個狀態(tài)變量和向量，X為m個輸入信號向量，即：Y = (y1(t), y2(t),q yt)，X = (X1 (t), X2 (t),應(t)則(2.6)式寫為：” dY/dt = f (X(t), Y)(2.6)'L Y(0) = Yo對于(2.6)式所表達的數學模型，與(2.3)式相對應的歐拉公式為：| y k, n+1 = y k, n + f(X 1(t n), X 2 (t n),X m(t n), y1, n, y2, n,yq, n) h(

37、2.7)1 k = 1,2,- ,q或寫成向量的形式：Yn+1 = Yn+f (X n , Y n) h(2.7)'例2.1設兩種物質A和B合到一起產生化學反應，生成新物質C，假設1克A和1克B結合能產生2克C，形成C的速率與A和B的數量乘積成正比。同樣C也可分解為 A和B，C分解的速率正比于C的數量。問題：在給定 A和B的數量后，模擬有多少C物質產生出來，以及達到穩(wěn)定的時間。解：在任何時刻，設a，b，c分別是A，B，C的數量，則它們增加和減少的速度服從下列微分方程：'da / dt = k2 C - k 1 a b< db /dt = k2 C - k1 a b(2.8

38、)、dc / dt = 2 k 1 a b - 2 k2 C其中k1和k2是比例常數(實際問題中k1和k2會隨溫度、壓力等發(fā)生變化，但在模擬過程中為簡化模型， 可 視為常數)，是模型中的參數，-k1 ab中的負號是因為 A和B是減少的過程。假設模擬從 t = to(一般取0) 開始，使t以厶t時間間隔增加(由步長厶t和模擬周期可定出所分時間段數 )，則相應的歐拉公式為(2.9)式。 給出常數k1和k2值以及A、B、C的初始數量，即得迭代公式：廠a n+1 = a n + (k2 Cn - k1 an bn) tb n+1 = b n + (k2 Cn - k 1 an bn) t4C n+1

39、= C n + (2k 1 an g - 2k2 Cn) t(2.9)a (0) = a 0 , b (0) = b o , c (0) = 0l k1 = k 1,0 , k2 = k 2, 0由(2.9)式，從 t = 0 開始，由 a 0、b°、C0 可算出 a 1、b1、c 1，又可算出 a 2、b2、C2, (a 2 = a(2 t)圖2-2例2-1模擬流程圖以厶t為間隔，進行N = T/次計算，就可算出周期T的系統(tǒng)狀態(tài)，從而得出模擬結果。根據數學模型(2.9)就可設計編寫模擬程序，可選用任何編程語言，在一些流行的語言中，對常用的模 擬算法(比如后面要介紹的龍格-庫塔(Ru

40、nge-kutta)方法)都有相應的用于模擬計算的軟件包。public class Ratestatic double k1,k2,h,t;開始static double A=new double53;static double B=new double53;static double C=new double53; static int i;public static void main(String args) A1=100.0;B1=50.0;C1=0.0; t=0;h=0.1;k1=0.008;k2=0.002; for(i=1;i<53;i+)System.out.pri nt

41、l n(i+" "+t+", "+Ai+","+Bi+","+Ci);strut(i); static void strut( int i)Ai+1=Ai+(k2*Ci-k1*Ai*Bi)*h; Bi+1=Bi+(k2*Ci-k1*Ai*Bi)*h;Ci+1=Ci+2.0*(k1*Ai*Bi-k2*Ci)*h; t=t+h;歐拉法是最簡單的數值積分法，在介紹更好的數值積分法之前，先介紹幾個關于數值積分的基本概念。1、單步法與多步法數值積分法是用遞推公式求解，如果僅由前一時刻的數值yn就能算出后一時刻的數值yn+1，

42、則稱為單步法，反之，如果求yn+1時需用到y(tǒng)n, yn-1 , yn-2 ,等多個值，則稱為多步法.。如，歐拉法是單步法，而 改進和歐拉法是多步法。單步法由遞推公式自身就能啟動運算(由初值即可算出y1，由y1可算出y2,)，所以它是能自啟動的算法；多步法在開始的時候要先用其它的方法計算該時刻前面的函數值，所以不能自啟動運算。一般來說，由于多步法更能充分利用多個時刻的信息，所以模擬速度快，精度高，但計算量要大一些，后面還要介紹多步法的算法。2、顯式與隱式在遞推公式中，計算yn+1時所用到的數據均已算出的計算公式稱為顯式公式.，相反，在算式中隱含未知量yn+1稱為隱式公式。如歐拉法中遞推公式是顯式

43、的，而梯形法中遞推公式是隱式的。在隱式公式中，必須先用顯式公式估計一個值，再用隱式公式迭代，即預估-校正法。隱式與顯式相比，有明顯的高精度和穩(wěn)定性。3、誤差在數值解法的過程中， 每一步計算都會產生誤差，誤差的來源有兩個方面： 一是計算誤差(即計算機計 算本身的誤差)，一個是公式誤差(用差分方程代替微分方程)。分別稱之為舍入誤差.和截斷誤差。舍入誤差：計算機的字長是有限的，數字不能表示的完全精確，實際上計算的結果是用有限精度的有 理數(如計算機中使用的浮點數)來近似無限精度的實數，所以在對動態(tài)系統(tǒng)模擬的過程中，舍入誤差是不 可避免的。通常舍入誤差的大小與積分步長成反比，步長h越小，計算次數多則舍

44、入誤差大，但不能隨意加大步長，否則將產生大的截斷誤差甚至影響系統(tǒng)穩(wěn)定性。在給定運算序列的條件下，唯一可減少舍入誤 差的方法是增加數字表示的精度，如將單精度浮點數(有7位數的精度)表示改為雙精度數(有15位數的精度)。實際的動態(tài)系統(tǒng)模擬時，多采用雙精度數據類型，盡管這樣變量占用的存儲空間大，處理時間稍長， 但對提高模擬的精度來講是值得的。截斷誤差：是用差分方程代替微分方程產生的誤差，即用數值解代替微分方程的精確解產生的誤差， 所以是公式本身的誤差，一般用臺勞級數來分析積分公式的精度：假設在t n時積分(精確值)已經算出，則用臺勞級數可求得t n+1時的精確解：y(t n+1)= y(t n +

45、h) = y(t n) + y ')(+ h / 2! y II (t) + , +h / r! y ( t n) + O(h )(2.10)如果一個數值解法是用前r+1項來近似的計算y(t n+1)，則后面的各項(記為)0(h r+1)是在這一步計算中引進的附加誤差，稱為這個算法的(局部)截斷誤差。誤差0(h r+1)與h r+1同階(ht0)，即計算的精度保特了r階,此時稱這個方法是r階的精度。一個數值方程的階數可視為衡量這一方法的精確度的重要標志，不同的數 值解法有不同的精度。歐拉法只取(2.10)中的前兩項作近似計算，所以是一階精度的算法；梯形法相當于取(2.10)式中的前三項

46、，故是二階精度的方法。4、穩(wěn)定性如果系統(tǒng)是穩(wěn)定的(系統(tǒng)的狀態(tài)隨時間的推移逐步穩(wěn)定在某個水平上)，則在模擬的迭代過程中，數值積分的解也應該是穩(wěn)定的，但由于初始數據的誤差及在迭代運算中的舍入誤差會對后面的計算結果產生影 響。當步長h選擇不合理時，可能使模擬結果不穩(wěn)定。對于歐拉法，可用下面的檢驗方程(其中入為方程的 特征根)：-y' =y(2.11)丫 y(0) = y 0(注意其精確解為y = yoe"t)來檢驗步長對數值解穩(wěn)定性的影響：對(2.11)表示的方程，歐拉公式為：y n+1 = y n + 入 y n h = (1+ 入 h)y n(2.12)(y(0) = y 0所

47、以，要使數值解穩(wěn)定，必須使：| 1+入h| < 1 ，解得|入h|<2或h<2T (T = 1/入是系統(tǒng)時間常數)。所以要保證歐拉方法計算的穩(wěn)定性，步長h必須小于系統(tǒng)時間常數的2倍。在(2.11)中特征根入在一定數量級范圍變動，現令 入=1，來作一個具體的模擬，此時 h應小于2，否 貝U將不穩(wěn)定。取 h = 0.01 , h = 0.05 , h = 1.0, h = 1.9 , h = 2.0 , h = 2.1六個不同的步長，yo = 1。模擬結果 為：h = 0.01 :表現出較好精度h = 0.05 :雖然近似解接近精確解，但存在誤差h = 1.0 :解在一步后趨于零

48、，并一直保持，雖然穩(wěn)定，但明顯精度不高h = 1.9 :解振蕩，幅度值逐衰減，并趨于穩(wěn)定h = 2.0 :解不衰減，等幅振蕩h = 2.1 :解的振蕩幅度值遞增，表明系統(tǒng)不穩(wěn)定，數值解發(fā)散 數值解圖如圖2-3所示。2.2 龍格-庫塔(Run ge-Kutta)法R-K方法的基本思想是用臺勞展開式的前幾項來對微分方程求近似解。再以模型(2.1)為例：y ' = f (t, y)(2.1)I y (t o) = yo假設從t o開始，以h增長，hi = t o + h，yi = y (t o + h)，在to附近展開成臺勞級數，保留h2嘰 則有:2yi = y o + f (t o , y

49、 o) h + h /2 ( S f / S y dy / dt + S f / S t)(2.13)(此式括號中的導數是在(t o , y o)處的導數值)為了求(2.13 )的解，假設(2.13 )的解可寫成如下的形式:yi = y o + (b1 k1 + b2 k2 ) h ,其中k1 = f (t o , y o),(2.14)k 2 = f (t o +C2 h , y o + a1 k1 h)(注意(2.13 )式是關于函數y的導數的算式，而(2.14)式中k1和k2都是y在某點處的導數，相差的只是常量級的系數不同，問題正是要定出這些系數，從而得到數值解表達式，為此)對k2式右端

50、的函數在(to , yo)處展開臺勞級數，保留 h項，得：k 2 f (t o , yo) + (C2 S f / S t + a1 k1 S f / S y dy / dt) h把燈、k2代入(2.14)中，得yi = y o + bif (t o , y 0) h + b2 h f (t 0 , yo) + (C2 S f / S t + ai ki S f / S y dy / dt) h (2.15)將(2.15)式與(2.14)式比較，得關于系數的方程：(2.16)bi = b2，得一個解:(2.17). bi + b2 = 1b2 C2 = 1/2-b2 ai = 1/2方程組(2

51、.16)中四個未知量，三個方程，故有無窮多解，求出一個解，比如令b i = b 2 = 1/2ai = C 2 = 1代入(2.14)，得如下一個公式：yi = y o + 1/2 (K i + K 2 ) h其中ki = f (t 0 , y 0)jI k 2 = f (t 0 + h , y 0 + ki h)這是從t 0計算t 1時刻的公式，寫成一般的形式，得如下遞推公式:其中yn+i = y n + h/2 (K1 + K2 ) ki = f (t n , y n)k 2 = f (t n + h , y n + ki h)模型(2.1)的數值解公式(2.17)即稱為R-K公式，這種數

52、值解法即稱為 R-K方法。在(2.13)式中，只保留了h2項，故公式(2.17)的精度是2階的，公式(2.17)稱為二階R-K公式。目前在實際模擬問題中，四階R-K公式用的最為普遍。在推導四階R-K公式時，在臺勞公式中保留到h4項，推導過程與前面類似。一個四階R-K公式可以是下面的形式：yn+i = y n + h/6 (K1 + 2K2 +2K3 +K4)飛其中f ki = f (t n , y n)k 2 = f (t n + h /2 , y n + ki h/2)卜(2.18)k 3 = f (t n + h /2 , y n + k2 h/2)x k 4 = f (t n + h ,

53、 y n + k3 h)-問題：關于R-K方法還有下面一些問題需了解。(1) 多種形式：方程組(2.16)有無窮多解，我們取了一個解，得到了公式 (2.17)，也可以取其它的解， 所以每一階R-K公式都有多種形式。二階R-K公式常用的除(2.17)式外，還有yn+1 = y n + K2 hki = f (t n , y n)' k 2 = f (t n + h /2 , y n + ki h/2)(對應(2.16)的解bi = 0 , b2 = 1 , C2 = a 1 = 1/2 )。四階R-K公式常用的除(2.18)式外，還有廠 yn+i = y n + h/8 (K1 + 3K

54、2 +3K3 +K4) ki = f (t n, yn)Y k 2 = f (t n+h /3 ,y n + ki h/3)k 3 = f (t n+h 2/3, y n + kih/3 +k2 h)'k 4 = f (t n+h , yn + hki -hk2 + hk3)(2) 精度：也可推導3階、5階或更高階的R-K公式，但在一般工程中，四階公式就完全能達到要求 的精度，而且四階公式是最常用的，所以一般不使用更高階的公式。另一個特殊情況：一階 R-K公式，只含有h的1次項，即：yn+i = y n + h f (t n , y n)，這就是歐拉公式, 所以歐拉公式也可看面一階R-

55、K公式，其精度最低。R-K方法的精度取決于步長 h及求解方法。一般來說，為達到同樣的精度，四階方法的步長可以比二 階方法的步長大10倍，而四階方法每步的計算量僅比二階方法大一倍，所以總的計算量仍比二階方法小。 R-K方法可使用較大的步長也是其一個特點。(3) 單步法：不論幾階的R-K公式都是單步法，在計算 yn+i時只用到y(tǒng) n，可自啟動，使用的存儲量也小。另外，無論幾階的R-K公式，算式中均有兩部分組成：一部分是上一步的結果y n，第二部分是h=t n -n-1中對f(t,y)的積分，它是步長h乘以各點斜率的加權平均。如在四階公式(2.18)中，取四點斜率ki、k?、k?、k4，對k2、k3各取兩份，而ki和k4各取一份進行加權平均。(4) 變步長的R-K方法：在前面討論的這些數值方法中，都需要在模擬之前選擇步長，且假設步長是 固定的。這樣作有一定的缺點：如果 h選的太小，會增加計算量，增長模擬的時間；而太大又會達不到精 度要求。另一種方法是在模擬過程中動態(tài)調整步長，這樣，當狀態(tài)變量變化緩慢時，步長可取大