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

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

描述對象特征隨時間（空間）的演變過程求解方法變量的變化率或導數及變量之間的等式關系式叫微分方程模型微分方程4建立微分方程模型的三種方法

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

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

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

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

交叉路口的寬度

b

典型的車身長度

w

汽車重量μ摩擦系數

x

汽車行駛的距離

T駕駛員的反應時間

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

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

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

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

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

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

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

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

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

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

問題二的模型由題設條件，圓桶受到的阻力應改為，類似問題一的模型，可得到圓桶的速度應滿足如下的微分方程

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

那么經計算可得圓桶的運動時間就不能超過T=13.0025s

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

世界人口增長情況年

1625183019301960197419871999人口（億）5102030405060中國人口增長概況年

19081933195319641982199019952000人口（億）3.04.76.07.210.311.312.013.0

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

建立人口數學模型

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

時刻t的人口數量g

人口自然增長率1.3模型的建立與求解在Δt時間內，人口增長的數量為：令

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

x(t)

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

人口自然增長率,為x(t)的減函數

人口的固有增長率

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

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

時，(3)、當

與

相比很大時，與

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

時，31

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

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

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

,x

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

指數增長模型和阻滯增長模型都是針對人口總數和總的增長率，不涉及年齡結構。事實上，在人口預測中人口按年齡的分布狀況是十分重要的，因為不同年齡人的生育率和死亡率有著很大的差別，兩個國家或地區(qū)目前人口一樣，如果一個國家或地區(qū)年輕人的比例明顯高于另一個國家或地區(qū)，那么二者人口的發(fā)展狀況將大不一樣。在考慮年齡結構的人口模型中，除了時間變量外，年齡是另外一個變量。353.1模型的建立與求解

t時刻小于r歲的人口總數

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

人口年齡密度t時刻單位時間內年齡為r的死亡人數

死亡率函數

令

故

有

令36

t時刻小于r歲的人口總數

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

人口年齡密度

t時刻單位時間內年齡為r的死亡人數

死亡率函數

t時刻單位時間內出生的嬰兒總數

變換得

若已知

有

又

故37：平均一個r歲婦女單位時間生育嬰

兒數：t時刻年齡為r歲的婦女人數與年齡為r歲總人數的比值

：人口年齡密度

：t時刻每個婦女平

均生育胎兒

：生育模式與t無關

：t時刻單位時間內出生的嬰兒總數

令

有

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

是第

個地區(qū)的人口密度函數，

分別表示第個地區(qū)的相對死亡率函數、人口出生率函數和相對擾動密度函數。是第個地區(qū)向第個地區(qū)的移民率。案例分析：美國人口的預報模型年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

認識人口數量的變化規(guī)律，建立人口模型，做出較準確的預報，是有效控制人口增長的前提。利用表1給出的近兩個世紀的美國人口統(tǒng)計數據（以百萬為單位），建立人口預測模型，最后用它預報2010年美國的人口。

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

使用logistic模型預測美國口變化

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

%提出人口數據，提取2,4,6行（人口）數據,并且轉置x=nonzeros(x);%去掉后面的零，并變成列向量t=[1790:10:2000]’;t0=t(1);x0=x(1);%取1970年的人口和年份作為初始數據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))%非線性最小二乘估計得出cs

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

利用向后差分，得到差分方程,其中步長，下面擬合其中的參數和。

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

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

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

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

i0t

1i056由左式可以看出接觸數=1是一個閾值

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

是否大于

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

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

=/γ1）總人數仍為N，病人、健康人和移除者的比例分別為，因此有4.3模型的建立與求解t到t+Δt時刻病人的增加人數為：58t到t+Δt時刻健康人的人數變化為：通過上式無法求出s(t)和i(t)的解析解，我們先做數值計算，設利用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)的數值計算結果圖1s(t),i(t)圖形圖2i-s圖形60消去dt相軌線

的定義域11si0D可得相軌線在數值計算和圖形觀察的基礎上，利用相軌線進行討論

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

相軌線的方向P1:s0>

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

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

提高閾值1/

降低

(=

/

)

,

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

(日治愈率)

醫(yī)療水平

降低s0

由提高r0知可以群體免疫

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

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

得實根

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

從這個領域中的某一點使得則稱點

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

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

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

t時刻漁場的魚量

h(x)

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

自然情況下，單位時

E

單位時間的捕撈率

間魚量的增長量

r

固有增長率

N

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

F(x)單位時間內魚量的變化

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

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

P

魚的單價

T

單位時間的收入

c

單位捕撈率的成本

S

單位時間的支出

R

單位時間的利潤

E

單位捕撈率h(x)

單位時間的捕撈量

3.1.2符號定義3.1.3模型建立及求解在平衡點穩(wěn)定的情況下，有魚量，故可得：利潤最大時，捕撈強度為：3.3計劃捕撈時利潤最大模型76分析：與求捕撈量最大的模型相比，在最大利潤的情況下，捕撈強度和產量都有所減少，而穩(wěn)定的魚量有所增加。將代入得到穩(wěn)定魚量將穩(wěn)定魚量代入自然產量公式可得：77當時，，故在盲目捕撈下，經營者會加大捕撈強度；當時，，故經營者會減小捕撈強度。將

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

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

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

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

n階常系數線性差分方程：2.1對于一階線性常系數差分方程平衡點求解：由解得