微分方程模型

微分方程模型第一頁，共四十三頁，2022年，8月28日第3講微分方程建模方法及案例在研究實際問題時，常常會聯系到某些變量的變化率或導數，這樣所得到變量之間的關系式就是微分方程模型。微分方程模型反映的是變量之間的間接關系，因此，要得到直接關系，就得求微分方程。第二頁，共四十三頁，2022年，8月28日求解微分方程的方法求解微分方程有三種方法1）求精確解；2）求數值解（近似解）；3）定性理論方法。第三頁，共四十三頁，2022年，8月28日建立微分方程模型的方法（1）根據規(guī)律列方程利用數學、力學、物理、化學等學科中的定理或經過實驗檢驗的規(guī)律等來建立微分方程模型。（2）微元分析法利用已知的定理與規(guī)律尋找微元之間的關系式，與第一種方法不同的是對微元而不是直接對函數及其導數應用規(guī)律。第四頁，共四十三頁，2022年，8月28日（3）模擬近似法在生物、經濟等學科的實際問題中，許多現象的規(guī)律性不很清楚，即使有所了解也是極其復雜的，建模時在不同的假設下去模擬實際的現象，建立能近似反映問題的微分方程，然后從數學上求解或分析所建方程及其解的性質，再去同實際情況對比，檢驗此模型能否刻畫、模擬某些實際現象。第五頁，共四十三頁，2022年，8月28日3.1年代鑒定問題在巴基斯坦一個洞穴里，發(fā)現了具有古代某種動物特征的骨碎片，科學家把它帶到實驗室，作碳14年代測定，分析表明，與的比例僅僅是活組織內的6.24%，能否判斷此動物生活在多少年前？第六頁，共四十三頁，2022年，8月28日

年代測定：活體中的碳有一小部分是放射性同位素，這種放射性碳是由于宇宙射線在高層大氣中的撞擊引起的，經過一系列交換過程進入活組織內，直到在生物體內達到平衡濃度，這意味著在活體中，的數量與穩(wěn)定的的數量成定比，生物體死亡后，交換過程就停止了，放射性碳便以每年八千分之一的速度減少。背景第七頁，共四十三頁，2022年，8月28日設t為死后年數，第八頁，共四十三頁，2022年，8月28日用Matlab求解對數log(x)x的自然對數log2(x)x的以2為底的對數log10(x)x的以10為底的對數以a為底x的對數，根據換底公式可表示為：log(x)/log(a)第九頁，共四十三頁，2022年，8月28日Matlab中的求解方法-8000*log(0.0624)ans=2.2194e+004第十頁，共四十三頁，2022年，8月28日年代測定的修訂

1966年，耶魯實驗室的MinzeStuiver和加利福尼亞大學圣地亞哥分校的HansE.Suess在一份報告中指出：在2500到10000年前這段時間中測得的結果有差異，其根本原因在于那個年代，宇宙射線的放射性強度減弱了，偏差的峰值發(fā)生在大約6000年以前。他們提出了一個很成功的誤差公式，用來校正根據碳測定出的2300年到6000年前這期間的年代：

真正的年代＝第十一頁，共四十三頁，2022年，8月28日3.2人口演化模型第十二頁，共四十三頁，2022年，8月28日預備知識:最小二乘原理當由實驗提供了大量數據時,不能要求擬合函數在數據點處的偏差,即(i=1,2,…,m)

嚴格為零,但為了使近似曲線盡量反映所給數據點的變化趨勢,需對偏差有所要求.通常要求偏差平方和最小,此即稱為最小二乘原理第十三頁，共四十三頁，2022年，8月28日Matlab的最優(yōu)化工具箱中提供了求非線性最小二乘擬合的函數：lsqcurvefit，調用格式為：

x=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數F(x,xdata)的M-文件,自變量為x和xdata迭代初值已知數據點用Matlab作非線性最小二乘擬和注意:1.如無合理初值，那就只能給出一個猜想初值2.擬合結果是初值敏感的，因為找到的不一定是全局最優(yōu)而可能是初值附近的局部最優(yōu).第十四頁，共四十三頁，2022年，8月28日

例

用下面一組數據擬合中的參數a，b，k該問題即解最優(yōu)化問題：第十五頁，共四十三頁，2022年，8月28日

1）編寫M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('curvefun1',x0,tdata,cdata)解法第十六頁，共四十三頁，2022年，8月28日問題的提出人口的增長情況是當前世界引起普遍關注的問題，早在18世紀人們就開始進行人口預報工作了。20世紀90年代，我們經常可以在報刊上看見關于人口增長的預報，說到20世紀末，或21世紀中葉，全世界（或某地區(qū)）的人口將達到多少多少億。這些人口預報的數值是從哪里來的？準確不準確？你能不能對某地人口數目的演化進行一下估算？第十七頁，共四十三頁，2022年，8月28日年（公元）人口（百萬）17903.918005.318107.218209.6183012.9184017.1185023.2年（公元）人口（百萬）186031.4187038.6188050.2189062.9190076.0191092.01920106.5年（公元）人口（百萬）1930123.21940131.71950150.71960179.31970204.01980226.51990251.4下表是近兩百年的美國人口統(tǒng)計數據，試依此建立美國人口增長的數學模型，最后用它預報2000年、2010年美國人口.第十八頁，共四十三頁，2022年，8月28日問題的分析從整體來說，人口的變化由兩個因素決定：出生和死亡。出生使得人口增加，死亡使得人口減少。對于局部地區(qū)來說，除了出生和死亡外，影響人口的變化還有兩個因素：遷入和遷出。遷入使得局部人口增加，遷出使得局部人口減少。在遷入、遷出人口的差別不大時，人口的變化主要由出生率和死亡率決定。根據上面的分析，不難建立起人口演化模型。第十九頁，共四十三頁，2022年，8月28日求解步驟[1]假設[2]建立模型[3]模型求解[4]模型的參數估計第二十頁，共四十三頁，2022年，8月28日[1]

假設人口增長率r(即出生率b-死亡率d)是常數（或單位時間內人口的增長量與當時的人口成正比）.第二十一頁，共四十三頁，2022年，8月28日[2]

建立模型記時刻t=0時人口數為x0,時刻t的人口為x(t),由于量大,可視為連續(xù)、可微函數.t到時間內人口的增量為：于是x(t)滿足微分方程：

(1)第二十二頁，共四十三頁，2022年，8月28日[3]

模型求解解微分方程得

(2)

(2)被稱為指數增長模型——馬爾薩斯提出(1798)，表示隨著時間增加，人口按指數規(guī)律無限增長被稱為等比數列模型第二十三頁，共四十三頁，2022年，8月28日[4]

模型的參數估計要用模型的結果來預報人口，必須對其中的參數r進行估計，這可以用表中所給的數據通過擬合得到.

第二十四頁，共四十三頁，2022年，8月28日擬合的關鍵問題1.t的變換指數函數exp(t)當t很大時可能會溢出,為了減小數據誤差,首先將時間域變換至[0,20],所用的變換為:t=1800+(t-1800)/102.x0和r初值的確定x0初始值自然應取t=0時的x的值3.9r初始值取增長率的平均值mean(diff(x)./diff(t)./x(1:20))第二十五頁，共四十三頁，2022年，8月28日Matlab求解過程首先編寫M文件fun.m，其中參數c(1)表示x0,c(2)表示r。functionx=fun(c,t)x=c(1)*exp(c(2)*t);第二十六頁，共四十三頁，2022年，8月28日再編寫一個M文件clearall;closeall;t=0:1:20;x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5251.4];plot(t,x,'o');holdon;c(1)=3.9;c(2)=mean(diff(x)./diff(t)./x(1:20));k=lsqcurvefit(@fun,c,t,x)tt=[2122];xx=fun(k,tt);plot(tt,xx,'r*');tt=0:0.1:22;xx=fun(k,tt);plot(tt,xx,'r');holdoff第二十七頁，共四十三頁，2022年，8月28日holdon和holdoffholdon作用:當前軸及圖形保持而不被刷新，準備接受此后的繪制例如，在作圖時，想在一張圖上同時顯示多組數據(便于觀察等等原因)的時候用holdon(就是等一下,我還想在本圖上畫個東西)holdoff自然就是取消這個功能第二十八頁，共四十三頁，2022年，8月28日指數增長模型的應用及局限性與19世紀以前歐洲一些地區(qū)人口統(tǒng)計數據吻合適用于19世紀后遷往加拿大的歐洲移民后代可用于短期人口增長預測不符合19世紀后多數地區(qū)人口增長規(guī)律不能預測較長期的人口增長過程19世紀后人口數據人口增長率r不是常數(逐漸下降)第二十九頁，共四十三頁，2022年，8月28日美國實際人口與按指數增長模型計算的人口比較年實際人口（百萬）指數增長模型預測人口（百萬）誤差（%）17903.918005.318107.27.31.418209.610.04.2183012.913.76.2184017.118.79.4185023.225.610.3186031.435.010.8187038.647.823.8188050.265.530.5189062.989.642.4190076.0122.561.2191092.0167.682.11920106.5229.3115.3第三十頁，共四十三頁，2022年，8月28日阻滯增長模型(Logistic模型)馬爾薩斯模型為什么不能預測未來的人口呢？這主要是地球上的各種資源只能供一定數量的人生活,隨著人口的增加,自然資源、環(huán)境條件等因素對人口增長的阻滯作用越來越顯著如果當人口較少時,人口的自然增長率可以看作常數的話,那么當人口增加到一定數量以后,這個增長率就要隨人口的增加而減小.因此,應對馬爾薩斯模型中關于凈增長率為常數的假設進行修改.r是x的減函數第三十一頁，共四十三頁，2022年，8月28日[1]

假設最簡單假定（線性函數），r為固有增長率(x很小時).1838年,荷蘭生物數學家韋爾侯斯特(Verhulst)引入常數,用來表示自然環(huán)境條件所能容許的最大人口數（一般說來,一個國家工業(yè)化程度越高,它的生活空間就越大,食物就越多,從而就越大）

第三十二頁，共四十三頁，2022年，8月28日[2]建立模型當時，增長率應為0代入得：

(3)

第三十三頁，共四十三頁，2022年，8月28日將（3）式代入（1）得：模型為：（4）

第三十四頁，共四十三頁，2022年，8月28日[3]

模型求解解方程組(4)得(5)

根據方程（4）作出曲線圖，見下圖，由該圖可看出人口增長率隨人口數的變化規(guī)律.

ox第三十五頁，共四十三頁，2022年，8月28日根據結果（5）作出x~t曲線，見下圖，由該圖可看出人口數隨時間的變化規(guī)律.

S形曲線,x增加先快后慢第三十六頁，共四十三頁，2022年，8月28日也可將方程（4）離散化，得

t=0,1,2,…(6)第三十七頁，共四十三頁，2022年，8月28日[4]

模型的參數估計利用表中數據對r和xm擬合作業(yè)1!第三十八頁，共四十三頁，2022年，8月28日模型檢驗用模型計算2000年美國人口，與實際數據比較實際為281.4(百萬)阻滯增長模型(Logistic模型)第三十九頁，共四十三頁，2022年，8月28日年實際人口阻滯增長模型公式（5）公式（6）預測人口誤差（%）預測人口誤差（%）17903.918005.35.90250.11373.90000.264218107.27.26140.00856.50740.096218209.68.93320.06958.68100.0957183012.910.98990.148111.41530.1151184017.113.52010.209415.12320.1156185023.216.63280.283119.81970.1457186031.420.46210.348326.52280.1553187038.625.17310.347835.45280.0815188050.230.96870.383143.53290.1328189062.938.09860.394356.18840.1067190076.046.86990.383370.14590.0770191092.057.66070.373384.73050.07901920106.570.93590.3339102.46260.03791930123.287.26740.2917118.95090.03451940131.7107