1、第四章 微分方程建模方法4.1 導言在許多實際問題中，例如物理中的速率問題，人口的增長問題，放射性衰變問題，經(jīng)濟學中的邊際問題等，常常涉及到兩個變量之間的變化規(guī)律。微分方程是研究上述問題的一種機理分析方法，它在科技、工程、生態(tài)、環(huán)境、人口以及經(jīng)濟管理等領域中有著十分廣泛的應用。在應用微分方程解決實際問題時，必須經(jīng)過兩個階段。一是微分方程的建立，建立一個微分方程的實質就是構建函數(shù)、自變量以及函數(shù)對自變量的導數(shù)之間的一種平衡關系。而正確地構建這種平衡關系，需要對實際問題的深入淺出的刻畫，根據(jù)物理的和非物理的原理、定律或定理，作出合理的假設和簡化并將它升華成數(shù)學問題。另一個是方程的求解和結果分析。對

2、一些常系數(shù)的或特殊函數(shù)形式的微分方程，往往能得到解析解，這對實際問題的分析和應用都是有利的,但是大多數(shù)變系數(shù)的、非線性函數(shù)形式的微分方程都是求不出解析解的,此時就需要應用求解微分方程的另一個重要方法數(shù)值解法。本章簡要介紹有關微分方程模型的概念，微分方程的數(shù)值解法和圖解法，主要介紹若干建模實例，通過它們展示微分方程模型的建模步驟及解決實際問題的全過程。4. 微分方程模型及其解法4.2.1微分方程模型由物理意義可知，速率就是變化率，而變化率即數(shù)學中的導數(shù)。如果一個實際問題中涉及到一個變量的變化速率，那么就有可能用微分方程模型來描述該問題。例如某國人口的增長率與其人口的總數(shù)成正比，現(xiàn)需要預測未來某年

3、該國的人口總數(shù)。這就是一個涉及到變量的變化速率（人口的增長率）的實際問題。要解決該問題，就必須將它轉化為數(shù)學問題。如果設表示時間，表示時刻的人口總數(shù)，假設是的可導函數(shù)，則人口增長率就是，根據(jù)假設可建立微分方程一般地，包含變量、函數(shù)及其導數(shù)的等式稱為微分方程。微分方程的一般形式是或在上述問題中，如果已知時刻該國的人口總數(shù)為，則這個人口預測問題就可完整地表示為 (4.1)關系式(4.1)就是上述人口預測問題的數(shù)學模型.可以求得（.1）的解析解。下面我們考察葫蘆形狀容器壁上的刻度問題.對于圓柱形狀容器,由于體積與相對于容器底部的任意高度的函數(shù)關系明確:其中直徑是常數(shù),因此可以非常方便地在容器壁上標示

4、容積刻度.但是對于非簡單規(guī)則幾何形狀的容器,在其壁上標示容積刻度則比較困難.現(xiàn)在假設有一葫蘆形狀的容器,人們通過測量得到高度與直徑的離散數(shù)據(jù)如表4-1.表4-100.20.40.60.81.01.20.640.851.000.580.200.120.12 根據(jù)表4-1中的數(shù)據(jù),可以大致畫出該容器的狀,如圖4.1所示。0圖4.1 葫蘆形狀容器 如何標出任意高度處容器體積的刻度？由微元法分析知 其中是任意高度，直徑是高度的函數(shù)，記為.因此 （.）只要求解上述微分方程，就可以求出體積與高度之間的函數(shù)關系，從而標出任意高度處容器體積。（4.2）就是上述物理問題的數(shù)學模型。但是，由于無解析表達式，故無法

5、求得（.2）的解析解。這就需要利用微分方程的數(shù)值解法。下節(jié)我們主要簡述關于微分方程求解的系列方法。4.2.2微分方程解法 在高等數(shù)學中，介紹了微分方程的基本知識以及一些比較特殊類型的微分方程的解析解法，但是大量的復雜微分方程，無法或難于求出其解析解。因此，在實際問題中經(jīng)常需要用數(shù)值解法求微分方程和微分方程組。為了數(shù)學建模的需要，本節(jié)我們主要補充微分方程的數(shù)值解法。設有微分方程問題 （.）數(shù)值解法的基本方法是:根據(jù)需要給定自變量點列：定義步長為.常給定等間距自變量點列,此時步長為常數(shù)。數(shù)值解法一般只能求得微分方程的近似解。設微分方程(4.3)在點列上的精確解是,而其近似解為。根據(jù)一定的原理,用差

6、商代替導數(shù),結合初始條件,推出計算的迭代公式。1 歐拉方法當較小時，在小區(qū)間上用差商代替導數(shù)，即用近似代替方程（4.3）左端導數(shù)，而方程右端函數(shù)中取,由方程(4.3)得近似表達式 (4.4)在(4.4)中近似地取,則由(4.4)得 (4.5)其中是已知的初始點。由(4.5)和已知的初始點,我們就可以一步一步地推算出:因此,我們稱公式（4.5）為顯式歐拉公式。如果在上述過程中，在中取，那末可得另一個歐拉公式 (4.6)大多數(shù)情況下是非線性函數(shù),從而公式（4.6）是非線性方程,無法直接計算,因此,我們稱公式（4.6）為隱式歐拉公式。一般而言,歐拉方法計算精度低,收斂速度慢。如果把顯式歐拉公式（4.

7、5）和隱式歐拉公式（4.6）加以平均，得到 （4.7）我們稱公式（4.7）為梯形公式。可以證明，梯形公式比歐拉公式精度高,收斂速度快。例4.1 設微分方程問題分別用顯式歐拉法、隱式歐拉法、梯形法求解，取步長和0.001，并利用其解析解（精確解），對數(shù)值解誤差進行分析比較。解 這是一個一階線性非齊次微分方程，用解析解法得到的精確解是下面用數(shù)值解法：取h=0.1時，1） 顯示歐拉法的迭代公式為： ， ，.2） 隱式歐拉法的迭代公式為： ， ， 將它變形為： ， 3） 梯形法的迭代公式為：， 將它變形為：， 計算結果列入表4-2。表4-2 解析解及時三種數(shù)值法計算結果解析解顯式歐拉解隱式歐拉解梯形解

8、011110.11.004811.00911.00480.21.01871.01001.02641.01860.31.04081.02901.05131.04060.41.07031.05611.08301.07010.51.10651.09051.12091.10630.61.14881.13141.16451.14850.71.19661.17831.21321.19630.81.24931.23051.26651.24900.91.30661.28741.32411.306311.36791.34871.38551.3673取時，容易得到與上述類似的公式，其計算結果列入表4-3。表4-3

9、 解析解及時三種數(shù)值法計算結果解析解顯式歐拉解隱式歐拉解梯形解011110.11.00481.00481.00491.00480.21.01871.01861.01881.01870.31.04081.04071.04091.04080.41.07031.07021.07051.07030.51.10651.10641.10671.10650.61.14881.14861.14901.14880.71.19661.19641.19681.19660.81.24931.24911.24951.24930.91.30661.30641.30681.306611.36791.36771.36811.

10、3679由表4-2和表4-3首先可見，取步長時，計算結果具有兩位有效數(shù)字，而 時，計算結果具有四位有效數(shù)字。（設為一準確值，為的一個近似值，則稱為的誤差。若的誤差不超過的某一位上的半個單位，而且從該位到的左邊第一位非零數(shù)字共有位，則稱具有位有效數(shù)字 ）一般說來，在迭代中，步長越小，計算結果越精確（當步長在數(shù)值方法的絕對穩(wěn)定區(qū)域內時）。其次可以發(fā)現(xiàn)，計算點離開初始點越遠，的誤差越大。最后還可知，梯形法顯然優(yōu)于顯式、隱式歐拉法，例如取步長時，梯形法的計算結果基本具有四位有效數(shù)字（僅在處，只有三位有效數(shù)字）。從上例看好象用梯形法計算非常容易，其實不然，梯形公式是一種隱式公式，一般情況下計算比較困難。

11、于是產(chǎn)生了一種改進的方法：第一步，由顯式歐拉公式（4.5）算出的預測值，第二步，將代入梯形公式（4.7），進行校正，即 （4.8）我們稱公式（4.8）為改進的歐拉公式。在程序編寫中，它常常寫為 （4.9）改進的歐拉公式是顯式公式，易于計算，精度較高，收斂速度快，是人們最常用的方法之一。歐拉法和改進的歐拉法還可以推廣到微分方程組的情形。例4.2 用改進歐拉公式求例4.1中的初值問題。解 由公式（4.8）知由于故化簡可得計算結果列于表4-4。表4-400.10.20.30.40.511.005001.0190251.0412181.0708021.107076對比例4.1與例4.2可見改進歐拉公式

12、的精確性高于歐拉公式，與梯形公式相當。2.龍格-庫塔方法利用泰勒公式可以得到龍格-庫塔方法（簡稱方法）是由上式產(chǎn)生迭代公式進行計算。若則稱以上迭代公式為階公式，越大，截斷誤差越小，精度越高。要得到一個階公式，關鍵是選擇函數(shù)，使之滿足截斷誤差為的要求。 經(jīng)常使用的公式有1） 2階公式 中點公式 （4.10）2）3階公式 （4.11）3）4階公式 經(jīng)典龍格-庫塔公式 （4.12）關于微分方程組也有類似的數(shù)值解法。在MATLAB軟件中含有數(shù)值求解的系統(tǒng)函數(shù)，其實現(xiàn)原理就是龍格-庫塔方法。例4.3 用經(jīng)典龍格-庫塔法求初值問題在處的解。 解 ，由公式（4.12）得 ， ，本題的精確解是，現(xiàn)將計算結果列

13、表如下：表4-500.20.40.60.81.011.183231.341671.483281.612511.7321411.183221.341641.483241.612451.73205 下面我們簡要介紹圖解法。 高等數(shù)學中介紹的解析解法，雖然可以得到精確解，但是應用范圍狹窄。上面介紹的數(shù)值解法適應面廣，雖然得到的是離散點上解的近似值，但可以控制其誤差，因此得到了廣泛地應用。下面簡要介紹微分方程的圖解法。該法的特點是可以將微分方程解的全局信息直觀地、形象地展現(xiàn)出來。如果一個一階微分方程可以寫成如下形式其中是已知的二元函數(shù)，那么我們就能夠確定積分曲線（即解函數(shù)）在任意點處的斜率。從圖象上看

14、，給出平面上一個點，我們就可以畫出一條通過該點的短直線，這條短直線的斜率等于積分曲線在該點處準確的斜率。適當設定在坐標平面上一定數(shù)量的點處畫這種短直線以及這種短直線的長度，這種短直線布滿整個坐標面所形成的圖形稱為斜率場。另一方面如果我們求出了微分方程的通解那么我們就得到一蔟積分曲線，不同的初始值對應于不同的解曲線。在這些曲線蔟中任意取定一條，進行如下處理：1） 曲線分段: 將曲線分割成若干弧線段。2） 以直代曲： 將每一弧線段用其左端點處的短切線段代替。3） 相鄰不連： 短切線段足夠地短，相鄰短切線不相連。經(jīng)過上述處理，由微分方程的解可以得到斜率場，反之，由微分方程直接得到斜率場后（不需求微分

15、方程的解析解），從斜率場可以隱約看出積分曲線的形狀。可以將斜率場看作一簇路標，它指示短切線左端點處要走的方向。要得到微分方程初值問題的解曲線，就可以從坐標平面上的初始點出發(fā)沿該點斜率場所指示的方向前進，按斜率場所指示的走向一段一段地走下去，就會走出一條“路”。這條“路”就是微分方程的積分曲線在坐標平面上的圖形。大家知道圖形也是函數(shù)的一種表現(xiàn)形式，因此，能畫出上述“路”就相當于確定了解函數(shù)。 4. 新產(chǎn)品的推銷模型一種新產(chǎn)品上市后，經(jīng)銷者自然十分關心它的銷售情況，尤其是銷售一段時間后，經(jīng)銷者往往需要根據(jù)已有銷售情況，預測該產(chǎn)品在本地區(qū)的總銷售量，從而恰當?shù)亟M織貨源。下面我們根據(jù)兩種不同的假設建立

16、兩種產(chǎn)品銷售速度模型。第一種模型 假設產(chǎn)品是以自然推銷的方式賣出，即被賣出的產(chǎn)品起著實物性廣告宣傳作用，自然地吸引著暫未購買的潛在消費者。設在時刻已售出的產(chǎn)品總數(shù)為，并假設售出的每一產(chǎn)品在單位時間內平均可以吸引個顧客，使其購買個該產(chǎn)品，則滿足微分方程 （4.13）設初始條件為 （4.14）上述微分方程的解是 （4.15）此即為著名的模型。關于這個模型，我們作如下幾點說明：（1）經(jīng)過與實際銷售情況對比，發(fā)現(xiàn)（4.15）式的計算結果與實際銷售量在初始階段的增長情況比較符合。（2）在產(chǎn)品即將賣出、尚未賣出之時刻為，顯然，這時由（4.15）式知，這一結果自然與事實不符。產(chǎn)生這一錯誤結果的原因在于我們假

17、定產(chǎn)品是自然推銷的，然而，自然推銷方式無法自動啟動，因為在最初產(chǎn)品還未賣出之時，在消費者中還沒有作為實物性廣告的產(chǎn)品，從而無法吸引消費者來購買。為此，我們這樣來解決上述矛盾。將第一個單位時間內看作是以非自然推銷方式進行銷售（商家做廣告），而將第一個單位時間之后看作以自然推銷方式進行銷售，因此，我們將時刻的銷售情況作為初始條，即 （4.16）微分方程初值問題（4.13）、（4.16）的解為 （4.17）(4) 由（4.15）或（4.17）得，一般對耐用消費品而言，這顯然與事實不符。事實上，耐用消費品，比如冰箱、電視機、電腦和汽車，需求量總是有限的，因此往往是有上界的。從上面的分析可見，模型不宜用

18、于銷售量的中、長期預測。下面我們來介紹一種有界模型。第二種模型 假設：1）需求量的上界為，2）經(jīng)銷商可以通過自然推銷方式和其它方式推銷產(chǎn)品，3）一個消費者僅購買一件該種產(chǎn)品。根據(jù)這三個假設，在某時刻時，產(chǎn)品銷量的增長既與到時刻為止的已經(jīng)購買該種產(chǎn)品消費者數(shù)目成正比，也與尚未購買該種產(chǎn)品的潛在消費者數(shù)目成正比，即設比例系數(shù)為，則 （4.18）注意到初始條件用分離變量法可求得上述微分方程的解 （4.19）0上式為模型。當時如果則由（4.19）式得到另外，在（4.19）中令（見圖4.2）。可見上述模型是有界的。圖4.2由（4.19）式知， 即是關于時刻的單調函數(shù)，這與實際情況吻合，賣出的產(chǎn)品累積量自

19、然越賣越多。此外，對（4.19）式兩端關于求導，得故令，得到。如圖4.2，當時，由函數(shù)圖象為上凹弧；同理，當時，單調減小，函數(shù)圖形為上凸弧。這說明，在銷售量小于最大需求量的一半時，銷售速度是不斷增加的；銷售量達到最大需求量的一半時，該產(chǎn)品正處于最暢銷時刻，此后銷售速度開始下降。實際情況表明，產(chǎn)品銷售情形與模型十分相似，特別在銷售后期更加吻合。如果時，可用初始條件此時可得解 （4.20）問題一 某公司生產(chǎn)的冰箱在某個城鎮(zhèn)10年的銷售記錄為（單位：萬臺）年序號12345678910銷量量1.01.02.25.57.58.59.08.57.85.7試預測該種冰箱在該城鎮(zhèn)的總銷售量。解 設在時刻已售出

20、的冰箱總數(shù)為，滿足模型（4.20），即此處需要求出最大需求量。由的定義知， 一般地，可得，用差商近似代替導數(shù)，注意到，可將微分方程化為處（）的差分方程即寫成矩陣形式為記,上式可寫為這是一個超定線性代數(shù)方程組,在上式兩端左乘得如果的逆存在，用左乘上式兩端，就得到將具體數(shù)據(jù)代入得 ， ， ， ， 故。因此，我們預測該種冰箱在該城鎮(zhèn)的總銷售量為69.29萬臺。另外，我們還可利用模型預測以后某年的該種冰箱在該城鎮(zhèn)的銷售量。4.4抵押貸款購房問題問題 隨著改革開放的不斷深入，我國開始實行住房貨幣化分配制度，一般人家都無力用現(xiàn)金買下自己滿意的住房，從而面臨貸款購房問題。下面是1991年1月1日某大城市晚報

21、上刊登的一則商品房廣告。 名流 名 流 花 園 之 高 尚 住 宅 公 寓 花園 首 期 隆 重 發(fā) 售 用薪金，買高品質住房 對于大多數(shù)工薪階層的人士來說，想買房，簡直是天方夜談。現(xiàn)在有這樣一棟：自備款只需七萬元人民幣，其余由銀行貸款，分五年還清。相當于每月只需付1200元人民幣。那么，這對于您還有什么問題呢？ 任何人看了這則廣告都會產(chǎn)生許多疑惑，比如廣告中沒有說明住房面積、房型、設施以及所處地段等等，下面我們的數(shù)學建模僅針對人們更關心問題：如果一次付款買這棟房要多少錢呢？銀行貸款的利息是多少呢？為什么每月要付1200元呢？是如何算出來的？因為人們都知道，如果知道了一次付款買房的價格，自己暫

22、時只能支付一部分，其余的房款就不得不找銀行貸款，只要知道了利息，就應該可以算出五年內每月應付多少錢才能按時還清貸款了，從而就可以對是否要去買該廣告中所說的房子作出決策。現(xiàn)在我們分四個部分討論。一 明確變量、參數(shù)，顯然下面的量是需要考慮的： 表示需要貸款量； 表示月利息（貸款通常按復利計）； 表示借期，按月數(shù)計；二 建立變量之間的數(shù)學關系。若用表示第個月時尚欠的款數(shù)，則一個月后，本錢加上利息為，假如每月還款元，則第個月時尚欠的款數(shù)為因此，我們得到如下數(shù)學模型 （4.21）三 數(shù)學模型（4.21）的求解。由 應用數(shù)學歸納法，易知 ， （4.22） 四 根據(jù)廣告中的信息，我們可知，年個月，每月還款元

23、，但是需要貸款數(shù)（即一次性付款購買價減去首次付款70000萬元后剩下的款數(shù)）和銀行貸款利率，廣告未予說明，這造成了我們決策的困難。然而，由（4.22）可知60個月還清，即，從而得從中將解出的得 （4.23）（4.23）表示給定時與之間的關系式，如果我們已知銀行的貸款利率，就可以算出。例如，若，則由（4.23）可算得元。因此，我們就計算出了該房地產(chǎn)公司該種房的一次性付款的房價應該是70000+53946=123946元。如果銀行對房地產(chǎn)公司和個人所發(fā)放的貸款，其貸款利率相同，都是，如果該房地產(chǎn)公司說一次性付款的房價小于123946元的話，就說明找房地產(chǎn)公司分期付款買房不劃算，你就應該自己銀行貸款

24、。利用（4.22）式，以下我們進一步討論兩個問題。例1 某單位一對年青夫婦為買房要用銀行貸款60000元，月利率0.01，貸款期25年=300月，這對夫婦希望知道每月要還多少錢，25年就可還清。假設這對夫婦每月可有節(jié)余900元，是否可以去買房呢？解： 根據(jù)上面的討論，現(xiàn)在的問題就是要求使的，由（4.22）式知將代入，算得元，這說明這對夫婦有能力貸款買房。例2 接著例1的背景討論，這對夫婦又看到某借貸公司的一則廣告：“若借款60000元，22年還清，只要（i）每半個月還316元；(ii)由于增加了文書工作量，你要預付三個月的款，即6323=1896元。”這對夫婦算了一下：每半個月跑一趟去交款31

25、6元，一個月正好交款632元，雖然預付三個月的款1896元，這使人不快，但提前三年還清可省下632123=22752元，這是1896元的12倍！面對這一計算結果，這對夫婦又產(chǎn)生了一個疑問：這家借貸公司是慈善機構呢還是仍然要賺我們的錢呢？現(xiàn)在我們來給他們一個滿意的回答。解 我們先就（i）,(ii)兩條作一個初步分析，先孤立地分析（i）和(ii)看看能否縮短還款期。分析（i），這時不變，月利率變?yōu)榘朐吕剩纱致缘卣J為半月利率正好是月利率的一半，即，于是由（4.22）式可求出使的還款期，即由求得 （4.24）將，代入（4.24）式得（半月）年，即只能提前大約一個月還清貸款。由此可見，該借貸公司如果

26、在廣告中只有條件（i）的話，它一定會虧本，而這是不可能的，除非它是一個慈善機構。分析（ii）,此時元，即你只借了58104元而不是60000元，按問題1中銀行貸款的條件算一下，即令（月息），（每月還款），求使的，來看看能否提前還清。將，代入（4.24）得（月）（年），這表明實際上不需通過借貸公司，直接向銀行借貸，就可以提前接近四年還清貸款。因此，該借貸公司只要去同樣的銀行借款，即使半個月收來的316元不動再過半個月合在一起去交給銀行，它還可以坐收第22年的用戶交款（22-21.09）（年）12（月）632（元）=6901.44（元），更何況它可以利用每上半個月收到的還款去做短期（半個月內）投資

27、賺取額外的錢。根據(jù)以上分析，我們現(xiàn)在可以告訴這對夫婦，這家借貸公司不是慈善機構，它還是要賺你們的錢的。所以你們應該直接找銀行貸款。4.5 利用微跡元素反演膏鹽的沉積環(huán)境在天然地質體系中，微跡元素在平衡共存的兩相間的分配服從能斯特(Nernst)分配定律,即當兩種相處于化學平衡時，某種元素i在兩相間的濃度比為常數(shù)Kd。在一定的濃度范圍內的Kd與i的濃度無關，當兩種相均為聚集相時，Kd受壓力影響較小，而與溫度的關系較為顯著。本文試圖根據(jù)微跡元素在溶液中由于蒸發(fā)而結晶沉淀過程中服從能斯特分配定律，從而建立沉積盆地演化環(huán)境的數(shù)學模型。一、數(shù)學模型的建立蒸發(fā)巖層中鹽巖（NaCl）總是含有微量的Br。這是

28、由于Br取代鹽巖晶格中的Cl。Braitsch和Herrmann（1963）測定了Br在鹽巖與海水之間的分配系數(shù)Kd，發(fā)現(xiàn)在25鹽巖沉淀開始時，Kd0.15。這意味著在鹽巖沉積過程海水中的Br含量就迅速增長，隨著溴化物濃度增高，Kd值可近似地保持不變。 因為在溶液中，只有當溶質的含量超過它的溶解度時，才發(fā)生結晶沉淀。因此可以假設在沉淀過程中，海水中NaCl溶液總是飽和的，沉淀的原因是由于海水的蒸發(fā)而引起的。假設在沉淀發(fā)生過程中溫度不發(fā)生變化，則1m3飽和的海水中含NaCl為（0.36kg），重量為11.36kg。、都可以看成是固定不變的。另一方面，假設海水的蒸發(fā)量為 (m/a)，海水的體積為

29、(m3)，海水中Br的濃度為 (ppm)，、是沉積開始時海水的體積和海水中Br的濃度，假設注入盆地的注水量為 (m/a)，其中1m3的注入水中含NaCl為 (kg)，注入水中Br的濃度為，為盆地的面積。 從時刻到時刻，由于海水的蒸發(fā)而晶出NaCl晶體，從而使海水中Br的含量減少了 (4.29)同時，由于注入了非飽和的含有NaCl的補充水，上式可以改寫為 (4.30)式（4.30）滿足的條件是。即海水的蒸發(fā)量要大于使注入水成為飽和溶液的量。 另一方面，從時刻到時刻，海水中Br的含量減少 (4.31)因此有 (4.32)令，有 (4.33)式（4.33）是的一階微分方程令 (4.34) (4.35)則 (4.36) 如果盆地每年蒸發(fā)率與注水率均為常數(shù)，計為與，那么可以計算出盆地海水的體積為 ， (4.37)則（4.36）式可以簡化為： 令 ， 當時 (4.38.1) 當時 (4.38.2)由此可以計算出NaCl晶體中Br的濃度 (4.39) 另一方面，由于海水的蒸發(fā)，NaCl結晶沉積，其沉積厚度可表示為 (4.40) (4.41)其中為晶體NaCl的密度(kg/m)。 從(4.38)、(4.39)、(4.41)可以消去時間t得出沉積盆地中晶體NaCl中Br的濃度（或含量）與其沉積厚度的關系式。 二、實例 假設某膏鹽沉積盆地中晶體NaC