計(jì)算機(jī)模擬技術(shù)

計(jì)算機(jī)模擬技術(shù)

課程名：計(jì)算機(jī)模擬技術(shù)

讓算機(jī)摸擬是在科學(xué)研究中常采用的一種技術(shù)，特別是在科學(xué)試驗(yàn)環(huán)節(jié)，利用計(jì)算機(jī)模擬非常有效。

所謂計(jì)算機(jī)模擬就是用計(jì)算機(jī)來模仿真實(shí)的事物，用一個模型（物理的-實(shí)物模擬；數(shù)學(xué)的-計(jì)算機(jī)模擬）來

模擬真實(shí)的系統(tǒng)，對系統(tǒng)的內(nèi)部結(jié)構(gòu)、外界影響、功能、行為等進(jìn)行實(shí)驗(yàn)，通過實(shí)驗(yàn)使系統(tǒng)達(dá)到優(yōu)良的性

能，從而獲得良好的經(jīng)濟(jì)效益和社會效益。

計(jì)算機(jī)模擬方面的研究始于六十年代，早期的研究主要用于國防和軍事領(lǐng)域（如航空航天、武器研制、

核試驗(yàn)等），以及自動控制等方面。隨著計(jì)算機(jī)應(yīng)用的普及，應(yīng)用范圍也在擴(kuò)大，現(xiàn)在已遍及自然科學(xué)和社

會科學(xué)的各個領(lǐng)域。在農(nóng)業(yè)方面,我國從80年代開始進(jìn)行作物生長發(fā)育模擬模型和生產(chǎn)管理系統(tǒng)的研究，

目前有一定基礎(chǔ)的:在小麥方面有北農(nóng)大、中科院；棉花方面有中國農(nóng)業(yè)大學(xué)、中國棉花所；水稻方面有

江西農(nóng)科院；在土壤水份、水資源及灌溉方面西北農(nóng)業(yè)科技大學(xué)。目前影響較大的有比較成形的有江蘇省

農(nóng)科院。目前的主要成果有：我國主要農(nóng)作物栽培模擬優(yōu)化決策系統(tǒng)RCSODS（水稻）和WCSODS（小麥-江

蘇省農(nóng)科院）、MCSODS（玉米-河南省農(nóng)科院）、CCSODS（棉花-中國農(nóng)業(yè)大學(xué)）等。

計(jì)算機(jī)模擬特別適合于實(shí)驗(yàn)條件苛刻、環(huán)境惡劣（如真空、高溫、高壓、有毒有害的場所）、試驗(yàn)周期

長，花費(fèi)大的場合。

農(nóng)作物的生產(chǎn)系統(tǒng)就很適合于計(jì)算機(jī)模擬：農(nóng)作物的生產(chǎn)受各種條件的影響，不同作物、不同品種也

有差異。比如，要想提高一種作物的產(chǎn)量，就先要作試驗(yàn)，通過試驗(yàn)了解這種作物的特性：抗旱性、耐寒

性、對氮、磷、鉀哪種肥更有效等。但農(nóng)業(yè)的田間實(shí)驗(yàn)不能保證精度（除人為可控條件外，還有許多隨機(jī)因

素）、周期長（周期一年），耗費(fèi)大?？赏ㄟ^計(jì)算機(jī)模擬來實(shí)現(xiàn)：先建立這種作物生產(chǎn)系統(tǒng)的數(shù)學(xué)模型（依靠專

業(yè)知識或試驗(yàn)數(shù)據(jù)。一般來說，諸如作物產(chǎn)量和農(nóng)業(yè)環(huán)境的關(guān)系可用微分方程或其它方程來描述），通過計(jì)

算機(jī)模擬來找出這種作物的生長與農(nóng)業(yè)環(huán)境相互作用的關(guān)系，以及各種條件之間的協(xié)迫情況。不僅可大大

節(jié)省實(shí)驗(yàn)經(jīng)費(fèi)、加快研究進(jìn)度（周期一年的實(shí)險(xiǎn)結(jié)果幾秒鐘內(nèi)即可得到），這種模擬軟件的開發(fā)還可與農(nóng)業(yè)

生產(chǎn)管理系統(tǒng)，決策系統(tǒng)相聯(lián)系，實(shí)現(xiàn)對農(nóng)作物生產(chǎn)的預(yù)測、分析、調(diào)控、設(shè)計(jì)的數(shù)字化和科學(xué)化。

作為一門課程，不是研究某個特定系統(tǒng)的模擬問題，而是了解計(jì)算機(jī)模擬的一般過程、基本原則，掌

握基礎(chǔ)知識，掌握建模及動態(tài)模擬的一般方法。

第一章計(jì)算機(jī)模擬概述

1.1計(jì)算機(jī)模擬技術(shù)

?研究對象在一個計(jì)算機(jī)模擬問題中，我們研究的對象是一個羞纏。

系統(tǒng)：一些具有特定功能的、相互之間按一定規(guī)律聯(lián)系著的實(shí)體的集合。如作物的生產(chǎn)系統(tǒng)可看作由

作物、環(huán)境、技術(shù)、經(jīng)濟(jì)等要素構(gòu)成的。各要素之間相互影響、相互聯(lián)系，稱為系統(tǒng)的相關(guān)性；一個系統(tǒng)

是一個整體，整體內(nèi)的各個部分不能分割，各因素之間必須相互協(xié)調(diào)，不能在任何一個環(huán)節(jié)出問題，才能

使系統(tǒng)達(dá)到優(yōu)良的狀態(tài)，稱為系統(tǒng)的完整性。

?目標(biāo)計(jì)算機(jī)模擬的目標(biāo)是了解系統(tǒng)的各個實(shí)體之間的相互制約關(guān)系，從而使系統(tǒng)在預(yù)定的目標(biāo)

下達(dá)到最優(yōu)和完善。如在作物生產(chǎn)系統(tǒng)中，怎樣控制、實(shí)施各水、肥、栽培技術(shù)等，從而使產(chǎn)量最高，以

獲得最優(yōu)的經(jīng)濟(jì)效益。

?方法模擬的方法是先建立系統(tǒng)與環(huán)境相互作用的數(shù)學(xué)模型，用數(shù)學(xué)模型來類比、模仿現(xiàn)實(shí)系統(tǒng)

（一個數(shù)學(xué)模型就是從數(shù)學(xué)上表達(dá)系統(tǒng)各因素之間的數(shù)量關(guān)系，或各因素之間協(xié)調(diào)的規(guī)則；從整個模擬過程

來看就是一個算法，或一系列數(shù)據(jù)，這些數(shù)據(jù)綜合描述一個系統(tǒng)過程或現(xiàn)象的重要行為），然后在數(shù)學(xué)模型

和對系統(tǒng)深刻了解的基礎(chǔ)上，開發(fā)模擬軟件，用影響系統(tǒng)目標(biāo)的因素作為輸入，通過計(jì)算機(jī)技術(shù)來表達(dá)系

統(tǒng)各因素作用的狀態(tài)。從數(shù)學(xué)的角度來看，模擬的過程就是對數(shù)學(xué)模型求解的過程，并把系統(tǒng)過程演示出

來。

?基礎(chǔ)知識可見，對一個系統(tǒng)進(jìn)行計(jì)算機(jī)模擬，

(1)要對模擬對象有深刻的了解。如：一個公交車調(diào)試系統(tǒng)，要編制一個好的調(diào)度程序，必須先對現(xiàn)

行系統(tǒng)作周密的調(diào)查，搞清哪些是影響調(diào)度程序的成分及實(shí)體，如現(xiàn)有車輛數(shù)、每車載客數(shù)、每瓶車花費(fèi)

的時(shí)間、沿途乘客的密集程度、乘客的一般去向、乘客高峰期的人數(shù)等.只有經(jīng)過周密的調(diào)查研究，才能

形成一個完整的模擬系統(tǒng)；作物生長模擬系統(tǒng)中，也要搞清影響作物生長的各因素，具備豐富的專業(yè)知識，

否則不能建立起精確的模型。在對一個系統(tǒng)的動態(tài)特性不完全清楚的情況下，有必要通過實(shí)驗(yàn)獲取數(shù)據(jù)，

以用于數(shù)學(xué)模型的建立。

(2)要有數(shù)學(xué)知識。一般研究的系統(tǒng)較復(fù)雜，不能用簡單的函數(shù)或方程來描述，要綜合使用各種數(shù)學(xué)

方法，才能使模型準(zhǔn)確、可靠。

(3)計(jì)算機(jī)知識和編程技巧。軟件應(yīng)完整地實(shí)現(xiàn)系統(tǒng)的數(shù)學(xué)描述，輸出應(yīng)直觀、形象，如三維可視化

輸出等。軟件的開發(fā)是項(xiàng)目的重要工作，也直接影響模擬的效果。軟件的編寫可使用任何編程語言，如C、

VB、Java等。專門的語言如GPSS,GPSS是面向?qū)ο髥栴}的離散事件的專用模擬語言，優(yōu)其適用于排隊(duì)

系統(tǒng)。1961年，IBM公司發(fā)表GPSS的第一個版本，后來又有其它公司的各種版本。標(biāo)準(zhǔn)的版本有52個

模塊，每個模塊用特定的名稱和圖形來表示其功能。現(xiàn)在一般常用的語言都有模擬庫(已編好的用于實(shí)現(xiàn)模

擬功能的函數(shù))專門用于模擬軟件的開發(fā)。

1.2系統(tǒng)的分類

可模擬的系統(tǒng)各種各樣，不同類型的系統(tǒng)用不同的模型來描述。系統(tǒng)的分類方法很多，重要的分類方

法是按系統(tǒng)的狀態(tài)是否隨時(shí)間變化來分：一個行為與時(shí)間有關(guān)的系統(tǒng)稱為動態(tài)系統(tǒng)

待研系統(tǒng)：靜態(tài)系統(tǒng)：系統(tǒng)的行為與時(shí)間無關(guān)；用靜態(tài)模型來描述，一般為數(shù)學(xué)方程、邏輯表達(dá)式等.

如，電路的布爾表達(dá)式；電路中電壓與電流的關(guān)系；系統(tǒng)的穩(wěn)態(tài)解公式等

動態(tài)系統(tǒng)：連續(xù)系統(tǒng)：系統(tǒng)的狀態(tài)隨時(shí)間連續(xù)變化；常用微分方程來描述；方程對所有

時(shí)間點(diǎn)有效。如，衛(wèi)星運(yùn)行軌道，作物生長量等

確定性系統(tǒng)：系統(tǒng)的輸出完全由其輸入來描述，即系統(tǒng)輸入與

輸出按某種規(guī)則—對應(yīng)

集中參數(shù)模型：用常微分方程來描述，即方程中的

導(dǎo)數(shù)不是偏導(dǎo)數(shù)

分布參數(shù)模型：用偏微分方程來描述，但一般用

集中參數(shù)模型近擬表示

隨機(jī)系統(tǒng)：系統(tǒng)的輸出是隨機(jī)的，有規(guī)律的存在一族隨機(jī)變量，

且隨機(jī)變量序列與時(shí)間有關(guān)(隨機(jī)過程)。

在確定性系統(tǒng)的模擬中使用隨機(jī)變量的研究方法

稱為蒙特卡羅方法

離散系統(tǒng)：系統(tǒng)狀態(tài)的變化只在離散的時(shí)間發(fā)生；動態(tài)方程只在離散點(diǎn)上有效；

一般為隨機(jī)系統(tǒng)。如，庫存問題；企業(yè)的管理系統(tǒng)等

離散時(shí)間系統(tǒng)：時(shí)間步長固定；常用差分方程來描述

離散事件系統(tǒng)：用事件來表示系統(tǒng)在時(shí)間間隔內(nèi)的變化；常用

概率模型來描述

1.3建立數(shù)學(xué)模型

對一個系統(tǒng)，確定其類型后有助于選擇合適的方法建模，建模的一般步驟可分兩個大的階段：

(1)實(shí)質(zhì)內(nèi)容模型階段：首先對模擬對象進(jìn)行調(diào)查(了解系統(tǒng)，搜索模擬所需的信息)、實(shí)驗(yàn)(參數(shù)估計(jì)

等統(tǒng)計(jì)推斷方法，確定參數(shù)及參數(shù)的敏感性)、分析(將信息分類、量化，確定描述系統(tǒng)的規(guī)則)，盡可能全

面地掌握系統(tǒng)的基本特性、運(yùn)動規(guī)律以及中間狀態(tài)(最好是有對系統(tǒng)有深入了解的專業(yè)人員參與)，通過分

析和邏輯推理，揭示系統(tǒng)內(nèi)的規(guī)律。

(2)形式數(shù)量階段：在調(diào)查、實(shí)驗(yàn)、分析的基礎(chǔ)上，進(jìn)一步揭示系統(tǒng)內(nèi)部的數(shù)量關(guān)系，并對其進(jìn)行數(shù)

學(xué)處理，即對系統(tǒng)用數(shù)學(xué)形式來描述：用變量描述系統(tǒng)狀態(tài)，用各種數(shù)學(xué)方程定量表示各變量之間的相互

聯(lián)系，用遞歸方程描述系統(tǒng)狀態(tài)的發(fā)展趨勢。

多數(shù)情況下，建立一個好的(能真實(shí)的描述系統(tǒng)，有代表性，能準(zhǔn)確的模擬系統(tǒng)的數(shù)量信息，與實(shí)際系

統(tǒng)較吻合)數(shù)學(xué)模型不易，優(yōu)其是復(fù)雜多變的系統(tǒng)或系統(tǒng)本身的特性尚不完全清楚的情況(如對農(nóng)作物開發(fā)

新品種)。所以在數(shù)學(xué)模型建立后，必須進(jìn)行模擬驗(yàn)證，如與真實(shí)系統(tǒng)相差較大，則要重新建模或修改模型。

所以一個好的數(shù)學(xué)模型必須經(jīng)過多次模擬，不斷修改、完善，才能得到。

一個數(shù)學(xué)模型是描述系統(tǒng)行為的一個算法或一系列方程，工程中對系統(tǒng)建模應(yīng)先建立系統(tǒng)的需求規(guī)格

說明，在制定模擬規(guī)劃之前予以充分討論，過程中需考慮以下因素：

(1)模型中需考慮哪些有效的因素：一個系統(tǒng)可能有不同的行為(如一個作物生產(chǎn)系統(tǒng)中有作物的長

勢、產(chǎn)量、質(zhì)量、病蟲害等)，但不一定所有這些因素都要建立在模型中。模型中需要考慮的因素是：真正

能簡化系統(tǒng)的建模；對系統(tǒng)易于建模、測試和維護(hù)；使用較少的計(jì)算資源；對研究的系統(tǒng)有直接作用的。

(2)建模細(xì)化到什么程度：根據(jù)系統(tǒng)的需求確定細(xì)化的程度，是只須建立一個簡單的模型，還是要對

系統(tǒng)行為精確描述。在模型的準(zhǔn)確性和花費(fèi)之間求得平衡。

(3)與系統(tǒng)有相互作用的哪些外郵不學(xué)考慮在建模中：如在作物生長模型中，氣候、水、肥等。

(4)在建模中采取什么技術(shù)：首先，是基于物理的方程，還是基于測試數(shù)據(jù)。如果基于物理的方程，

是用微分方程，還是差分方程，考慮不考慮隨機(jī)因素等。這個問題常由專業(yè)人員根據(jù)系統(tǒng)的專業(yè)知識來決

定；如果基于試驗(yàn)數(shù)據(jù)，則建立經(jīng)驗(yàn)方程。

(5)建模時(shí)必須獲得什么樣的數(shù)據(jù)：如林木生長模型中，需要胸徑、樹高、材積等數(shù)據(jù)。另外還需考

慮這些數(shù)據(jù)的方便的輸出形式，以及模型需要的計(jì)算資源，如模型需要占用的內(nèi)存、磁盤空間、消耗的

CPU時(shí)間等。

(6)建模和測試模型需多長時(shí)間，多少人力、物力、財(cái)力：隨著模型復(fù)雜性的增加，成本也會增加。

一般可從簡單開始，隨著應(yīng)用逐步完善。

(7)如何驗(yàn)證和確認(rèn)模型：必須確認(rèn)建立的模型被正確實(shí)現(xiàn)，模型所描述的行為與真實(shí)系統(tǒng)匹配到可

接收的程度，才能有價(jià)值。那么以什么標(biāo)準(zhǔn)來衡量。現(xiàn)在已有驗(yàn)證、確認(rèn)與簽定(VV&A)技術(shù)來證明和驗(yàn)

證模擬的精確性。

以上問題應(yīng)列在系統(tǒng)需求說明書中，并應(yīng)用于最高層次：在實(shí)際的模擬問題中，可能將系統(tǒng)分解為子

系統(tǒng)、組件，在對各組件、子系統(tǒng)建模時(shí)仍須遵從以上規(guī)則。

模擬步驟：

(1)明確系統(tǒng)

(2)建立模型

(3)模型變換

(4)軟件設(shè)計(jì)開發(fā)

(5)測試檢驗(yàn)

(6)評估、對結(jié)果的評價(jià)和分析

一個模擬項(xiàng)目中各項(xiàng)工作的過程應(yīng)該是一

個迭代過程，如圖所示。下面通過實(shí)例來

說明模擬過程。

1.4應(yīng)用舉例

1.4.1兩物體追逐問題。設(shè)有一架殲擊機(jī)追蹤一架敵方轟炸機(jī)，假設(shè)兩機(jī)相距10公里以內(nèi)可實(shí)施攻擊,

且須在12分鐘內(nèi)完成追擊任務(wù)，否則認(rèn)為追擊失敗。設(shè)兩機(jī)初始位置如圖1-2所示。問題：對轟炸機(jī)的任

一條特定航線，模擬殲擊機(jī)的追擊過程。

圖1-2兩機(jī)追蹤模擬

分析：在這個系統(tǒng)中，是對指定的轟炸機(jī)的??條航線而言，模擬殲擊機(jī)的追擊情況：殲擊機(jī)按什么航

線飛行，何時(shí)完成追擊任務(wù)。能否殲滅敵機(jī)由在規(guī)定的時(shí)間內(nèi)兩機(jī)隨時(shí)間變動的距離而定，所以實(shí)施功擊

的時(shí)間、距離是要輸出的數(shù)據(jù)。在時(shí)刻t兩機(jī)距離由殲擊機(jī)在t時(shí)的位置決定，而t時(shí)的位置依賴于其速度

和航向。

為使模型簡單，作如F簡化：

(1)設(shè)兩機(jī)在同一平面飛行。于是三維問題轉(zhuǎn)化為二維問題。

(2)設(shè)殲擊機(jī)的速度(必須考慮的因素)VF是常數(shù)(20km/mim)。變速須用微分方程描述，而常速即可用

?般方程表達(dá)，求解更簡單。

(3)設(shè)殲擊機(jī)的航向(必須考慮的因素)在At(設(shè)1分鐘)改變一次，而在1分鐘以內(nèi)操作不變。這樣曲

線運(yùn)動即變?yōu)檎劬€運(yùn)動。

(4)轟炸機(jī)航向(航線)可任意，比如現(xiàn)給定一條航線，如表1-1所示。

表1-1轟炸機(jī)航線

t012345678901112

XB(t)809099108116125133141151160169179180

YB(050484541353227212225293033

(5)殲擊機(jī)初始位置：XF(O)=O,YF(O)=O(初始條件)

建立數(shù)學(xué)模型：

設(shè)在t時(shí)刻兩機(jī)位置為(XF()YF⑴)、%⑴,YB(。),兩機(jī)連線與水平線夾角為0,則1分鐘后殲擊機(jī)

的位置為：

XF(t+l)=XF(t)+VFCosO(1.1)

Yh(t+l)=Yl.(t)+VbSin0(1.2)

由(XF(t),YF(t))計(jì)算(XF(t+l),YF(t+l))時(shí)涉及角0,從圖中看出:

SinO=[YB(t)-Y[.(t)]/D(t)(13)

CosO=[XB(t)-XF(t)J/D(t)(1.4)

而

2

D(t)=SQR{[YB(t)-YF(t)?+[XB(t)-XF(t)]}(1.5)

式(1.5)-(1.1)即是這個問題的數(shù)學(xué)模型。先算出兩機(jī)之距離，不斷判斷是否在12分鐘內(nèi)到達(dá)可攻擊的

距離之內(nèi)。程序流程圖如圖1-3所示。

此例是連續(xù)系統(tǒng)，確定性模型，用一組方程來描述。下例是一個離散事件系統(tǒng)，是隨機(jī)系統(tǒng)，用概率

模型來描述。

開始

圖1-3追擊模擬流程圖

VB程序見實(shí)例。

1.4.2庫存問題。商業(yè)部門為了合理利用有限的流動資金，每項(xiàng)商品都要在庫存與銷售之間求得平衡:

庫存量太大會增加管理費(fèi)用、枳壓資金；庫存量太小又可能造成缺貨，也會造成銷售損失、信譽(yù)損失。所

以，當(dāng)庫存量不滿足某一時(shí)段的顧客需求時(shí)，就要到廠家訂貨。這就需要采取一種策略：當(dāng)庫存量（比如布

匹）降到P匹布時(shí)（稱為車蛻鄉(xiāng)）則向廠家訂貨，定貨量為Q匹（稱為本袋拿）。假設(shè)現(xiàn)在有5種策略（方案）。

從中選擇一種使花費(fèi)最少。

表1-2庫存策略

策略PQ

1125150

2125250

3150250

4175250

5175300

已知條件：

（1）從發(fā)出訂貨單到收到貨物需3天，即第i天訂貨，第i+3天收到。

（2）每匹布的保管費(fèi)為0.75元，缺貨損失為1.8元/匹，訂貨費(fèi)（包括手續(xù)費(fèi)、采購差旅費(fèi)及其他費(fèi)用）

為750元。

（3）需求量為一個0-99之間的均勻分布隨機(jī)數(shù)。

（4）原始庫存量為115匹，并設(shè)第?天沒發(fā)出訂貨。

分析：庫存問題是商業(yè)上的一個重:要問題在數(shù)學(xué)上有專門的模型研究，存儲論也是運(yùn)籌學(xué)的一個重要

分支。庫存問題用計(jì)算機(jī)模擬最合適，若通過銷售來找最優(yōu)方案，勢必造成經(jīng)濟(jì)損失。這是典型的離散事

件系統(tǒng)，用概率模型來描述。這里已給定概率模型:

X-U[0,99]即密度函數(shù)為：p(x)=1/990WxW99

r0x<0

分布函數(shù)為：F(x)=?x/990WxW99

.1x>99

以F框圖對每種策略模擬180天,選出效益最好的方案，如圖1-5所示。

圖1-5庫存問題模擬程序框圖

C程序見實(shí)例。

以上兩個例子一個連續(xù)系統(tǒng)，狀態(tài)隨時(shí)間連續(xù)變化，用方程來描述；一個離散系統(tǒng)，到貨和銷售都按

一些離散的步驟進(jìn)行，存在隨機(jī)性，只能用概率模型來描述。但解決問題的過程、分析方法類似：先分析

系統(tǒng)，使對系統(tǒng)有充分的了解；建立數(shù)學(xué)模型，設(shè)定一些初始條件(如t=0時(shí)的狀態(tài)，開始時(shí)的庫存量)；

系統(tǒng)狀態(tài)的變化對應(yīng)一組方程或一組規(guī)則，隨著時(shí)是的變化，系統(tǒng)狀態(tài)改變，當(dāng)一個周期結(jié)束時(shí)，收集模

擬過程的統(tǒng)計(jì)數(shù)據(jù)(即問題的解)。如果程序設(shè)計(jì)的好，就會使模擬非常逼真，就象一個真正的系統(tǒng)在演示

一樣。

從理論上講，任何問題都可用計(jì)算機(jī)模擬。計(jì)算機(jī)模擬具有經(jīng)濟(jì)、安全可靠、周期短等優(yōu)點(diǎn)。對任何

問題，只要建立起數(shù)學(xué)模型，改變參數(shù)值及變量值，就可模擬各種情況卜的系統(tǒng)運(yùn)行情況，從模擬輸出的

結(jié)果，可分析系統(tǒng)內(nèi)各因素的權(quán)重及其制約關(guān)系，幫助決策者作出合理的決策，克服盲目性，使系統(tǒng)在實(shí)

際運(yùn)行過程中取得最好的效益。

1.4.3傳染病傳播問題。假設(shè)某一地區(qū)有一種傳染病正在流行，那么政府、醫(yī)務(wù)部門就要采取措施來

控制這種疾病的傳播，要使采取的措施能夠有效的控制傳攀速度，就必須先搞清被傳染的人數(shù)跟哪些因素

有關(guān)，被傳染的人數(shù)是一個什么樣的發(fā)展趨勢，從而來有效的預(yù)測和控制，下面建立描述被傳染人數(shù)的數(shù)

學(xué)模型。

傳染病的傳播涉及因素很多，不可能通過一兩次簡單的假設(shè)就能建好完善的數(shù)學(xué)模型。這里的作法是:

先作出最簡單的假設(shè)，看會得到什么樣的結(jié)果，然后對不合理的地方再行修改，逐步得到較滿意的模型。

先討論一個粗略的模型。

模型I：假設(shè)：

(1)每個病人在單位時(shí)間內(nèi)傳染其他人的人數(shù)為一個常數(shù)k0?

(2)一人得病后，經(jīng)久不愈，該人在傳染期內(nèi)不會死亡。

記時(shí)刻t的得病人數(shù)為y(t),開始模擬時(shí)有y0個傳染病人，則在At時(shí)間內(nèi)增加的病人人數(shù)為

y(t+At)-y(t)=koy(t)At

由導(dǎo)數(shù)定義得(在假設(shè)(1)、(2)下的數(shù)學(xué)模型)：

rdy/dt=koy(t)(1.6)

?y(0)=yo

初值問題(1.6)的解為：y(t)=yoek。'

(分析：)這個結(jié)果表明，得病人數(shù)將按指數(shù)形式無限增加，當(dāng)t-8時(shí)，y(t)-8,顯然與實(shí)際不符，

說明上面的假設(shè)條件不合理。事實(shí)上，一個地區(qū)的總?cè)藬?shù)大致可視為常數(shù)(不考慮疾病傳播期間出生的、遷

出的、死亡的)，所以一個病人在單位時(shí)間內(nèi)能傳播的人數(shù)k0是在改變的：在初期，k0較大，隨著病人的

增多，健康人數(shù)的減少，被傳染的機(jī)會也將減少，所以ko將變小。

對上述模型進(jìn)行改進(jìn)：記t時(shí)刻的健康人數(shù)人S(t),當(dāng)總?cè)藬?shù)不變時(shí)，k0隨S(t)減少而減少。

模型H：假設(shè)：

(1)總?cè)藬?shù)為常數(shù)n,t時(shí)刻的健康人數(shù)為S⑴，得病人數(shù)為Y(t),則

Y(t)+S(t)=n(1.7)

(2)單位時(shí)間內(nèi)一個病人傳染的人數(shù)與當(dāng)時(shí)的健康人數(shù)成正比，比例系數(shù)為k(醫(yī)學(xué)上稱為,學(xué)學(xué)層)

(3)同模型1的假設(shè)(2)。

由假設(shè)(2),(1.6)式中的Q應(yīng)為kS⑴,即:

rdy(t)/dt=kS(t)Y(t)(1.8)

ty(0)=y0

將(1.7)式代入(1.8)式，得(上述假設(shè)下的數(shù)學(xué)模型)：

'dy(t)/dt=kY(t)(n-Y(t))(1.9)

.y(O)=yo

初值問題(1.9)的解為：

-kn,

Y(t)=n/[l+(n/y0-l)e](1.10)

(分析：)由(1.10)式得:

2-kn,-knl2

dy/dt=[kn(n/y0-l)e]/[1+(n/y0-I)e](1.11)

(1.10)式是被傳染人數(shù)隨時(shí)間變化的關(guān)系；(1/1)式是被傳染病人的變化率與時(shí)間的關(guān)系，如圖1-6和

圖1-7所示。

圖1-6病人人數(shù)變化曲線圖1-7病人變化曲線

這個模型可預(yù)報(bào)疾病高峰到來的時(shí)間：令d2y/dt2=O,得極大值點(diǎn)：

tj=ln(n/yo-l)/kn(1.12)

由(12)式可知，當(dāng)傳染強(qiáng)度k或人數(shù)n較大時(shí)，L變化較小，表明傳染高峰到來的快，這與實(shí)際情況相吻

合。但由(L10)式知當(dāng)t-8時(shí)，丫⑴一必這意味著最終人人都被傳染，這與實(shí)際不符，原因在于假設(shè)(3)

不合理。

再改進(jìn)：可將人員分為三類：第一類為傳染者(y)；第二類為易受傳染者(s),即這一類是非傳染者，但

能得病而成為傳染者；第三類為除前兩者之外的人(r),包括患病死去的、痊愈后具有長期免疫力的、被隔

離的。

用y(t)、s(t)、r(t)分別表示這三類人數(shù)。

模型HI：假設(shè)：

(1)總?cè)藬?shù)為常數(shù)n,貝IJ：

y(t)+s(t)+r(t)=n(1.13)

(2)同模型H的假設(shè)(2)

(3)單位時(shí)間內(nèi)病愈免疫的人數(shù)與當(dāng)時(shí)病人的人數(shù)成正比，比例系數(shù)為m(恢復(fù)系數(shù))

由假設(shè)⑶,有

dr/dt=my(t)(1.14)

由于引入了r(t),則模型II的方程(1.8)應(yīng)改為：

dy/dt=kS(t)y(t)-dr/dt(1.15)

(1.15)式表示單位時(shí)間內(nèi)病人人數(shù)的增加應(yīng)等于被傳染的人數(shù)減去病愈的人數(shù)。從(1.13)-(1.15)中消去dy,

并設(shè)S(0)=So,r(0)=r0得

dS/dt=-kS(t)y(t)、

dr/dt=my(t)

y⑴+s(t)+r(t)=n>(1.16)

S(0)=So

r(0)=r0.

模型(1.16)較好地描述了傳染病傳播問題。

通過以上實(shí)例，對數(shù)學(xué)模型的建立就有了一個基本的思路，對計(jì)算機(jī)模擬技術(shù)、模擬的過程、問題的

方法步驟也有了一個概括的了解。后面兩章下分別對連續(xù)系統(tǒng)和離散系統(tǒng)討論基本的模擬方法。

習(xí)題：

(1)對例1.4.1中的飛機(jī)追擊問題采用極坐標(biāo)(r,?),相應(yīng)的方程為：

dr/dt=VBCos9-VF

r,d0/dt=-VBSin0

其中r為兩機(jī)距離，0為兩機(jī)連線的角度。編程模擬兩機(jī)追擊情況。

（2）球擺：如圖1-8所示，設(shè)繩長為L,夾角為。，球質(zhì)量為m,初始速///（///

度。（0）=0,初始偏角0°,確定擺動周期。

提示：影響擺動周期的因素：|9\

①擺球重力；|\

②擺球質(zhì)量：、

③擺球尺寸相對于繩很小，故可視為質(zhì)點(diǎn)建模；

④繩子質(zhì)量忽略；圖1-8球擺

⑤磨擦力忽略；

⑥空氣阻力不考慮。

建模：影響擺球運(yùn)動的重力（使擺球回到中心位置的力）F的方向與繩子垂直；對繩子產(chǎn)生拉緊力的重

力與繩子平行，它不影響擺球運(yùn)動，忽略。F=-mgSinO,將牛頓定律F=ma代入上式得：a=-gSin0,

其中a是擺球的切線加速度，由于L=0”,故得運(yùn)動方程：

0"=-L/g?Sin9

e(0)=fio*0'(0)=0

第二章連續(xù)系統(tǒng)的計(jì)算機(jī)模擬技術(shù)

連續(xù)系統(tǒng)的狀態(tài)隨時(shí)間連續(xù)的動狀變化，這種變化依賴相關(guān)的因素，且遵從一定的規(guī)律，所以一般來

說可用數(shù)學(xué)方程描述(代數(shù)方程、微分方程，此外還有傳遞函數(shù)、結(jié)構(gòu)圖及狀態(tài)方程等)。問題：對數(shù)學(xué)模

型怎樣求解,求解過程也就是模擬的過程，即程序的算法。對前面的飛機(jī)追擊的問題，其模型是一組代數(shù)

方程，而且是亞4的，所以，作法是：假設(shè)1分鐘改變一次航向，每隔1分鐘計(jì)算一次系統(tǒng)的狀態(tài)，而且

山模型知道了殲擊機(jī)前一分鐘的狀態(tài)就可以算出后一分鐘的位置。這樣，通過迭代，就可求出問題的解。

對一般的模擬問題，模型用微分方程表達(dá)(方程中除含有自變量外，還有導(dǎo)數(shù)或偏導(dǎo)數(shù)一分布參數(shù)型可

通過一些變換轉(zhuǎn)換為集中參數(shù)型)，而對一個常微分方程(只有線性的或幾種特殊的能求出通解及特解，即

解析解)?般來說求解析解是不可能的，(即使能求出解析解)在計(jì)算機(jī)上求解要用數(shù)值積分法：將連續(xù)的系

統(tǒng)離散化，把微分方程轉(zhuǎn)化為差分方程，通過迭代運(yùn)算，求出問題的數(shù)值解。通過正確的控制步長、選擇

合理的算法即能達(dá)到要求的精度，這就是連續(xù)系統(tǒng)模擬的主要技術(shù)。

數(shù)值積分法種類很多，本章介紹幾種簡單常用的方法。

2.1歐拉(Euler)法

設(shè)連續(xù)系統(tǒng)的數(shù)學(xué)模型為

rdy/dl=f(t,y)(2.1)

Iy(0)=y0

為了模擬系統(tǒng)狀態(tài)y隨時(shí)間t的變化，需求解微分方

程(2.1)的數(shù)值解(不是解析解)，為此，把模擬周期分為

若干小區(qū)間，比如分為N個相等的小區(qū)間，如圖2-1

所示。只須在每個離散的點(diǎn)上計(jì)算系統(tǒng)的狀態(tài)。

在區(qū)間9,t向)上求積分，得

y(tn+i)-y(tn)=f(t,y)dt圖2-1歐拉法的兒何意義

枳分的幾何意義是小曲邊梯形的血枳。

如果小區(qū)間取的足夠小，則在區(qū)間(t.,tm)之間的f(t,y)可近似的看成常數(shù)f(tn,yn),這樣可用小矩形的

面積代替小曲邊梯形的面積，于是得在t用時(shí)的積分值為

y(tn+i)yn+f(tn,yn),h=yn+,(2.2)

其中h=T/N,將(2.2)式寫為差分方程形式，得

yn+l=yn+f(tn,yn)?hn=0,1,2,—(2.3)

這就是歐拉公式：它是一個遞推的差分方程，任何一個新的數(shù)值解yn+l均基于前一個數(shù)值解yn以及導(dǎo)

數(shù)值f(tn,yn)求得,只要給出初始條件yo及步長h,就可算出力，由力可算出y2,如此迭代計(jì)算y3,y4,…,

直到滿足所需計(jì)算的范圍。

歐拉法也叫折線法，特點(diǎn)是方法簡單，計(jì)算量小，但計(jì)算精度也底。

擴(kuò)展I：改進(jìn)的歐拉法(梯形法)

為了提高計(jì)算的精度，可對歐拉法進(jìn)行改進(jìn)：從幾何意義上看到，用小梯形的面積來代替小矩形的面

積，必能提高精度。這樣(2.3)式即可寫成：

yn+l=yn+h/2?[f(tn,yn)+f(tn+byn+i)]

=yn+h/2-[fn+fn+IJ(2.4)

但(2.4)式是隱式公式，即公式右端含有yn+1,而這是未知的待求量，故梯形法不能自行啟動運(yùn)算，要借助

于其它算法：如用歐拉法算出y1+i作為梯形法中yC用的預(yù)估值，再進(jìn)行計(jì)算，這樣，公式可寫為：

預(yù)估：p

yn+i=yn+f(tn,y?),h

校正：cp

yn+i=y?+h/2?[f(tn,yn)+f(tn+l,yn+l)])(2.5)

p

=yn+h/2-[fn+fn+1]

(2.5)式也稱為預(yù)估校正法，計(jì)算量稍大(需要附加的計(jì)算)，但精度要比歐拉法高，穩(wěn)定性好。

擴(kuò)展H：多個輸入的多階系統(tǒng)

在(2.1)式所表示的數(shù)學(xué)模型中，y是系統(tǒng)的狀態(tài)，稱之為不春變量，它隨時(shí)間t而連續(xù)變化。(2.1)式是

最簡單的情況，只有一個狀態(tài)變量，而它直接依賴時(shí)間而變化，在實(shí)際模擬中，狀態(tài)變量可能不只一個，

而影響系統(tǒng)狀態(tài)的因素更多，這些因素是隨時(shí)間變化的同，而系統(tǒng)狀態(tài)的改變也正是由于這些因素隨時(shí)間

變化而變化。如在林木的生長系統(tǒng)中，要考察的狀態(tài)可能有樹高、胸徑、材積等；在作物的生長系統(tǒng)中，

考察的狀態(tài)可能有株高、葉子的片數(shù)、葉子的寬度等，而影響系統(tǒng)狀態(tài)的因素如水、肥、氣候、土質(zhì)、病

蟲害等，而這些因素都隨時(shí)間而變化，如水份會隨時(shí)間而被蒸發(fā)，肥會隨時(shí)間被作物吸收，氣候等自然環(huán)

境更是隨時(shí)間變化無常。所以與(2.1)式相應(yīng)的數(shù)學(xué)模型應(yīng)為：

Jdyk/dt=f(x,(t),x2(t),..?xm(t),yby2,...,yq)k=1,2,...,q(2.6)

Iyk(O)=yk.o

其中yk為系統(tǒng)狀態(tài)變量，Xi為系統(tǒng)輸入信號，這樣的系統(tǒng)稱為縣直皿加揄羽統(tǒng)。(2.6)式是q個微

分方程的方程組，q個初始條件。

可以把(2.6)式表示成向量的形式，就和(2.1)式在形式上一致了：記Y為q個狀態(tài)變量和向量，X為m

個輸入信號向量，即：

Y=(y1(t),y2(t),…，yq⑴)，X=(x)(t),x2(t)，…，xm(t))

則(2.6)式寫為：

rdY/dt=f(X(t),Y)(2.6)，

LY(0)=Yo

對于(2.6)式所表達(dá)的數(shù)學(xué)模型，與(2.3)式相對應(yīng)的歐拉公式為：

fyk.n+1=丫k.n+f(X](tn),X2(tn),…，Xm(tn),y1.n,丫2,n,…，Yq,n)*h(2.7)

Ik=12…,q

或?qū)懗上蛄康男问剑?/p>

Yn+l=Yn+f(Xn,Yn)-h(2.7),

例2.1設(shè)兩種物質(zhì)A和B合到一起產(chǎn)生化學(xué)反應(yīng)，生成新物質(zhì)C,假設(shè)1克A和1克B結(jié)合能產(chǎn)生

2克C,形成C的速率與A和B的數(shù)量乘積成正比。同樣C也可分解為A和B,C分解的速率正比于C

的數(shù)量。問題：在給定A和B的數(shù)量后，模擬有多少C物質(zhì)產(chǎn)生出來，以及達(dá)到穩(wěn)定的時(shí)間。

解：在任何時(shí)刻，設(shè)a,b,c分別是A,B,C的數(shù)量，則它們增加和減少的速度服從下列微分方程:

'da/dt=k2C-k|a,b

■db/dt=k2C-kIa?b(2.8)

.de/dt=2k|a,b-2k2C

其中k1和k2是比例常數(shù)(實(shí)際問題中心和k2會隨溫度、壓力等發(fā)生變化，但在模擬過程中為簡化模型，可

視為常數(shù))，是模型中的參數(shù)，-k|a?b中的負(fù)號是因?yàn)锳和B是減少的過程。假設(shè)模擬從1=啾一般取0)

開始，使t以時(shí)間間隔增力口(由步長at和模擬周期可定出所分時(shí)間段數(shù))，則相應(yīng)的歐拉公式為(2.9)式。

給出常數(shù)％和k2值以及A、B、C的初始數(shù)量，即得迭代公式：

'an+i=an+(k2Cn-klan?bn)?At

bn+l=bn+(卜2Cn-k|an*t>n)*

?cn+1=cn+(2k|a??bn-2k2Cn)?At(2.9)

a(0)=ao,b(0)=bo,c(0)=0

、k|=ki,o,k2=k2,o

由式，從開始，由。、、可算出、、又可算出

(2.9)t=0ab0Coa?b]c”a2'b2>c2,(a2=a(2At)...),

以為間隔，進(jìn)行N=T/4次計(jì)算，就可算出周期T的系統(tǒng)狀態(tài)，從而得出模擬結(jié)果。

根據(jù)數(shù)學(xué)模型（2.9）就可設(shè)計(jì)編寫模擬程序，可選用任何編程語言，在一些流行的語言中，對常用的模

擬算法（比如后面要介紹的龍格?庫塔（Runge?kutta）方法）都有相應(yīng)的用于模擬計(jì)算的軟件包。

publicclassRate{

staticdoublekl,k2,h,t;

staticdoubleA[J=newdoublel53J;

staticdoubleB[]=newdouble[53];

staticdoubleC[J=newdoublel53J;

staticinti;

publicstaticvoidmain(Stringargs[J){

A[ll=l00.0;B[l]=50.0;C[l]=0.0;

t=O;h=O.l;

kl=0.008;k2=0.002;

for(i=l;i<53;i++){

System.out.println(i+""+t+",

n+A[i]+';"+B[i]+”，n+C[i]);

strut(i);

staticvoidstrut(inti){

A[i+l]=A[i]+(k2*C[i]-kl*A[i]*B[i])*h;

B[i+l]=B[i]+(k2*C[i]-kl*A[i]*B[i])*h;

C[i+l]=C[i]+2.0*(kl*A[i]*B[i]-k2*C[i])*h;

t=t+h;

)

圖2-2例2-1模擬流程圖

)

歐拉法是最簡單的數(shù)值積分法，在介紹更好的數(shù)值積分法之前，先介紹幾個關(guān)于數(shù)值積分的基本概念。

1、單步法與多步法

數(shù)值積分法是用遞推公式求解，如果僅由前一時(shí)刻的數(shù)值yn就能算出后一時(shí)刻的數(shù)值yn+l，則稱為單

步法，反之，如果求y用時(shí)需用到y(tǒng)n，yn.1，yn々，…等多個值，則稱為多井渚。如，歐拉法是單步法，而

改進(jìn)和歐拉法是多步法。

單步法由遞推公式自身就能啟動運(yùn)算（由初值即可算出yi，由力可算出y2,…），所以它是能自啟動

的算法；多步法在開始的時(shí)候要先用其它的方法計(jì)算該時(shí)刻前面的函數(shù)值，所以不能自啟動運(yùn)算。

一般來說，由于多步法更能充分利用多個時(shí)刻的信息，所以模擬速度快，精度高，但計(jì)算量要大一些,

后面還要介紹多步法的算法。

2、顯式與隱式

在遞推公式中，計(jì)算y用時(shí)所用到的數(shù)據(jù)均已算出的計(jì)算公式稱為顯式公式，相反，在算式中隱含未

知量稱為用耳公耳。如歐拉法中遞推公式是顯式的，而梯形法中遞推公式是隱式的。

在隱式公式中，必須先用顯式公式估計(jì)一個值，再用隱式公式迭代，即預(yù)估-校正法。隱式與顯式相比,

有明顯的高精度和穩(wěn)定性。

3、誤差

在數(shù)值解法的過程中，每一步計(jì)算都會產(chǎn)生誤差，誤差的來源有兩個方面：一是計(jì)算誤差（即計(jì)算機(jī)計(jì)

算本身的誤差)，一個是公式誤差(用差分方程代替微分方程)。分別稱之為舍入誤差和截?cái)嗾`差。

舍入誤差：計(jì)算機(jī)的字長是有限的，數(shù)字不能表示的完全精確，實(shí)際上計(jì)算的結(jié)果是用有限精度的有

理數(shù)(如計(jì)算機(jī)中使用的浮點(diǎn)數(shù))來近似無限精度的實(shí)數(shù)，所以在對動態(tài)系統(tǒng)模擬的過程中，舍入誤差是不

可避免的。通常舍入誤差的大小與積分步長成反比，步長h越小，計(jì)算次數(shù)多則舍入誤差大，但不能隨意

加大步長，否則將產(chǎn)生大的截?cái)嗾`差甚至影響系統(tǒng)穩(wěn)定性。在給定運(yùn)算序列的條件下，唯一可減少舍入誤

差的方法是增加數(shù)字表示的精度，如將單精度浮點(diǎn)數(shù)(有7位數(shù)的精度)表示改為雙精度數(shù)(有15位數(shù)的精

度)。實(shí)際的動態(tài)系統(tǒng)模擬時(shí)，多采用雙精度數(shù)據(jù)類型，盡管這樣變量占用的存儲空間大，處理時(shí)間稍長，

但對提高模擬的精度來講是值得的。

截?cái)嗾`差：是用差分方程代替微分方程產(chǎn)生的誤差，即用數(shù)值解代替微分方程的精確解產(chǎn)生的誤差，

所以是公式本身的誤差，一般用臺勞級數(shù)來分析積分公式的精度：假設(shè)在tn時(shí)積分(精確值)已經(jīng)算出，則

用臺勞級數(shù)可求得t用時(shí)的精確解：

2,5(r)r+l

y(t11+1)=y(tn+h)=y(tn)+y'(tn)+h/2!?y(tn)+…+h「/r!?y(tn)+0(h)(2.10)

如果一個數(shù)值解法是用前r+1項(xiàng)來近似的計(jì)算y(t用)，則后面的各項(xiàng)(記為)O(h向)是在這一步計(jì)算中引進(jìn)

的附加誤差，稱為這個算法的(局部)截?cái)嗾`差。誤差O(h"i)與h向同階(h-0),即計(jì)算的精度保特了r階,

此時(shí)稱這個方法是r階的精度。一個數(shù)值方程的階數(shù)可視為衡量這一方法的精確度的重要標(biāo)志，不同的數(shù)

值解法有不同的精度。

歐拉法只取(2.10)中的前兩項(xiàng)作近似計(jì)算，所以是一階精度的算法；梯形法相當(dāng)于取(2.10)式中的前三

項(xiàng)，故是二階精度的方法。

4、穩(wěn)定性

如果系統(tǒng)是穩(wěn)定的(系統(tǒng)的狀態(tài)隨時(shí)間的推移逐步穩(wěn)定在某個水平上)，則在模擬的迭代過程中，數(shù)值

積分的解也應(yīng)該是穩(wěn)定的，但由于初始數(shù)據(jù)的誤差及在迭代運(yùn)算中的舍入誤差會對后面的計(jì)算結(jié)果產(chǎn)生影

響。當(dāng)步長h選擇不合理時(shí)，可能使模擬結(jié)果不穩(wěn)定。對于歐拉法，可用下面的檢驗(yàn)方程(其中人為方程的

特征根)：

jy'=Ay(2.11)

'y(o)=yo

(注意其精確解為y=y。e”)來檢驗(yàn)步長對數(shù)值解穩(wěn)定性的影響：對(2.11)表示的方程，歐拉公式為：

x

(-yn+i=yn+ynh=(i+xh)yn(2.12)

Iy(0)=y0

所以，要使數(shù)值解穩(wěn)定，必須使：11+Ahl<1,解得|Ahl<2或h<2T(T=1/入是系統(tǒng)時(shí)間常數(shù))。所

以要保證歐拉方法計(jì)算的穩(wěn)定性，步長h必須小于系統(tǒng)時(shí)間常數(shù)的2倍。

在(2.11)中特征根人在一定數(shù)量級范圍變動，現(xiàn)令人=1,來作?個具體的模擬，此時(shí)h應(yīng)小于2,否

則將不穩(wěn)定。取h=0.01,h=0.05,h=1.0,h=1.9,h=2.0,h=2.1六個不同的步長，y0=1?模擬結(jié)果

為：

h=0.01：表現(xiàn)出較好精度

h=0.05：雖然近似解接近精確解，但存在誤差

h=1.0：解在一步后趨于零，并一直保持，雖然穩(wěn)定，但明顯精度不高

h=1.9：解振蕩，幅度值逐衰減，并趨于穩(wěn)定

h=2.0：解不衰減，等幅振蕩

h=2.1：解的振蕩幅度值遞增，表明系統(tǒng)不穩(wěn)定，數(shù)值解發(fā)散

數(shù)值解圖如圖2-3所示。

22

11\

00----

-1-1

-2-2

22

00

-1-1

-2-2

22

00

-1-1

-2-2

圖2-3不同步長的數(shù)值解

2.2龍格-庫塔(Runge-Kutta)法

R-K方法的基本思想是用臺勞展開式的前幾項(xiàng)來對微分方程求近似解。再以模型(2.1)為例：

ry'=f(t,y)(2.1)

[y(to)=yo

假設(shè)從to開始，以h增長，h|=to+h,yi=y(to+h),在t0附近展開成臺勞級數(shù)，保留1?項(xiàng)，則有：

yi=yo+f(to,yo)*h+h2/2(8f/5y?dy/dt+8f/8t)(2.13)

(此式括號中的導(dǎo)數(shù)是在(to,yo)處的導(dǎo)數(shù)值)為了求(2.13)的解，假設(shè)(2.13)的解可寫成如下的形式：

yi=y0+(b|kl+b2k2)?h-

其中k,=f(to,yo)\(2.14)

J

k2=f(to+C2h,yo+alklh)

(注意(2.13)式是關(guān)于函數(shù)y的導(dǎo)數(shù)的算式，而(2.14)式中心和k2都是y在某點(diǎn)處的導(dǎo)數(shù)，相差的只是

年拿季的系數(shù)不同，問題正是要定出這些系數(shù)，從而得到數(shù)值解表達(dá)式，為此)對k2式右端的函數(shù)在(to,y°)

處展開臺勞級數(shù)，保留h項(xiàng)，得：

k2—f(to,yo)+(C2?8f/5t+aikj?5f/6y?dy/dt)?h

把ki、k2代入(2.14)中，得

yi=yo+bjf(to,yo)h+b2h[f(t0,yo)+(C2§f/5t+ajkj5f/6ydy/dt)h](2.15)

將(2.15)式與(2.14)式比較，得關(guān)于系數(shù)的方程：

*bj+62=1

vb2c2=1/2(2.16)

〔b2al=1/2

方程組(2.16)中四個未知量，三個方程，故有無窮多解，求出一個解，比如令b|=b2,得一個解：

bj=62=1/2ai=C2=1

代入(2.14),得如下一個公式：

yi=yo+l/2(K|+K2)h-

其中jk|=f(to,y())”

*■k2=f(to+h,y()+kih))

這是從to計(jì)算t1時(shí)刻的公式，寫成一般的形式，得如下遞推公式：

yn+l=yn+h/2(K|+K2)'

其中rk!=f(tn,yn)?(2.17)

Ik2=f(tll+h,yn+k,h),

模型(2.1)的數(shù)值解公式(2.17)即稱為R-K公式，這種數(shù)值解法即稱為R-K方法。

在(2.13)式中，只保留了I?項(xiàng),故公式(2.17)的精度是2階的，公式(2.17)稱為二階R-K公式。

目前在實(shí)際模擬問題中，四階R-K公式用的最為普遍。在推導(dǎo)四階R-K公式時(shí)，在臺勞公式中保留

到。項(xiàng)，推導(dǎo)過程與前面類似。一個四階R-K公式可以是下面的形式：

Yn+l=yn+h/6(K|+2K2+2K3+K4)、

其中rk|=f(tn,yn)

,k2=f(tn+h/2,yn+klh/2)}(2.18)

k3=f(tn+h/2,yk2h/2)

Ik4=f(tn+h,yn+k3h)J

問題：關(guān)于R-K方法還有下面一些問題需了解。

(1)多種形式：方程組(2.16)有無窮多解，我們?nèi)×艘粋€解，得到了公式(2.17),也可以取其它的解，

所以每一階R-K公式都有多種形式。

二階R-K公式常用的除(2.17)式外，還有

?yn+i=yn+K2h

■ki=f(tn,y?)

'k2=f(tn+h/2,yn+k|h/2)

(對應(yīng)(2.16)的解b|=0,b2=l,C2=a1=1/2)。四階R-K公式常用的除(2.18)式外,還有

Lyn+l=yn+h/8(K|+3K2+3K3+K4)

k|=1(tn5yn)

<k2=f(tn+h/3,yn+k|h/3)

k3=f(tn+h2/3,yn+k!h/3+k2h)

k4=f(tn+h,yn+hk|-hk2+hk3)

(2)精度：也可推導(dǎo)3階、5階或更高階的R-K公式，但在一般工程中，四階公式就完全能達(dá)到要求

的精度，而且四階公式是最常用的，所以?般不使用更高階的公式。

另一個特殊情況：一階R-K公式，只含有h的1次項(xiàng)，即：y?+i=yn+hf(tn,yn)，這就是歐拉公式,

所以歐拉公式也可看面一階R-K公式，其精度最低。

R-K方法的精度取決于步長h及求解方法。一般來說，為達(dá)到同樣的精度，四階方法的步長可以比二

階方法的步長大10倍，而四階方法每步的計(jì)算量僅比二階方法大一倍，所以總的計(jì)算量仍比二階方法小。

R-K方法可使用較大的步長也是其一個特點(diǎn)。

(3)單步法：不論幾階的R-K公式都是單步法，在計(jì)算y