微分方程建模

微分方程建模當(dāng)我們描述實(shí)際對(duì)象的某些特性隨時(shí)間（空間）而演變的過(guò)程、分析它的變化規(guī)律、預(yù)測(cè)它的未來(lái)形態(tài)、研究它的控制手段時(shí)。通常要建立對(duì)象的動(dòng)態(tài)模型。

在許多實(shí)際問(wèn)題中，當(dāng)直接導(dǎo)出變量之間的函數(shù)關(guān)系較為困難，但導(dǎo)出包含未知函數(shù)的導(dǎo)數(shù)或微分的關(guān)系式較為容易時(shí)，可用建立微分方程模型的方法來(lái)研究該問(wèn)題。一、數(shù)學(xué)建模的基本思維過(guò)程

轉(zhuǎn)化實(shí)際問(wèn)題

1、對(duì)要討論的問(wèn)題所涉及的重要特征進(jìn)行合理的數(shù)學(xué)刻畫(huà)（轉(zhuǎn)化），即用數(shù)學(xué)語(yǔ)言對(duì)問(wèn)題涉及到的重要特征進(jìn)行表述.

2、尋求的實(shí)際問(wèn)題的文字?jǐn)⑹?，利用一些原則或定律，將其轉(zhuǎn)化為數(shù)學(xué)描述。解數(shù)學(xué)問(wèn)題用數(shù)學(xué)工具求解得到的數(shù)學(xué)問(wèn)題。二、微分方程模型涉及“改變”、“變化”、“增加”、“減少”、“衰變”、“邊際”、“速度”、“運(yùn)動(dòng)”、“追趕”、“逃跑”、、、等等詞語(yǔ)的確定性連續(xù)問(wèn)題。b、微分方程建模的基本手段微元法等a、微分方程建模的對(duì)象

1、尋找改變量一般說(shuō)來(lái)微分方程問(wèn)題都遵循這樣的文字等式變化率（微商）=單位增加量--單位減少量等式通常是利用已有的原則或定律。c、微分方程建模的基本規(guī)則2、對(duì)問(wèn)題中的特征進(jìn)行數(shù)學(xué)刻畫(huà)3、用微元法建立微分方程；4、確定微分方程的定解條件（初邊值條件）；5、求解或討論方程（數(shù)值解或定性理論）；

6、模型和結(jié)果的討論與分析。我們來(lái)建立如下的一些問(wèn)題的模型：

我們通過(guò)一些最簡(jiǎn)單的實(shí)例來(lái)說(shuō)明微分方程建模的一般方法。在連續(xù)變量問(wèn)題的研究中，微分方程是十分常用的數(shù)學(xué)工具之一。1、Malthus模型2、Logistic模型3、車(chē)間空氣的清潔

為了保持自然資料的合理開(kāi)發(fā)與利用，人類(lèi)必須保持并控制生態(tài)平衡，甚至必須控制人類(lèi)自身的增長(zhǎng)。首先我們建立兩個(gè)簡(jiǎn)單的單種群增長(zhǎng)模型，以簡(jiǎn)略分析一下這方面的問(wèn)題。

種群的數(shù)量本應(yīng)取離散值，但由于種群數(shù)量一般較大，為建立微分方程模型，可將種群數(shù)量看作連續(xù)變量,由此引起的誤差將是十分微小的。

世界人口數(shù)量統(tǒng)計(jì)數(shù)據(jù)：年1625183019301960197419871999人口億5102030405060中國(guó)人口數(shù)量統(tǒng)計(jì)數(shù)據(jù)：年1908193319531964198219902000人口3.04.76.07.210.311.312.951

馬爾薩斯（Malthus）模型

馬爾薩斯在分析人口出生與死亡情況的資料后發(fā)現(xiàn)，人口凈增長(zhǎng)率r基本上是一常數(shù)，（r=b-d,b為出生率，d為死亡率），因而提出了著名的人口指數(shù)增長(zhǎng)模型

。

分析與建模：人口的凈增長(zhǎng)率是一個(gè)常數(shù)，也就是單位時(shí)間內(nèi)人口增長(zhǎng)量與當(dāng)時(shí)人口數(shù)成正比。設(shè)t時(shí)刻人口數(shù)為N(t)，t=t0時(shí)，N(t0)=N0，則這個(gè)方程的解為：

馬爾薩斯模型的一個(gè)顯著特點(diǎn)：種群數(shù)量翻一番所需的時(shí)間是固定的。令種群數(shù)量翻一番所需的時(shí)間為T(mén)，則有：

故即Malthus模型模型檢驗(yàn)

比較歷年的人口統(tǒng)計(jì)資料，可發(fā)現(xiàn)人口增長(zhǎng)的實(shí)際情況與馬爾薩斯模型的預(yù)報(bào)結(jié)果基本相符，例如，1961年世界人口數(shù)為30.6（即3.06×109），人口增長(zhǎng)率約為2%，人口數(shù)大約每35年增加一倍。檢查1700年至1961的260年人口實(shí)際數(shù)量，發(fā)現(xiàn)兩者幾乎完全一致，且按馬氏模型計(jì)算，人口數(shù)量每34.6年增加一倍，兩者也幾乎相同。模型預(yù)測(cè)假如人口數(shù)真能保持每34.6年增加一倍，那么人口數(shù)將以幾何級(jí)數(shù)的方式增長(zhǎng)。例如，到2510年，人口達(dá)2×1014個(gè)，即使海洋全部變成陸地，每人也只有9.3平方英尺的活動(dòng)范圍，而到2670年，人口達(dá)36×1015個(gè)，只好一個(gè)人站在另一人的肩上排成二層了。故馬爾薩斯模型是不完善的。幾何級(jí)數(shù)的增長(zhǎng)Malthus模型實(shí)際上只有在群體總數(shù)不太大時(shí)才合理，到總數(shù)增大時(shí)，生物群體的各成員之間由于有限的生存空間，有限的自然資源及食物等原因，就可能發(fā)生生存競(jìng)爭(zhēng)等現(xiàn)象。所以Malthus模型假設(shè)的人口凈增長(zhǎng)率不可能始終保持常數(shù)，它應(yīng)當(dāng)與人口數(shù)量有關(guān)。2

Logistic模型

人口凈增長(zhǎng)率應(yīng)當(dāng)與人口數(shù)量有關(guān)，即：r=r(N)

從而有：（1）r(N)是未知函數(shù)，但根據(jù)實(shí)際背景，它無(wú)法用擬合方法來(lái)求。為了得出一個(gè)有實(shí)際意義的模型，我們不妨采用一下工程師原則。工程師們?cè)诮?shí)際問(wèn)題的數(shù)學(xué)模型時(shí)，總是采用盡可能簡(jiǎn)單的方法。r(N)最簡(jiǎn)單的形式是常數(shù)，此時(shí)得到的就是馬爾薩斯模型。對(duì)馬爾薩斯模型的最簡(jiǎn)單的改進(jìn)就是引進(jìn)一次項(xiàng)（競(jìng)爭(zhēng)項(xiàng)）此時(shí)得到微分方程：

或（2）

（2）被稱(chēng)為L(zhǎng)ogistic模型或生物總數(shù)增長(zhǎng)的統(tǒng)計(jì)籌算律，是由荷蘭數(shù)學(xué)生物學(xué)家弗赫斯特（Verhulst）首先提出的。一次項(xiàng)系數(shù)是負(fù)的，因?yàn)楫?dāng)種群數(shù)量很大時(shí)，會(huì)對(duì)自身增大產(chǎn)生抑制性，故一次項(xiàng)又被稱(chēng)為競(jìng)爭(zhēng)項(xiàng)。（2）可改寫(xiě)成：

（3）

(3)式還有另一解釋?zhuān)捎诳臻g和資源都是有限的，不可能供養(yǎng)無(wú)限增長(zhǎng)的種群個(gè)體，當(dāng)種群數(shù)量過(guò)多時(shí)，由于人均資源占有率的下降及環(huán)境惡化、疾病增多等原因，出生率將降低而死亡率卻會(huì)提高。設(shè)環(huán)境能供養(yǎng)的種群數(shù)量的上界為K（近似地將K看成常數(shù)），N表示當(dāng)前的種群數(shù)量，K-N恰為環(huán)境還能供養(yǎng)的種群數(shù)量，（3）指出，種群增長(zhǎng)率與兩者的乘積成正比，正好符合統(tǒng)計(jì)規(guī)律，得到了實(shí)驗(yàn)結(jié)果的支持，這就是（3）也被稱(chēng)為統(tǒng)計(jì)籌算律的原因。對(duì)（3）分離變量：兩邊積分并整理得：令N(0)=N0，求得：故（3）的滿(mǎn)足初始條件N(0)=N0的解為：（4）易見(jiàn)：

N(0)=N0

，N(t)的圖形請(qǐng)看右圖

模型檢驗(yàn)

用Logistic模型來(lái)描述種群增長(zhǎng)的規(guī)律效果如何呢？1945年克朗比克（Crombic）做了一個(gè)人工飼養(yǎng)小谷蟲(chóng)的實(shí)驗(yàn)，數(shù)學(xué)生物學(xué)家高斯（E·F·Gauss）也做了一個(gè)原生物草履蟲(chóng)實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果都和Logistic曲線十分吻合。

大量實(shí)驗(yàn)資料表明用Logistic模型來(lái)描述種群的增長(zhǎng)，效果還是相當(dāng)不錯(cuò)的。例如，高斯把5只草履蟲(chóng)放進(jìn)一個(gè)盛有0.5cm3營(yíng)養(yǎng)液的小試管，他發(fā)現(xiàn)，開(kāi)始時(shí)草履蟲(chóng)以每天230.9%的速率增長(zhǎng)，此后增長(zhǎng)速度不斷減慢，到第五天達(dá)到最大量375個(gè)，實(shí)驗(yàn)數(shù)據(jù)與r=2.309，a=0.006157，N(0)=5的Logistic曲線：

幾乎完全吻合，見(jiàn)右圖

Malthus模型和Logistic模型的總結(jié)

Malthus模型和Logistic模型均為對(duì)微分方程（1）所作的模擬近似方程。前一模型假設(shè)了種群增長(zhǎng)率r為一常數(shù)，（r被稱(chēng)為該種群的內(nèi)稟增長(zhǎng)率）。后一模型則假設(shè)環(huán)境只能供養(yǎng)一定數(shù)量的種群，從而引入了一個(gè)競(jìng)爭(zhēng)項(xiàng)。

用模擬近似法建立微分方程來(lái)研究實(shí)際問(wèn)題時(shí)必須對(duì)求得的解進(jìn)行檢驗(yàn)，看其是否與實(shí)際情況相符或基本相符。相符性越好則模擬得越好，否則就得找出不相符的主要原因，對(duì)模型進(jìn)行修改。

Malthus模型與Logistic模型雖然都是為了研究種群數(shù)量的增長(zhǎng)情況而建立的，但它們也可用來(lái)研究其他實(shí)際問(wèn)題，只要這些實(shí)際問(wèn)題的數(shù)學(xué)模型有相同的微分方程即可。3、車(chē)間空氣的清潔

問(wèn)題：已知一個(gè)車(chē)間體積為V立方米，其中有一臺(tái)機(jī)器每分鐘能產(chǎn)生r立方米的二氧化碳（CO2），為清潔車(chē)間里的空氣，降低空氣中的CO2含量，用一臺(tái)風(fēng)量為K立方米/分鐘的鼓風(fēng)機(jī)通入含CO2為m%的新鮮空氣來(lái)降低車(chē)間里的空氣的CO2含量。假定通入的新鮮空氣能與原空氣迅速地均勻混合，并以相同的風(fēng)量排出車(chē)間。又設(shè)鼓風(fēng)機(jī)開(kāi)始工作時(shí)車(chē)間空氣中含x0%的CO2.問(wèn)經(jīng)過(guò)t時(shí)刻后，車(chē)間空氣中含百分之幾的CO2?最多能把車(chē)間空氣中CO2的百分比降到多少？設(shè)t時(shí)刻（單位為分鐘）車(chē)間每立方米空氣含CO2的百分比為x(t)%,考慮時(shí)間區(qū)間并利用質(zhì)量守恒定律：內(nèi)車(chē)間空氣含CO2量的“增加”等于時(shí)間內(nèi)進(jìn)入的新鮮空氣中含CO2的量加上機(jī)器產(chǎn)生的CO2的量減去排出空氣中CO2的量。用數(shù)學(xué)公式表示出來(lái)就是分析和建模于是，令得其中，解為這就是t時(shí)刻空氣中含CO2的百分比。

通常否則含CO2的量只會(huì)增加。

令得這表明車(chē)間空氣中含CO2的量最多只能降到討論：如果設(shè)V=10000立方米，r=0.3立方米/分鐘，K=1500立方米/分鐘，m=0.04%，x0=0.12%。試問(wèn)：（1）需多少分鐘后，車(chē)間空氣中含CO2的百分比低于0.08%?（2）車(chē)間空氣中含CO2的百分比最多只能降到多少？

古尸年代鑒定問(wèn)題

在巴基斯坦一個(gè)洞穴里，發(fā)現(xiàn)了具有古代尼安德特人特征的人骨碎片，科學(xué)家把它帶到實(shí)驗(yàn)室，作碳14年代測(cè)定，分析表明，與的比例僅僅是活組織內(nèi)的6.24%，能否判斷此人生活在多少年前？小知識(shí)

年代測(cè)定：活體中的碳有一小部分是放射性同位素，這種放射性碳是由于宇宙射線在高層大氣中的撞擊引起的，經(jīng)過(guò)一系列交換過(guò)程進(jìn)入活組織內(nèi)，直到在生物體內(nèi)達(dá)到平衡濃度，這意味著在活體中，的數(shù)量與穩(wěn)定的的數(shù)量成定比，生物體死亡后，交換過(guò)程就停止了，放射性碳便以每年八千分之一的速度減少。背景設(shè)t為死后年數(shù)，年代測(cè)定的修訂：

1966年，耶魯實(shí)驗(yàn)室的Minze

Stuiver和加利福尼亞大學(xué)圣地亞哥分校的HansE.Suess在一份報(bào)告中指出：在2500到10000年前這段時(shí)間中測(cè)得的結(jié)果有差異，其根本原因在于那個(gè)年代，宇宙射線的放射性強(qiáng)度減弱了，偏差的峰值發(fā)生在大約6000年以前。他們提出了一個(gè)很成功的誤差公式，用來(lái)校正根據(jù)碳測(cè)定出的2300年到6000年前這期間的年代：

真正的年代＝4、導(dǎo)彈跟蹤問(wèn)題一、實(shí)驗(yàn)?zāi)康谋緦?shí)驗(yàn)主要涉及常微分方程。通過(guò)實(shí)驗(yàn)復(fù)習(xí)微分方程的建模和求解；介紹兩種微分方程的數(shù)值方法：Euler法和改進(jìn)的Euler法；還介紹了仿真方法。二、實(shí)際問(wèn)題某軍的一導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處海面上有敵艦一艘以90km/h的速度向正東方向行駛。該基地立即發(fā)射導(dǎo)彈跟蹤追擊敵艇，導(dǎo)彈速度為450km/h，自動(dòng)導(dǎo)航系統(tǒng)使導(dǎo)彈在任一時(shí)刻都能對(duì)準(zhǔn)敵艇。試問(wèn)導(dǎo)彈在何時(shí)何處擊中敵艦？三、數(shù)學(xué)模型微分方程建模的方法主要是依據(jù)守恒律來(lái)建立等量關(guān)系式。對(duì)于這個(gè)問(wèn)題，尋求等量關(guān)系是比較簡(jiǎn)單的。設(shè)坐標(biāo)系如圖3.1所示，取導(dǎo)彈基地為原點(diǎn)O(0,0)，x軸指向正東方,y軸指向正北方。當(dāng)t=0時(shí)，導(dǎo)彈位于點(diǎn)O，敵艇位于點(diǎn)A(0,H)，其中H=120（km）。設(shè)導(dǎo)彈在t時(shí)刻的位置為P(x(t),y(t))，由題意,其中

另外在t

時(shí)刻,敵艇位置為,其中。由于導(dǎo)彈軌跡的切線方向必須指向敵艇，即直線PM的方向就是導(dǎo)彈軌跡上點(diǎn)P的切線方向,故有或?qū)憺榉匠?3.1),(3.3)連同初值條件構(gòu)成了一個(gè)關(guān)于時(shí)間變量t的一階微分方程組的初值問(wèn)題。

為了尋求x與y的關(guān)系,要設(shè)法消去變量t,由式(3.2)得兩邊對(duì)t

求導(dǎo)即有把式(3.1)寫(xiě)為代入上式,就得到軌跡方程.這是一個(gè)二階非線性微分方程,加上初值條件,則初值問(wèn)題就是導(dǎo)彈軌跡的數(shù)學(xué)模型。值得注意的是，前面導(dǎo)出的一階微分方程組(3.1)，(3.3)和(3.4)實(shí)際上已經(jīng)是一個(gè)數(shù)學(xué)模型了，不過(guò)多一個(gè)變量（或說(shuō)參數(shù)）而已。四、解析方法方程（3.5）可以降階.令則（3.5）化為一階可分離變量方程

即易得

由初值條件(3.7)即

注意到上式可改寫(xiě)為于是有這樣我們又得到一個(gè)可分離變量方程積分得利用于是導(dǎo)彈軌跡方程為設(shè)導(dǎo)彈擊中敵艇于B(L,H),以Y=H代入(3.9)式,得而導(dǎo)彈擊中敵艇的時(shí)刻將數(shù)據(jù)代入(3.10),(3.11)式,得五、數(shù)值方法

將初值問(wèn)題(3.5)~(3.7)化為一階微分方程組取自變量y的步長(zhǎng)為,于是得分點(diǎn)相應(yīng)點(diǎn)上的x的值和p的值記為顯然,有初值條件我們將介紹兩種近似算法來(lái)進(jìn)行數(shù)值處理.ⅠEuler

方法

Euler方法十分簡(jiǎn)單,就是用差商代替微商,即以代之以而代之以這樣有設(shè)導(dǎo)彈到達(dá)處的時(shí)刻為那么得到計(jì)算的迭代格式于是表3.1是取n=4時(shí)的計(jì)算結(jié)果,讀者可以用來(lái)檢驗(yàn)程序或應(yīng)用軟件的正確性.

表3.1

kyk

xk

pk000013000.052601.50.123905.00.22412011.50.42此時(shí)表3.2是對(duì)于不同的n值所對(duì)應(yīng)的計(jì)算結(jié)果.顯然,n越大(即h越小),結(jié)果就越精確.

表3.2n4812244896120240L11.5215.9617.9720.5522.2523.3323.5824.15T0.1280.1770.2000.2280.2470.2590.2620.268

注意，由問(wèn)題（3.1）,(3.3),(3.4)消去t推導(dǎo)出問(wèn)題（3.5）～(3.7)是較為巧妙和偶然的.一般而言，一個(gè)微分方程組未必能消去一些變?cè)鴾p少方程的個(gè)數(shù)。

那么，我們能否直接對(duì)初值問(wèn)題(3.1),(3.3),(3.4)進(jìn)行數(shù)值處理呢？答案是肯定的。實(shí)際上，只要由方程(3.1),(3.3)解出和的表達(dá)式，這樣問(wèn)題變?yōu)槿r(shí)間步長(zhǎng)對(duì)應(yīng)時(shí)導(dǎo)彈軌跡上點(diǎn)的坐標(biāo)為則Euler格式為當(dāng)計(jì)算到即停止,于是,表3.3和表3.4分別列出了取步長(zhǎng)為0.1和0.05時(shí)的計(jì)算結(jié)果:

表3.3ktk

xk

yk

00.00.000000.0000010.10.0000045.0000020.25.3615489.6794630.322.67495131.21553此時(shí)取表3.4ktk

xk

yk

10.050.0000022.5000020.101.0373642.9760730.153.4120567.3504140.207.6461589.4484350.2514.86790110.7579660.3029.19480128.10702此時(shí)取表3.5是對(duì)應(yīng)不同的,用Euler法所得相應(yīng)的步長(zhǎng)推進(jìn)次數(shù)n和計(jì)算結(jié)果.

表3.50.10.050.0050.001n3656278L22.6749529.1948025.6673125.04935T0.251940.324390.285190.27833

Euler方法較為簡(jiǎn)單,但也較為粗糙,對(duì)形式較復(fù)雜的微分方程更易有較大的誤差.人們?cè)O(shè)計(jì)了不少更精確的近似算法,這里我們介紹其中的一種,進(jìn)一步研究可看參考書(shū).

Ⅱ.改進(jìn)的Euler方法(預(yù)報(bào)---校正法)

以一維情況為例,對(duì)問(wèn)題Euler迭代格式是其中而改進(jìn)的Euler迭代格式則是其中由積分表達(dá)式的幾何意義看,右邊為下方的曲邊梯形,從圖3.2我們可以看出Euler法是用矩形來(lái)代替曲邊梯形,而改進(jìn)的Euler法則是用梯形來(lái)代替曲邊梯形.

對(duì)問(wèn)題(3.18)~(3.20),我們寫(xiě)出相應(yīng)的改進(jìn)Euler迭代格式表3.6和表3.7分別列出了取步長(zhǎng)為0.1和0.05時(shí)的計(jì)算結(jié)果:

表3.6ktk

xk

yk

00.00.000000.0000010.12.6807744.8397320.212.5752488.2867930.322.07242130.25569此時(shí)取表3.7ktk

xk

yk

10.050.5186822.4880420.102.1059644.9219530.155.0982167.2021340.2010.1609689.0690650.2519.65646108.9789860.3024.24089130.99030此時(shí)取表3.8是對(duì)應(yīng)不同的,用改進(jìn)的Euler法所得相應(yīng)的步長(zhǎng)推進(jìn)次數(shù)n和計(jì)算結(jié)果.

表3.80.10.050.0050.001n3656278L27.0724224.2408925.1355224.98112T0.300800.269340.279280.27757

圖3.3畫(huà)出了導(dǎo)彈軌跡由解析式所給出的精確曲線以及由Euler法和改進(jìn)的Euler法進(jìn)行數(shù)值計(jì)算所給出的近似曲線.六、仿真方法

如果建立微分方程很困難，或者微分方程很復(fù)雜而較難作出數(shù)值處理，

常?？梢杂梅抡娣椒?。所謂仿真方法，顧名思義，指的是模仿真實(shí)行為和過(guò)程的方法。在這個(gè)具體問(wèn)題中，就是一步步地模擬導(dǎo)彈追蹤敵艇的實(shí)際過(guò)程。

而計(jì)算機(jī)仿真，則是在計(jì)算機(jī)上通過(guò)相應(yīng)的程序和軟件來(lái)實(shí)現(xiàn)對(duì)事件運(yùn)行的實(shí)際過(guò)程的模擬。設(shè)導(dǎo)彈和敵艇在初始時(shí)刻（即t=0時(shí)）分別位于P0(0,0)和M0(0,H)。此時(shí)，導(dǎo)彈指向M0。而在

t=時(shí)，導(dǎo)彈的位置P1(x1,y1)，其中，敵艇的位置則為