(3.8)-第八章 微分方程建模

1第八章微分方程建模1.微分方程簡(jiǎn)介2.物理原理建模3.人口模型4.傳染病模型6.差分方程模型5.平衡點(diǎn)理論及建模7.微分方程的數(shù)值解2第一節(jié)微分方程簡(jiǎn)介3求解微分方程有三種方法：1）求精確解；2）求數(shù)值解（近似解）；3）定性理論方法。動(dòng)態(tài)特性

描述對(duì)象特征隨時(shí)間（空間）的演變過(guò)程求解方法變量的變化率或?qū)?shù)及變量之間的等式關(guān)系式叫微分方程模型微分方程4建立微分方程模型的三種方法

1、微元法——變化量相等

將所討論的變量的變化過(guò)程，分割成許多微小單元，尋找某一個(gè)小區(qū)間的動(dòng)態(tài)平衡，從而得到一個(gè)近似關(guān)系，進(jìn)而建立出其數(shù)學(xué)模型。（1）在空間解析幾何上，可用微元法求曲線的弧長(zhǎng)、平面圖形的面積、旋轉(zhuǎn)曲面的體積、旋轉(zhuǎn)體體積；（2）代數(shù)方面求近似值及流體（3）物理上求變力做功、壓力、靜心距與重心等問(wèn)題如：半球形容器,水從它的底部小孔流出，容器里水面的高度h隨時(shí)間t的變化規(guī)律。

52、根據(jù)數(shù)學(xué)物理規(guī)律建模利用數(shù)學(xué)、力學(xué)、物理、化學(xué)等學(xué)科中的定理或經(jīng)過(guò)實(shí)驗(yàn)檢驗(yàn)的規(guī)律等來(lái)建立微分方程模型。如：動(dòng)力學(xué)中，跳傘運(yùn)動(dòng)員的安全問(wèn)題。原理：牛頓運(yùn)動(dòng)定律以及阻力和速度的關(guān)系。63、模擬近似法在生物、經(jīng)濟(jì)等學(xué)科的實(shí)際問(wèn)題中，許多現(xiàn)象的規(guī)律性不很清楚，即使有所了解也是極其復(fù)雜的，建模時(shí)在不同的假設(shè)下去模擬實(shí)際的現(xiàn)象，建立能近似反映問(wèn)題的微分方程，然后從數(shù)學(xué)上求解或分析所建方程及其解的性質(zhì)，再去同實(shí)際情況對(duì)比，檢驗(yàn)此模型能否刻畫(huà)、模擬某些實(shí)際現(xiàn)象。（例如人口、傳染病模型）

7第二節(jié)物理原理建模8一、紅綠燈問(wèn)題為使正行駛在交叉路口或離交叉路口太近而無(wú)法停下的車輛通過(guò)路口，紅綠燈轉(zhuǎn)換中間還要亮起一段時(shí)間的黃燈，黃燈應(yīng)亮多長(zhǎng)時(shí)間才最為合理。1.1符號(hào)定義法定的行車速度a

交叉路口的寬度

b

典型的車身長(zhǎng)度

w

汽車重量μ摩擦系數(shù)

x

汽車行駛的距離

T駕駛員的反應(yīng)時(shí)間

y黃燈狀態(tài)的時(shí)間91.2模型的建立與求解由牛頓第二定律可得：

即：有初始條件：對(duì)（1）式從0到t進(jìn)行一次積分，并代入初始條件可得：對(duì)（2）式從0到t進(jìn)行一次積分，并代入初始條件可得：（1）（2）10在式②中令，即可得剎車所用時(shí)間為：將③代入得剎車距離為：綜上可得黃燈狀態(tài)的時(shí)間為：即：（3）11二、彈簧振動(dòng)問(wèn)題

設(shè)質(zhì)量為m的質(zhì)點(diǎn)固定在彈簧上沿水平軸在有阻力的介質(zhì)中振動(dòng)，平衡位置是x=0，由胡克定律可知，質(zhì)點(diǎn)受彈力-bx，又因其所受阻力與質(zhì)點(diǎn)運(yùn)動(dòng)速度成正比，故阻力為,且設(shè)質(zhì)點(diǎn)受到外力f(t)，綜上，由牛頓第二定律可得：可將彈簧振動(dòng)問(wèn)題分成以下三種特殊的振動(dòng)模型：1、簡(jiǎn)諧振動(dòng)；2、衰減振動(dòng)；3、受迫振動(dòng)121、簡(jiǎn)諧振動(dòng)（無(wú)外力且無(wú)阻力作用)2、衰減振動(dòng)（無(wú)外力但有阻力作用)3、受迫振動(dòng)（受外力且受阻力作用，設(shè)外力為

)案例分析：放射性廢料的處理美國(guó)原子能委員會(huì)以往處理濃縮的放射性廢料的方法，一直是把它們裝入密封的圓桶里，然后扔到水深為90多米的海底。生態(tài)學(xué)家和科學(xué)家們表示擔(dān)心，怕圓桶下沉到海底時(shí)與海底碰撞而發(fā)生破裂，從而造成核污染。原子能委員會(huì)分辯說(shuō)這是不可能的。為此工程師們進(jìn)行了碰撞實(shí)驗(yàn)。發(fā)現(xiàn)當(dāng)圓桶下沉速度超過(guò)12.2m/s與海底相撞時(shí)，圓桶就可能發(fā)生碰裂。這樣為避免圓桶碰裂，需要計(jì)算一下圓桶沉到海底時(shí)速度是多少?已知圓桶質(zhì)量239.46kg，體積0.2058

m3，海水密度ρ=1035.71kg/m3，若圓桶速度小于12.2m/s就說(shuō)明這種方法是安全可靠的，否則就要禁止使用這種方法來(lái)處理放射性廢料。假設(shè)水的阻力與速度大小成正比例，其正比例常數(shù)k=0.6?，F(xiàn)要求建立合理的數(shù)學(xué)模型，解決如下實(shí)際問(wèn)題（1）判斷這種處理廢料的方法是否合理?（2）一般情況下，v大，k也大；v小，k也小。當(dāng)v很大時(shí)，常用kv來(lái)代替k，那么這時(shí)速度與時(shí)間關(guān)系如何?并求出當(dāng)速度不超過(guò)12.2m/s，圓桶的運(yùn)動(dòng)時(shí)間t和位移s應(yīng)不超過(guò)多少?(k的值仍設(shè)為0.6)1.問(wèn)題一的模型以海平面上的一點(diǎn)為坐標(biāo)原點(diǎn)，垂直向下為坐標(biāo)軸的正向建立坐標(biāo)系。首先要找出圓桶的運(yùn)動(dòng)規(guī)律，由于圓桶在運(yùn)動(dòng)過(guò)程中受到本身的重力以及水的浮力H和水的阻力f的作用，所以根據(jù)牛頓運(yùn)動(dòng)定律得到圓筒受到的合力F滿足（1）又因?yàn)?，，以及可得到圓桶的位移和速度分別滿足下面的微分方程

（2）（3）根據(jù)方程（2），加上初始條件，求得位移函數(shù)為（4）由方程（4），加上初始條件，求得速度函數(shù)為

（4）由，求得圓筒到達(dá)水深90m的海底需要時(shí)間再把它帶入方程（4），求出圓桶到達(dá)海底的速度為顯然此圓桶的速度已超過(guò)，S,可以得出這種處理廢料的方法不合理。因此，美國(guó)原子能委員會(huì)已經(jīng)禁止用這種方法來(lái)處理放射性廢料。計(jì)算的matlab程序如下：clc,clearsymsmVrhogks(t)v(t)%定義符號(hào)常數(shù)和符變量ds=diff(s);%定義s的一階導(dǎo)數(shù)，為了初值條件賦值s=dsolve(m*diff(s,2)-m*g+rho*g*V+k*diff(s),s(0)==0,ds(0)==0);

%使用dsolve解出s的關(guān)系式s=subs(s,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});%常數(shù)賦值s=simplify(s);%化簡(jiǎn)s=vpa(s,6)%顯示小數(shù)形式的位移函數(shù)v=dsolve(m*diff(v)-m*g+rho*g*V+k*v,v(0)==0);

%使用dsolve解出v的關(guān)系式v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});

%常數(shù)賦值v=simplify(v);%化簡(jiǎn)v=vpa(v,6)%顯示小數(shù)形式的速度函數(shù)y=s-90;tt=solve(y);tt=double(tt)%求到達(dá)海底90米處的時(shí)間vv=subs(v,tt);vv=double(vv)%求到底海底90米處的速度結(jié)果如下：

問(wèn)題二的模型由題設(shè)條件，圓桶受到的阻力應(yīng)改為，類似問(wèn)題一的模型，可得到圓桶的速度應(yīng)滿足如下的微分方程

（5）根據(jù)方程（6.34），加上初始條件，求出圓桶的速度，，利用位移，這時(shí)若速度要小于

那么經(jīng)計(jì)算可得圓桶的運(yùn)動(dòng)時(shí)間就不能超過(guò)T=13.0025s

計(jì)算得位移不能超過(guò)84.8439m。通過(guò)這個(gè)模型，也可以得到原來(lái)處理核廢料的方法是不合理的。計(jì)算的matlab程序如下：clc,clearsymsmVrhogkv(t)%定義符號(hào)變量v=dsolve(m*diff(v)-m*g+rho*g*V+k*v^2,v(0)==0);%用dsolve函數(shù)求出v的函數(shù)式v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});%代入數(shù)值v=simplify(v);v=vpa(v,6)%顯示小數(shù)形式的速度函數(shù)T=solve(v-12.2);T=double(T)%求時(shí)間的臨界值Ts=int(v,0,T)%求位移的臨界值計(jì)算結(jié)果如下：結(jié)果分析：由于在實(shí)際中K與V的關(guān)系很難確定，所以上面的模型有它的局限性，而且對(duì)不同的介質(zhì)，比如在水中與在空氣中K與V的關(guān)系也不同。如果假設(shè)K為常數(shù)的話，那么水中的這個(gè)K就比在空氣中對(duì)應(yīng)的V要大一些。在一般情況下，K應(yīng)是V的函數(shù)，即K=K（V），至于是什么樣的函數(shù)，這個(gè)問(wèn)題至今還沒(méi)有解決。22第三節(jié)人口模型23背景與問(wèn)題

世界人口增長(zhǎng)情況年

1625183019301960197419871999人口（億）5102030405060中國(guó)人口增長(zhǎng)概況年

19081933195319641982199019952000人口（億）3.04.76.07.210.311.312.013.0

研究人口變化規(guī)律

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

做出較準(zhǔn)確的預(yù)報(bào)控制人口過(guò)快增長(zhǎng)2417世紀(jì)末，英國(guó)神父Malthus發(fā)現(xiàn)，在人口自然增長(zhǎng)的過(guò)程中（忽略遷入率和遷出率），人口出生率與死亡率幾乎都可以看做常數(shù)，因而兩者之差也幾乎是常數(shù)。這就是說(shuō)，人口增長(zhǎng)率與當(dāng)時(shí)的人口數(shù)量成正比，比例常數(shù)g被稱為人口自然增長(zhǎng)率（它可以通過(guò)人口統(tǒng)計(jì)數(shù)據(jù)得到）。一、Malthus模型（人口增長(zhǎng)率是常數(shù)）1.1模型假設(shè)(1)假設(shè)在社會(huì)穩(wěn)定的前提下，生育和死亡率都比較穩(wěn)定；(2)假設(shè)只考慮本國(guó)內(nèi)部的遷移，而忽略國(guó)際之間的遷移；(3)忽略突發(fā)災(zāi)難性疾病、戰(zhàn)爭(zhēng)等對(duì)人口的影響；(4)假設(shè)國(guó)家對(duì)人口方面的政策基本穩(wěn)定；251.2符號(hào)定義x(t)

時(shí)刻t的人口數(shù)量g

人口自然增長(zhǎng)率1.3模型的建立與求解在Δt時(shí)間內(nèi)，人口增長(zhǎng)的數(shù)量為：令

可得：解得通解為：261.4模型結(jié)果分析比較歷年的人口統(tǒng)計(jì)資料，可發(fā)現(xiàn)人口增長(zhǎng)的實(shí)際情況與Malthus模型的預(yù)報(bào)結(jié)果基本相符。世界人口大約每35年增長(zhǎng)一倍，檢查1700年至1961年的260年中世界人口的實(shí)際數(shù)量，發(fā)現(xiàn)兩者幾乎完全一致。按照該模型計(jì)算，人口數(shù)量每34.6年增長(zhǎng)一倍，兩者幾乎也完全相同。此模型用于短期人口估算有較好的近似程度，但是當(dāng)t→∞時(shí)，有x(t)→∞,可見(jiàn)它不能用于對(duì)人口的長(zhǎng)期預(yù)報(bào)，導(dǎo)致這個(gè)后果的主要原因是在Malthus模型中做了如下假設(shè)：人口自然增長(zhǎng)率g僅與人口出生率和死亡率有關(guān)。且為常數(shù)。這一假設(shè)使得模型簡(jiǎn)化，但也隱含了人口無(wú)限增長(zhǎng)的缺陷，顯然用該模型來(lái)做長(zhǎng)期的人口預(yù)測(cè)會(huì)是不合理的，因此需要對(duì)此進(jìn)行改進(jìn)。27人口增長(zhǎng)到一定數(shù)量后，人口增長(zhǎng)將受到資源、環(huán)境的等因素的阻滯作用，從而導(dǎo)致人口自然增長(zhǎng)率下降，并且這種阻滯作用會(huì)隨著人口數(shù)量的增加而變大，即人口自然增長(zhǎng)率g是人口數(shù)量x(t)的減函數(shù)，從而可建立Logistic模型。二、Logistic模型（人口自然增長(zhǎng)率下降）2.1模型假設(shè)(1)假設(shè)只考慮本國(guó)內(nèi)部的遷移，而忽略國(guó)際之間的遷移；(2)忽略突發(fā)災(zāi)難性疾病、戰(zhàn)爭(zhēng)等對(duì)人口的影響；(3)假設(shè)國(guó)家對(duì)人口方面的政策基本穩(wěn)定。282.2符號(hào)定義2.3模型的建立與求解

x(t)

時(shí)刻t的人口數(shù)量g(x)

人口自然增長(zhǎng)率,為x(t)的減函數(shù)

人口的固有增長(zhǎng)率

自然資源和環(huán)境條件年容納的最大容量根據(jù)假設(shè)條件可知，當(dāng)

時(shí)，，故有：又故可得：解得：29Logistic模型說(shuō)明當(dāng)人口數(shù)量太大時(shí)，種群間會(huì)發(fā)生生存競(jìng)爭(zhēng)，并導(dǎo)致增長(zhǎng)率降低。競(jìng)爭(zhēng)的強(qiáng)弱既和當(dāng)前的種群數(shù)量x有關(guān)，又和環(huán)境的最大容納量有關(guān)。302.4模型結(jié)果分析(1)、(2)、當(dāng)

時(shí)，(3)、當(dāng)

與

相比很大時(shí)，與

相比可以忽略不計(jì)，從而Logistic模型就變成Malthus模型；當(dāng)與相比不是很大時(shí)，與相比就不可以忽略，其作用是使人口的增長(zhǎng)速度減緩下來(lái)。無(wú)論開(kāi)始時(shí)人口處于什么狀態(tài)，隨著時(shí)間的推移，人口總數(shù)最終將趨于其環(huán)境的最大容納量。人口數(shù)量超過(guò)環(huán)境容納量時(shí)，人口數(shù)量將減少；人口數(shù)量小于環(huán)境容納量時(shí)，人口數(shù)量將增加。當(dāng)

時(shí)，31

2.3中最后求得的阻滯增長(zhǎng)模型方程是荷蘭生物數(shù)學(xué)家Verhulst19世紀(jì)中葉提出的。它不僅能夠大體上描述人口及許多物種數(shù)量（如森林中的樹(shù)木、魚(yú)塘中的魚(yú)群等）的變化規(guī)律，而且在社會(huì)經(jīng)濟(jì)鄰域也有廣泛的應(yīng)用，例如耐用消費(fèi)品的售量就可以用它來(lái)描述。此模型還適用于生物種群繁衍，生物生長(zhǎng)，信息傳播，新技術(shù)推廣，傳染病擴(kuò)散，商品銷售，放射現(xiàn)象等?；谶@個(gè)模型能夠描述一些事物的符合邏輯的客觀規(guī)律，人們常稱它為logistic模型。2.5前景與應(yīng)用建模與求解過(guò)程上述兩個(gè)模型的時(shí)間均是連續(xù)的，但在實(shí)際生活中，研究對(duì)象有時(shí)候用離散化的時(shí)間更方便，因此，將微分方程化為式中：即為

；N表示環(huán)境最大容納量。結(jié)合一段時(shí)間的人口數(shù)據(jù)，可用最小二乘法擬合出r、N的數(shù)值。用差分形式表示有

式中，表示第k+1年的人口數(shù)量；表示第k年的人口數(shù)量。2.6差分形式的阻滯增長(zhǎng)模型t

,x

N,x=N是穩(wěn)定平衡點(diǎn)(與r大小無(wú)關(guān))將上式化簡(jiǎn)得令b=1+，=，則可將上式化簡(jiǎn)為一階(非線性)差分方程2.6差分形式的阻滯增長(zhǎng)模型34三、分齡人口模型(分析年齡對(duì)人口數(shù)目的影響，能精確地對(duì)大范圍人口模型進(jìn)行描述）

指數(shù)增長(zhǎng)模型和阻滯增長(zhǎng)模型都是針對(duì)人口總數(shù)和總的增長(zhǎng)率，不涉及年齡結(jié)構(gòu)。事實(shí)上，在人口預(yù)測(cè)中人口按年齡的分布狀況是十分重要的，因?yàn)椴煌挲g人的生育率和死亡率有著很大的差別，兩個(gè)國(guó)家或地區(qū)目前人口一樣，如果一個(gè)國(guó)家或地區(qū)年輕人的比例明顯高于另一個(gè)國(guó)家或地區(qū)，那么二者人口的發(fā)展?fàn)顩r將大不一樣。在考慮年齡結(jié)構(gòu)的人口模型中，除了時(shí)間變量外，年齡是另外一個(gè)變量。353.1模型的建立與求解

t時(shí)刻小于r歲的人口總數(shù)

研究地區(qū)t時(shí)刻的人口總數(shù)

人口年齡密度t時(shí)刻單位時(shí)間內(nèi)年齡為r的死亡人數(shù)

死亡率函數(shù)

令

故

有

令36

t時(shí)刻小于r歲的人口總數(shù)

研究地區(qū)t時(shí)刻的人總數(shù)

人口年齡密度

t時(shí)刻單位時(shí)間內(nèi)年齡為r的死亡人數(shù)

死亡率函數(shù)

t時(shí)刻單位時(shí)間內(nèi)出生的嬰兒總數(shù)

變換得

若已知

有

又

故37：平均一個(gè)r歲婦女單位時(shí)間生育嬰

兒數(shù)：t時(shí)刻年齡為r歲的婦女人數(shù)與年齡為r歲總?cè)藬?shù)的比值

：人口年齡密度

：t時(shí)刻每個(gè)婦女平

均生育胎兒

：生育模式與t無(wú)關(guān)

：t時(shí)刻單位時(shí)間內(nèi)出生的嬰兒總數(shù)

令

有

故38考慮一些不確定因素引起的人口擾動(dòng)以及不同地區(qū)之間人口特性的差異

是第

個(gè)地區(qū)的人口密度函數(shù)，

分別表示第個(gè)地區(qū)的相對(duì)死亡率函數(shù)、人口出生率函數(shù)和相對(duì)擾動(dòng)密度函數(shù)。是第個(gè)地區(qū)向第個(gè)地區(qū)的移民率。案例分析：美國(guó)人口的預(yù)報(bào)模型年1790180018101801830184018501860人口3.95.37.29.612.917.123.231.4年18701880189019001910192019301940人口38.650.262.976.092.0106.5123.2131.7年195019601970198019902000

人口150.7179.3204.0226.5251.4281.4

認(rèn)識(shí)人口數(shù)量的變化規(guī)律，建立人口模型，做出較準(zhǔn)確的預(yù)報(bào)，是有效控制人口增長(zhǎng)的前提。利用表1給出的近兩個(gè)世紀(jì)的美國(guó)人口統(tǒng)計(jì)數(shù)據(jù)（以百萬(wàn)為單位），建立人口預(yù)測(cè)模型，最后用它預(yù)報(bào)2010年美國(guó)的人口。

表1美國(guó)人口統(tǒng)計(jì)數(shù)據(jù)使用malthlus模型預(yù)測(cè)美國(guó)人口變化求解的matlab代碼如下：clc,cleart=[1790:10:2000]';a=textread('data4.txt’);%把原始數(shù)據(jù)保存在純文本文件data4.txt中x=a([2:2:6],:)';%提出人口數(shù)據(jù)x=nonzeros(x);%去0y=log(x);%令y=lnxp=polyfit(t,y,1)%擬合lnx=a*t+bg=p(1)%p（1）代表g：人口自然增長(zhǎng)率x0=exp(p(2))%因?yàn)閑^（a*t）*e^b=x,而x=x0*e^g*t,所以x0=e^bY=polyval(p,t);%y=polyval(p,t)為返回對(duì)應(yīng)自變量t在給定系數(shù)P的多項(xiàng)式的值。此時(shí)值為lnxX=exp(Y);%還原成預(yù)測(cè)的xplot(t,X,'*-',t,x,'+-')title('malthus預(yù)測(cè)美國(guó)人口變化')xlabel('年份')ylabel('人口（百萬(wàn)）')結(jié)果如下所示：結(jié)果分析由擬合的曲線可以知道m(xù)althus模型適用于短期內(nèi)的人口預(yù)測(cè)，并且模型是指數(shù)性的，與實(shí)際情況不符2.參數(shù)估計(jì)

使用logistic模型預(yù)測(cè)美國(guó)口變化

求解的matlab代碼如下：clc,cleara=textread('data4.txt’);%把原始數(shù)據(jù)保存在純文本文件data4.txt中，讀取data4里保存的數(shù)據(jù)，6*8矩陣x=a([2:2:6],:)’;

%提出人口數(shù)據(jù)，提取2,4,6行（人口）數(shù)據(jù),并且轉(zhuǎn)置x=nonzeros(x);%去掉后面的零，并變成列向量t=[1790:10:2000]’;t0=t(1);x0=x(1);%取1970年的人口和年份作為初始數(shù)據(jù)fun=@(cs,td)cs(1)./(1+(cs(1)/x0-1)*exp(-cs(2)*(td-t0)));%cs(1)=xm,cs(2)=rcs=lsqcurvefit(fun,rand(2,1),t(2:end),x(2:end),zeros(2,1))%非線性最小二乘估計(jì)得出cs

xhat=fun(cs,[t;2010])%預(yù)測(cè)已知年代和2010年的人口t1=[t;2010];plot(t1,xhat,'*-',t,x,'+-’)title('logistic預(yù)測(cè)美國(guó)人口變化’)xlabel('時(shí)間t’)ylabel('人口（百萬(wàn)）')3.結(jié)果顯示4.結(jié)果分析最后求得xm=342.4395，r=0.0274,2010年人口的預(yù)測(cè)值為282.6788百萬(wàn)人擬合的曲線也與理論的logistic曲線向吻合模型建立：把Logistic方程表示為

利用向后差分，得到差分方程,其中步長(zhǎng)，下面擬合其中的參數(shù)和。

求解的matlab代碼如下：clc,cleara=textread('data4.txt’);%把原始數(shù)據(jù)保存在純文本文件data4.txt中x=a([2:2:6],:)';x=nonzeros(x);%提取人口數(shù)據(jù)，去零t=[1790:10:2000]';%時(shí)間數(shù)據(jù)a=[ones(21,1),-x(2:end)];%構(gòu)造a矩陣作差分方程右端b=diff(x)./x(2:end)/10;%表示差分方程左端cs=a\b;%其實(shí)就是a*cs=b，得出cs兩個(gè)值r=cs(1),xm=r/cs(2)%cs(1)代表差分方程右端系數(shù)r；cs（2）代表差分方程右端s，因?yàn)閟=r/xm,所以xm=r/s結(jié)果顯示結(jié)果分析r=0.0247，Xm=373.5135上述方法與logistic模型擬合的參數(shù)相比，由于擬合方法不同，所得出的參數(shù)也不相同49第四節(jié)傳染病模型50傳染病模型背景與問(wèn)題基本方法·傳染病具有極大危害（如艾滋病、SARS等）·描述傳染病的傳播過(guò)程·分析受感染人數(shù)的變化規(guī)律·預(yù)報(bào)傳染病高潮到來(lái)的時(shí)刻·預(yù)防傳染病蔓延的手段不是從醫(yī)學(xué)角度分析各種傳染病的特殊機(jī)理，而是按照傳播過(guò)程的一般規(guī)律建立數(shù)學(xué)模型。51一、指數(shù)增長(zhǎng)模型1.1模型假設(shè)1.2符號(hào)定義i(t)時(shí)刻t已經(jīng)感染人數(shù)每個(gè)病人每天有效接觸（足以使人致?。┤藬?shù)1.3模型的建立與求解(1)、每個(gè)病人每天有效接觸（足以使人致?。┤藬?shù)為；(2)、一個(gè)人得病后，久治不愈，并在傳染期內(nèi)不會(huì)死亡。t到t+Δt時(shí)刻病人的增加人數(shù)為：設(shè)t=0時(shí)，有個(gè)病人，可得：其解為：隨著t的增加，病人人數(shù)會(huì)按指數(shù)函數(shù)無(wú)限增長(zhǎng)，之顯然不符合實(shí)際情況，故其只適用于傳染病傳播初期52二、SI模型(區(qū)分病人和健康人，病人不會(huì)再被感染）2.1模型假設(shè)(1)每個(gè)病人每天有效接觸人數(shù)為,且使接觸的健康人致病。(2)一個(gè)人得病后，久治不愈，并在傳染期內(nèi)不會(huì)死亡；(3)在疾病傳播期間內(nèi)所考察地區(qū)的總?cè)藬?shù)不變。2.2符號(hào)定義2.3模型的建立與求解i(t)時(shí)刻t已經(jīng)感染人數(shù)占總?cè)藬?shù)比例S(t)時(shí)刻t未被感染人數(shù)占總?cè)藬?shù)比例每個(gè)病人每天有效接觸（足以使人致?。┤藬?shù)N該考察地區(qū)總?cè)藬?shù)t到t+Δt時(shí)刻病人的增加人數(shù)為：53用分離變量法得其解：令

可得：ii010t1/2tmt=tm時(shí)，di/dt最大tm傳染病高潮到來(lái)時(shí)刻54三、SIS模型（病人可治愈成為健康人，健康人可再次被感染）3.1模型假設(shè)

在SI模型的基礎(chǔ)上，增加假設(shè)每天可治愈病人μ人3.2符號(hào)定義i(t)時(shí)刻t已經(jīng)感染人數(shù)占總?cè)藬?shù)比例S(t)時(shí)刻t未被感染人數(shù)占總?cè)藬?shù)比例每個(gè)病人每天有效接觸（足以使人致?。┤藬?shù)N該考察地區(qū)總?cè)藬?shù)1/μ平均傳染期3.3模型的建立與求解t到t+Δt時(shí)刻病人的增加人數(shù)為：55令可將上式化為：

~一個(gè)感染期內(nèi)每個(gè)病人的有效接觸人數(shù)，稱為接觸數(shù)。將上式改寫為idi/dt0>10ti>11-1/

i0t

1i056由左式可以看出接觸數(shù)=1是一個(gè)閾值

時(shí)i(t)最終趨于零，也就說(shuō)疾病的傳播被完全控制了如果，的增減性取決于

是否大于

。這時(shí)它將有一個(gè)非零的極限值，而不可能將疾病完全控制住。57四、SIR模型(人患病痊愈后有長(zhǎng)期免疫力)4.1模型假設(shè)假設(shè)患過(guò)傳染病而完全痊愈的任何人都具有長(zhǎng)期免疫力，不考慮反復(fù)受傳染的情形，且傳染病潛伏期可忽略不計(jì)，即一個(gè)人患病后立即成為傳染者。4.2符號(hào)定義

2）病人的日接觸率，日移除率γ,記

=/γ1）總?cè)藬?shù)仍為N，病人、健康人和移除者的比例分別為，因此有4.3模型的建立與求解t到t+Δt時(shí)刻病人的增加人數(shù)為：58t到t+Δt時(shí)刻健康人的人數(shù)變化為：通過(guò)上式無(wú)法求出s(t)和i(t)的解析解，我們先做數(shù)值計(jì)算，設(shè)利用MATLAB軟件編程可得：故可得：（通常約為0）59t012345678i(t)0.02000.03900.07320.12850.20330.27950.33120.34440.3247s(t)0.98000.95250.90190.81690.69270.54380.39950.28390.2027t91015202530354045i(t)0.28630.24180.07870.02230.00610.00170.00050.00010s(t)0.14930.11450.05430.04340.04080.04010.03990.03990.0398表1s(t),i(t)的數(shù)值計(jì)算結(jié)果圖1s(t),i(t)圖形圖2i-s圖形60消去dt相軌線

的定義域11si0D可得相軌線在數(shù)值計(jì)算和圖形觀察的基礎(chǔ)上，利用相軌線進(jìn)行討論

為閾值61相軌線i(s)及其分析s(t)單調(diào)減

相軌線的方向P1:s0>

i(t)先升后降至0傳染病蔓延P2:s0<

i(t)單調(diào)降至0傳染病不蔓延滿足62預(yù)防傳染病蔓延的手段傳染病不蔓延的條件——s0<1/

提高閾值1/

降低

(=

/

)

,

(日接觸率)衛(wèi)生水平

(日治愈率)

醫(yī)療水平

降低s0

由提高r0知可以群體免疫

的估計(jì)忽略可得五、SEIR模型(人被感染傳染病后要經(jīng)歷病毒潛伏期)假設(shè)患過(guò)傳染病而完全痊愈的任何人都具有長(zhǎng)期免疫力，不考慮反復(fù)受傳染的情形．并設(shè)傳染病的潛伏期不可忽略不計(jì)，即一個(gè)人患了病之后需要經(jīng)過(guò)一段時(shí)間后才能成為傳染者。在所考慮時(shí)期內(nèi)人口總數(shù)保持固定水平N不變，即不考慮出生及其他原因引起的死亡，以及遷出、遷入等情況．記t時(shí)刻已經(jīng)感染人數(shù)(病人)占總?cè)藬?shù)的比例為i(t)(第一類人)；t時(shí)刻未被感染人數(shù)(健康人)占總?cè)藬?shù)的比例為s(t)(第二類人)；t時(shí)刻第三類人數(shù)(潛伏者)為e(t)；t時(shí)刻第四類人數(shù)(移除者)占總?cè)藬?shù)的比例為r(t)；因此有

ω1表示疾病的潛伏周期，即潛伏者經(jīng)過(guò)時(shí)間ω1后后可能轉(zhuǎn)化為病人；ω2表示治愈者在經(jīng)過(guò)周期ω2后喪失免疫力，進(jìn)而轉(zhuǎn)化為易感者；λ表示個(gè)體與疾病的有效接觸率；α1表示個(gè)體從潛伏狀態(tài)轉(zhuǎn)化為發(fā)病狀態(tài)的比例；α2表示染病者轉(zhuǎn)化為治愈者的比例；b表示個(gè)體的直接免疫率；u表示個(gè)體從潛伏狀態(tài)轉(zhuǎn)化為免疫狀態(tài)的比例．5.1模型假設(shè)5.2符號(hào)定義5.3模型的建立傳染病SEIR模型的傳播機(jī)理為：當(dāng)易感者與染病者接觸并被傳染后即成為染病者，染病者恢復(fù)后就進(jìn)入移除者群體．當(dāng)移除者不具有永久免疫力時(shí)，經(jīng)過(guò)一段時(shí)間后又會(huì)成為易感者．有些疾病在潛伏期也具有傳染力，而且將接觸率推廣到一般形式的接觸率，具有一般形式的接觸率模型是對(duì)具體模型的一種抽象，對(duì)于有些傳染病在流行期間，易感者一旦被感染上病毒，在未發(fā)病之前(即潛伏期)就對(duì)外具有傳染性．根據(jù)上述描述的傳染病傳播機(jī)理，建立的SEIR微分方程模型如下66第五節(jié)平衡點(diǎn)理論及建模67對(duì)于某些實(shí)際問(wèn)題，建模的主要目的不是得到每個(gè)瞬間的動(dòng)態(tài)，而是更關(guān)注于系統(tǒng)在某種假設(shè)下穩(wěn)定狀態(tài)的特征，特別是當(dāng)時(shí)間趨于無(wú)窮時(shí)，系統(tǒng)的趨勢(shì)，此時(shí)，可利用穩(wěn)定性理論，直接研究系統(tǒng)平衡狀態(tài)的穩(wěn)定性。一、一階微分方程的平衡點(diǎn)及穩(wěn)定性1.1平衡點(diǎn)設(shè)一階微分自治方程:解方程

得實(shí)根

，稱其為原微分方程的平衡點(diǎn)1.2判斷平衡點(diǎn)是否穩(wěn)定1.2.1間接法若存在鄰域，使得x(t)從這個(gè)領(lǐng)域中的某個(gè)x(0)出發(fā)，使得，則稱平衡點(diǎn)是穩(wěn)定的。681.2.2直接法若，則對(duì)于原微分方程是穩(wěn)定的；若，則對(duì)于原微分方程是不穩(wěn)定的。二、二階微分方程的平衡點(diǎn)及穩(wěn)定性2.1平衡點(diǎn)對(duì)于二階自治微分方程:令:解得實(shí)根,,稱其為原微分方程的平衡點(diǎn)，記作692.2判斷平衡點(diǎn)是否穩(wěn)定2.2.1間接法若存在某個(gè)鄰域，使得

從這個(gè)領(lǐng)域中的某一點(diǎn)使得則稱點(diǎn)

是穩(wěn)定的；否則，則不是穩(wěn)定的。2.2.2直接法2.2.2.1當(dāng)為線性常系數(shù)方程

對(duì)線性常系數(shù)方程其系數(shù)矩陣70令若，則平衡點(diǎn)是穩(wěn)定的；若，則平衡點(diǎn)是不穩(wěn)定的。當(dāng)為一般的非線性線性方程

利用泰勒公式，得到原方程的近似線性方程其系數(shù)矩陣再判斷p，q的正負(fù)性，結(jié)論同上71三、二階衡點(diǎn)及穩(wěn)定性可再生資源的管理模型——以捕魚(yú)模型為例資源分為可再生資源（林業(yè)、漁業(yè)等）和非可再生資源（礦業(yè)等）。再生資源應(yīng)適度開(kāi)發(fā)——在持續(xù)穩(wěn)產(chǎn)前提下實(shí)現(xiàn)最大產(chǎn)量或最佳效益。背景723.1基本捕魚(yú)模型x(t)

t時(shí)刻漁場(chǎng)的魚(yú)量

h(x)

單位時(shí)間的捕撈量f(x)

自然情況下，單位時(shí)

E

單位時(shí)間的捕撈率

間魚(yú)量的增長(zhǎng)量

r

固有增長(zhǎng)率

N

環(huán)境允許的最大魚(yú)量

F(x)單位時(shí)間內(nèi)魚(yú)量的變化

3.1.1符號(hào)定義3.1.2模型建立及求解自然情況下，魚(yú)量的增長(zhǎng)滿足Logistic模型規(guī)律，故可得：又有：由于并不需要求解魚(yú)量的動(dòng)態(tài)變化過(guò)程，而只希望知道穩(wěn)定的魚(yú)量和保持穩(wěn)定魚(yú)量的條件，并進(jìn)一步確定最大可持續(xù)產(chǎn)量，所以我們可以利用微分方程的平衡點(diǎn)及穩(wěn)定性理論求解此類問(wèn)題。73單位時(shí)間內(nèi)魚(yú)量的變化＝單位時(shí)間內(nèi)自然情況下魚(yú)量的增長(zhǎng)量－單位時(shí)間內(nèi)捕撈量，故：根據(jù)一階微分方程的平衡點(diǎn)及穩(wěn)定性理論，令：得到兩個(gè)平衡點(diǎn)：易求得：故可知：當(dāng)，是穩(wěn)定的，不穩(wěn)定；

當(dāng)，是穩(wěn)定的，不穩(wěn)定。結(jié)論：當(dāng)捕撈適度（E<r)時(shí)，魚(yú)群總量最后將穩(wěn)定在；當(dāng)過(guò)度捕撈(E>r)時(shí)，魚(yú)群總量最后將趨近于0，即魚(yú)群將滅亡。743.2魚(yú)量穩(wěn)定且捕撈量最大的模型求解以漁場(chǎng)魚(yú)量為橫坐標(biāo)，自然情況下的增長(zhǎng)量或捕撈量為縱坐標(biāo)，將和做在同一個(gè)坐標(biāo)軸上，如右圖所示y=rxhPx0hmx0*=N/2P*y=E*xy0y=h(x)=ExxNy=f(x)保證魚(yú)量穩(wěn)定:找交點(diǎn);捕撈量最大：找最大縱坐標(biāo)值的交點(diǎn)。分析：魚(yú)量穩(wěn)定且捕撈適度，故有E<r，又y=f(x)在原點(diǎn)的切線斜率為r，故兩條線一定會(huì)有交點(diǎn)，交點(diǎn)表示自然情況下的增長(zhǎng)量和捕撈量相等，交點(diǎn)的橫坐標(biāo)為魚(yú)量的穩(wěn)定值，縱坐標(biāo)為此時(shí)的捕撈量。平衡點(diǎn)：最大捕撈量：捕撈率：75

P

魚(yú)的單價(jià)

T

單位時(shí)間的收入

c

單位捕撈率的成本

S

單位時(shí)間的支出

R

單位時(shí)間的利潤(rùn)

E

單位捕撈率h(x)

單位時(shí)間的捕撈量

3.1.2符號(hào)定義3.1.3模型建立及求解在平衡點(diǎn)穩(wěn)定的情況下，有魚(yú)量，故可得：利潤(rùn)最大時(shí)，捕撈強(qiáng)度為：3.3計(jì)劃捕撈時(shí)利潤(rùn)最大模型76分析：與求捕撈量最大的模型相比，在最大利潤(rùn)的情況下，捕撈強(qiáng)度和產(chǎn)量都有所減少，而穩(wěn)定的魚(yú)量有所增加。將代入得到穩(wěn)定魚(yú)量將穩(wěn)定魚(yú)量代入自然產(chǎn)量公式可得：77當(dāng)時(shí)，，故在盲目捕撈下，經(jīng)營(yíng)者會(huì)加大捕撈強(qiáng)度；當(dāng)時(shí)，，故經(jīng)營(yíng)者會(huì)減小捕撈強(qiáng)度。將

代入可得盲目捕撈下漁場(chǎng)穩(wěn)定魚(yú)量，由成本和價(jià)格比決定.3.4捕撈過(guò)度時(shí)模型求解令可以得其解：臨界捕撈強(qiáng)度S(E)T(E)0rEpNEEsS(E)，T(E)或求解的matlab代碼如下：clc,clearsymsrNExpc%r:固有增長(zhǎng)率N:環(huán)境允許的最大魚(yú)量E:捕撈強(qiáng)度

x:平衡狀態(tài)魚(yú)量p:魚(yú)的單價(jià)c:單位捕撈率的成本R=p*N*E*(1-E/r)-c*E;%單位時(shí)間的利潤(rùn)公式，已知平衡狀態(tài)魚(yú)量x=N*（1-E/r）ER=simplify(solve(diff(R,E),E))

%diff函數(shù)表示R對(duì)E的導(dǎo)數(shù)公式，然后求解出利潤(rùn)最大時(shí)捕撈強(qiáng)度ErxR=simplify(subs(N*(1-E/r),{E},{ER}))%求解出最大捕撈強(qiáng)度時(shí)平衡狀態(tài)魚(yú)量hR=simplify(subs(r*x*(1-x/N),{x},{xR}))%代入捕撈量公式求出hREs=simplify(solve(R,E));%求出過(guò)度捕撈時(shí)的捕撈強(qiáng)度Es=Es(2)%取非零的捕撈強(qiáng)度xs=simplify(subs(N*(1-E/r),{E},{Es}))%求解出過(guò)度捕撈時(shí)穩(wěn)定魚(yú)量

結(jié)果顯示80第六節(jié)差分方程模型81一、差分方程簡(jiǎn)介1.1差分定義對(duì)函數(shù)，規(guī)定t只取非負(fù)整數(shù)。記為變量y在t點(diǎn)的取值。一階差分：二階差分：n階差分：1.2差分方程對(duì)含有未知函數(shù)的差分的函數(shù)方程，稱其為差分方程如：或滿足差分方程的序列稱為該方程的解82

n階常系數(shù)線性差分方程：2.1對(duì)于一階線性常系數(shù)差分方程平衡點(diǎn)求解：由解得