數(shù)學建模講義-微分方程模型

數(shù)學建模講義微分方程模型動態(tài)模型

描述對象特征隨時間(空間)的演變過程

分析對象特征的變化規(guī)律

預報對象特征的未來性態(tài)

研究控制對象特征的手段

根據函數(shù)及其變化率之間的關系確定函數(shù)微分方程建模

根據建模目的和問題分析作出簡化假設

按照內在規(guī)律或用類比法建立微分方程主要內容生物單種群增長模型

3.1

人口增長模型

3.2

傳染病模型生物多種群增長模型

3.3

正規(guī)戰(zhàn)與游擊戰(zhàn)

3.4

捕食系統(tǒng)的Volterra方程

為了保持自然資料的合理開發(fā)與利用，人類必須保持并控制生態(tài)平衡，甚至必須控制人類自身的增長。本節(jié)將建立幾個簡單的單種群增長模型，以簡略分析一下這方面的問題。一般生態(tài)系統(tǒng)的分析可以通過一些簡單模型的復合來研究，大家若有興趣可以根據生態(tài)系統(tǒng)的特征自行建立相應的模型。

美麗的大自然

種群的數(shù)量本應取離散值，但由于種群數(shù)量一般較大，為建立微分方程模型，可將種群數(shù)量看作連續(xù)變量，甚至允許它為可微變量，由此引起的誤差將是十分微小的。離散化為連續(xù)，方便研究3.1

如何預報人口的增長

--Malthus模型與Logistic模型背景

年1625183019301960197419871999人口(億)5102030405060世界人口增長概況中國人口增長概況

年19081933195319641982199019952000人口(億)3.04.76.07.210.311.312.013.0研究人口變化規(guī)律控制人口過快增長指數(shù)增長模型——馬爾薩斯提出(1798)常用的計算公式x(t):時刻t的人口基本假設

:人口(相對)增長率r

是常數(shù),不考慮移民今年人口x0,年增長率rk年后人口隨著時間增加，人口按指數(shù)規(guī)律無限增長模型檢驗

比較歷年的人口統(tǒng)計資料，可發(fā)現(xiàn)人口增長的實際情況與馬爾薩斯模型的預報結果基本相符，例如，1961年世界人口數(shù)為30.6（即3.06×109），人口增長率約為2%，人口數(shù)大約每35年增加一倍。檢查1700年至1961的260年人口實際數(shù)量，發(fā)現(xiàn)兩者幾乎完全一致，且按馬氏模型計算，人口數(shù)量每34.6年增加一倍，兩者也幾乎相同。模型預測

假如人口數(shù)真能保持每34.6年增加一倍，那么人口數(shù)將以幾何級數(shù)的方式增長。例如，到2510年，人口達2×1014個，即使海洋全部變成陸地，每人也只有9.3平方英尺的活動范圍，而到2670年，人口達36×1015個，只好一個人站在另一人的肩上排成二層了。故馬爾薩斯模型是不完善的。幾何級數(shù)的增長Malthus模型實際上只有在群體總數(shù)不太大時才合理，到總數(shù)增大時，生物群體的各成員之間由于有限的生存空間，有限的自然資源及食物等原因，就可能發(fā)生生存競爭等現(xiàn)象。所以Malthus模型假設的人口凈增長率不可能始終保持常數(shù)，它應當與人口數(shù)量有關。指數(shù)增長模型的應用及局限性

與19世紀以前歐洲一些地區(qū)人口統(tǒng)計數(shù)據吻合

適用于19世紀后遷往加拿大的歐洲移民后代

可用于短期人口增長預測

不符合19世紀后多數(shù)地區(qū)人口增長規(guī)律

不能預測較長期的人口增長過程19世紀后人口數(shù)據人口增長率r不是常數(shù)(逐漸下降)阻滯增長模型(Logistic模型)人口增長到一定數(shù)量后，增長率下降的原因：資源、環(huán)境等因素對人口增長的阻滯作用且阻滯作用隨人口數(shù)量增加而變大假設r~固有增長率(x很小時)xm~人口容量（資源、環(huán)境能容納的最大數(shù)量）r是x的減函數(shù)dx/dtx0xmxm/2xmtx0x(t)~S形曲線,x增加先快后慢x0xm/2阻滯增長模型(Logistic模型)參數(shù)估計用指數(shù)增長模型或阻滯增長模型作人口預報，必須先估計模型參數(shù)r或r,xm

利用統(tǒng)計數(shù)據用最小二乘法作擬合例：美國人口數(shù)據（單位~百萬）186018701880……196019701980199031.438.650.2……179.3204.0226.5251.4專家估計阻滯增長模型(Logistic模型)r=0.2557,xm=392.1繼續(xù)最小二乘法設經實際測量已得到n組數(shù)據（xi,yi），i=1,…,n。將數(shù)據畫在平面直角坐標系中，見圖。如果建模者判斷這n個點很象是分布在某條直線附近，令該直線方程為y=ax+b，進而利用數(shù)據來求參數(shù)a和b。由于該直線只是數(shù)據近似滿足的關系式，故yi-(axi+b)=0一般不成立，但我們希望最小此式對a和b的偏導數(shù)均為0，解相應方程組，求得：y=ax+byO(xi,yi)x其中和分別為xi和yi的平均值

如果建模者判斷變量間的關系并非線性關系而是其他類型的函數(shù)，則可作變量替換使之轉化為線性關系或用類似方法擬合。用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)輸出擬合多項式系數(shù)a=[a1,…am,

am+1](數(shù)組)）輸入同長度的數(shù)組X，Y擬合多項式次數(shù)1.lsqcurvefit已知數(shù)據點：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合Matlab的提供了兩個求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得

輸入格式為:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據點選項見無約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）輸入格式為：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin（‘fun’，x0，options）；

3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin（‘fun’，x0，…）；

5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；說明：x=lsqnonlin

（‘fun’，x0，options）；fun是一個事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項有無約束優(yōu)化

例用下面一組數(shù)據擬合中的參數(shù)a，b，k該問題即解最優(yōu)化問題：

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)f=curvefun1(x,tdata)

F(x，tdata)=

，x=(a，b，k)解法1.用命令lsqcurvefit3）運算結果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結論：a=0.0063,b=-0.0034,k=0.2542返回模型檢驗用模型計算2000年美國人口，與實際數(shù)據比較實際為281.4(百萬)模型應用——預報美國2010年的人口加入2000年人口數(shù)據后重新估計模型參數(shù)Logistic模型在經濟領域中的應用(如耐用消費品的售量)阻滯增長模型(Logistic模型)r=0.2490,xm=434.0x(2010)=306.0模型檢驗

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

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

幾乎完全吻合，見圖3.6。

圖3-6Malthus模型和Logistic模型的總結

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

用模擬近似法建立微分方程來研究實際問題時必須對求得的解進行檢驗，看其是否與實際情況相符或基本相符。相符性越好則模擬得越好，否則就得找出不相符的主要原因，對模型進行修改。Malthus模型與Logistic模型雖然都是為了研究種群數(shù)量的增長情況而建立的，但它們也可用來研究其他實際問題，只要這些實際問題的數(shù)學模型有相同的微分方程即可，下面我們來看兩個較為有趣的實例。

年齡分布對于人口預測的重要性

只考慮自然出生與死亡，不計遷移人口發(fā)展方程人口發(fā)展方程一階偏微分方程例2:傳染病模型問題

描述傳染病的傳播過程

分析受感染人數(shù)的變化規(guī)律

預報傳染病高潮到來的時刻

預防傳染病蔓延的手段

按照傳播過程的一般規(guī)律，用機理分析方法建立模型

已感染人數(shù)(病人)i(t)

每個病人每天有效接觸(足以使人致病)人數(shù)為

模型1假設若有效接觸的是病人，則不能使病人數(shù)增加必須區(qū)分已感染者(病人)和未感染者(健康人)建模?模型2區(qū)分已感染者(病人)和未感染者(健康人)假設1）總人數(shù)N不變，病人和健康人的比例分別為2）每個病人每天有效接觸人數(shù)為

,且使接觸的健康人致病建模

~日接觸率SI模型模型21/2tmii010ttm~傳染病高潮到來時刻

(日接觸率)tm

Logistic模型病人可以治愈！?t=tm,di/dt最大模型3傳染病無免疫性——病人治愈成為健康人，健康人可再次被感染增加假設SIS模型3）病人每天治愈的比例為

~日治愈率建模

~日接觸率1/

~感染期

~一個感染期內每個病人的有效接觸人數(shù)，稱為接觸數(shù)。1-1/

idi/dt01>1i0i00ti>11-1/

i0i0t

1di/dt<0接觸數(shù)

=1~閾值如果感染期內有效接觸感染的健康者人數(shù)不超過病人數(shù),患者就會全部治愈。模型4傳染病有免疫性——病人治愈后即移出感染系統(tǒng)，稱移出者SIR模型假設1）總人數(shù)N不變，病人、健康人和移出者的比例分別為2）病人的日接觸率

,日治愈率

,

接觸數(shù)

=/建模需建立的兩個方程模型4SIR模型無法求出的解析解在相平面上研究解的性質模型4消去dtSIR模型相軌線的定義域相軌線11si0D在D內作相軌線的圖形，進行分析si101DSIR模型相軌線及其分析傳染病蔓延傳染病不蔓延s(t)單調減

相軌線的方向P1s0imP1:s0>1/

i(t)先升后降至0P2:s0<1/

i(t)單調降至01/~閾值P3P4P2S0模型4SIR模型預防傳染病蔓延的手段

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

(日治愈率)

醫(yī)療水平傳染病不蔓延的條件——s0<1/

降低s0提高r0

提高閾值1/

降低

(=

/

)

,

群體免疫

的估計模型4SIR模型被傳染人數(shù)的估計記被傳染人數(shù)比例x<<s0i0P1

i0

0,s0

1

小,s0

1提高閾值1/

降低被傳染人數(shù)比例xs0-1/

=

例3：正規(guī)戰(zhàn)與游擊戰(zhàn)戰(zhàn)爭分類：正規(guī)戰(zhàn)爭，游擊戰(zhàn)爭，混合戰(zhàn)爭只考慮雙方兵力多少和戰(zhàn)斗力強弱兵力因戰(zhàn)斗及非戰(zhàn)斗減員而減少，因增援而增加戰(zhàn)斗力與射擊次數(shù)及命中率有關建模思路和方法為用數(shù)學模型討論社會領域的實際問題提供了可借鑒的示例第一次世界大戰(zhàn)Lanchester提出預測戰(zhàn)役結局的模型一般模型

每方戰(zhàn)斗減員率取決于雙方的兵力和戰(zhàn)斗力

每方非戰(zhàn)斗減員率與本方兵力成正比

甲乙雙方的增援率為u(t),v(t)f,g

取決于戰(zhàn)爭類型x(t)~甲方兵力，y(t)~乙方兵力模型假設模型正規(guī)戰(zhàn)爭模型

甲方戰(zhàn)斗減員率只取決于乙方的兵力和戰(zhàn)斗力雙方均以正規(guī)部隊作戰(zhàn)

忽略非戰(zhàn)斗減員

假設沒有增援f(x,y)=

ay,a~乙方每個士兵的殺傷率a=rypy,ry~射擊率，

py~命中率0正規(guī)戰(zhàn)爭模型為判斷戰(zhàn)爭的結局，不求x(t),y(t)而在相平面上討論x與y的關系平方律模型乙方勝游擊戰(zhàn)爭模型雙方都用游擊部隊作戰(zhàn)

甲方戰(zhàn)斗減員率還隨著甲方兵力的增加而增加

忽略非戰(zhàn)斗減員

假設沒有增援f(x,y)=

cxy,c~乙方每個士兵的殺傷率c=rypyry~射擊率py~命中率py=sry/sxsx~甲方活動面積sry~乙方射擊有效面積0游擊戰(zhàn)爭模型線性律模型0混合戰(zhàn)爭模型甲方為游擊部隊，乙方為正規(guī)部隊乙方必須10倍于甲方的兵力設x0=100,rx/ry=1/2,px=0.1,sx=1(km2),sry=1(m2)微分方程的解析解

求微分方程（組）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)

結果：u=tg(t-c)

解

輸入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')結果為:y=3e-2xsin（5x）解輸入命令：

[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t')；

x=simple(x)%將x化簡

y=simple(y)z=simple(z)結果為：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

微分方程的數(shù)值解（一）常微分方程數(shù)值解的定義

在生產和科研中所處理的微分方程往往很復雜且大多得不出一般解。而在實際上對初值問題，一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達式。因此，研究常微分方程的數(shù)值解法是十分必要的。（二）建立數(shù)值解法的一些途徑1、用差商代替導數(shù)

若步長h較小，則有故有公式：此即歐拉法。2、使用數(shù)值積分對方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實際應用時，與歐拉公式結合使用：此即改進的歐拉法。故有公式：3、使用泰勒公式

以此方法為基礎，有龍格-庫塔法、線性多步法等方法。4、數(shù)值公式的精度

當一個數(shù)值公式的截斷誤差可表示為O（hk+1）時（k為正整數(shù)，h為步長），稱它是一個k階公式。k越大，則數(shù)值公式的精度越高。歐拉法是一階公式，改進的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達姆斯外插公式和內插公式。（三）用Matlab軟件求常微分方程的數(shù)值解[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫塔-芬爾格算法ode45：運用組合的4/5階龍格-庫塔-芬爾格算法自變量值函數(shù)值用于設定誤差限(缺省時設定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設定的相對誤差和絕對誤差.

1、在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-文件中的待解方程組應以x的分量形式寫成.

2、使用Matlab軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.注意:解:令y1=x，y2=y1’1、建立m-文件vdp1000.m如下：

functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);

2、取t0=0，tf=3000，輸入命令：

[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3、結果如圖解

1、建立m-文件rigid.m如下：

functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0，tf=12，輸入命令：

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、結果如圖圖中，y1的圖形為實線，y2的圖形為“*”線，y3的圖形為“+”線.穩(wěn)定性問題

在研究許多實際問題時，人們最為關心的也許并非系統(tǒng)與時間有關的變化狀態(tài)，而是系統(tǒng)最終的發(fā)展趨勢。例如，在研究某頻危種群時，雖然我們也想了解它當前或今后的數(shù)量，但我們更為關心的卻是它最終是否會絕滅，用什么辦法可以拯救這一種群，使之免于絕種等等問題。要解決這類問題，需要用到微分方程或微分方程組的穩(wěn)定性理論。在下兩節(jié)，我們將研究幾個與穩(wěn)定性有關的問題。一般的微分方程或微分方程組可以寫成：定義稱微分方程或微分方程組為自治系統(tǒng)或動力系統(tǒng)。（3.28）

若方程或方程組f(x)=0有解Xo，X=Xo顯然滿足（3.28）。稱點Xo為微分方程或微分方程組（3.28)的平衡點或奇點。例7Logistic模型

共有兩個平衡點：N=0和N=K，分別對應微分方程的兩兩個特殊解。前者為No=0時的解而后者為No=K時的解。

當No<K時，積分曲線N=N(t)位于N=K的下方；當No>K時，則位于N=K的上方。從圖3-17中不難看出，若No>0，積分曲線在N軸上的投影曲線（稱為軌線）將趨于K。這說明，平衡點N=0和N=K有著極大的區(qū)別。圖3-17

定義1自治系統(tǒng)的相空間是指以（x1,…,xn）為坐標的空間Rn。

特別，當n=2時，稱相空間為相平面。空間Rn的點集{(x1,…,xn)}|xi=xi(t)滿足(3.28)，i=1,…,n}稱為系統(tǒng)的軌線，所有軌線在相空間的分布圖稱為相圖。定義2設x0是（3.28）的平衡點，稱：

（1）x0是穩(wěn)定的，如果對于任意的ε>0，存在一個δ>0，只要|x(0)-x0|<δ，就有|x(t)-x0|<ε對所有的t都成立。（2）x0是漸近穩(wěn)定的，如果它是穩(wěn)定的且。

微分方程平衡點的穩(wěn)定性除了幾何方法，還可以通過解析方法來討論，所用工具為以下一些定理。（3）x0是不穩(wěn)定的，如果（1）不成立。根據這一定義，Logistic方程的平衡點N=K是穩(wěn)定的且為漸近穩(wěn)定的，而平衡點N=0則是不穩(wěn)定的。解析方法定理1設xo是微分方程的平衡點：若,則xo是漸近穩(wěn)定的若,則xo是漸近不穩(wěn)定的證由泰勒公式，當x與xo充分接近時，有：

由于xo是平衡點，故f(xo)=0。若，則當x<xo時必有f(x)>0，從而x單增；當x>xo時，又有f(x)<0，從而x單減。無論在哪種情況下都有x→xo，故xo是漸進穩(wěn)定的。的情況可類似加以討論。高階微分方程與高階微分方程組平衡點的穩(wěn)定性討論較為復雜，大家有興趣可參閱微分方程定性理論。為了下兩節(jié)的需要，我們簡單介紹一下兩階微分方程組平衡點的穩(wěn)定性判別方法?？疾靸呻A微分方程組：（3.29）

令，作一坐標平移，不妨仍用x記x’，則平衡點xo的穩(wěn)定性討論轉化為原點的穩(wěn)定性討論了。將f(x1,x2)、g(x1,x2)在原點展開，（3.29）又可寫成:考察（3.29）的線性近似方程組:（3.30）其中：記λ1、λ2為A的特征值則λ1、λ2是方程：det（A-λI）=λ2-(a+b)λ+(ad–bc)=0的根令p=a+d,q=ad-bc=|A|，則，記。討論特征值與零點穩(wěn)定的關系（1）若△>0，可能出現(xiàn)以下情形：

①若q>0，λ1λ2>0。當p>0時，零點不穩(wěn)定；當p<0時，零點穩(wěn)定

若q<0，λ1λ2<0

零點為不穩(wěn)定的鞍點③q=0，此時λ1=p，λ2=0，零點不穩(wěn)定。（2）△=0，則λ1=λ2:

當p>0時，零點不穩(wěn)定當p<0時，零點穩(wěn)定（2）△<0，此時若a<0，零點穩(wěn)定,a>0,零點不穩(wěn)定若a=0，有零點為中心的周期解

綜上所述：僅當p<0且q>0時，（3.30）零點才是漸近穩(wěn)定的；當p=0且q>0時（3.30）有周期解，零點是穩(wěn)定的中心（非漸近穩(wěn)定）；在其他情況下，零點均為不穩(wěn)定的。

非線性方程組（3.29）平衡點穩(wěn)定性討論可以證明有下面定理成立:定理2若（3.30）的零點是漸近穩(wěn)定的，則（3.29）的平衡點也是漸近穩(wěn)定的；若（3.30）的零點是不穩(wěn)定的，則（3.29）的平衡點也是不穩(wěn)定的。捕食系統(tǒng)的Volterra方程

問題背景：

意大利生物學家D’Ancona曾致力于魚類種群相互制約關系的研究，在研究過程中他無意中發(fā)現(xiàn)了一些第一次世界大戰(zhàn)期間地中海沿岸港口捕獲的幾種魚類占捕獲總量百分比的資料，從這些資料中他發(fā)現(xiàn)各種軟骨掠肉魚，如鯊魚、鰩魚等我們稱之為捕食者（或食肉魚）的一些不是很理想的魚類占總漁獲量的百分比。在1914~1923年期間，意大利阜姆港收購的魚中食肉魚所占的比例有明顯的增加：年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.810.7

他知道，捕獲的各種魚的比例近似地反映了地中海里各種魚類的比例。戰(zhàn)爭期間捕魚量大幅下降，但捕獲量的下降為什么會導致鯊魚、鰩魚等食肉魚比例的上升，即對捕食者有利而不是對食餌有利呢？他百思不得其解，無法解釋這一現(xiàn)象，就去求教當時著名的意大利數(shù)學家V.Volterra，希望他能建立一個數(shù)學模型研究這一問題。Volterra將魚劃分為兩類。一類為食用魚（食餌），數(shù)量記為x1(t)，另一類為食肉魚（捕食者），數(shù)量記為x2(t)，并建立雙房室系統(tǒng)模型。1、模型建立

大海中有食用魚生存的足夠資源，可假設食用魚獨立生存將按增長率為r1的指數(shù)律增長（Malthus模型），既設：

由于捕食者的存在，食用魚數(shù)量因而減少，設減少的速率與兩者數(shù)量的乘積成正比（競爭項的統(tǒng)計籌算律），即：對于食餌（Prey）系統(tǒng)：λ1反映了捕食者掠取食餌的能力對于捕食者（Predator）系統(tǒng)：捕食者設其離開食餌獨立存在時的死亡率為r2，即：但食餌提供了食物，使生命得以延續(xù)。這一結果也要通過競爭來實現(xiàn)，再次利用統(tǒng)計籌算律，得到：綜合以上分析，建立P-P模型（Volterra方程）的方程組：（3.31）方程組（3.31）反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的相互制約關系。下面我們來分析該方程組。2、模型分析

方程組（3.31）是非線性的，不易直接求解。容易看出，該方程組共有兩個平衡點，即：Po(0,0)是平凡平衡點且明顯是不穩(wěn)定，沒必要研究方程組還有兩組平凡解：

和和所以x1、x2軸是方程組的兩條相軌線。

當x1(0)、x2(0)均不為零時，，應有x1(t)>0且x2(t)>0，相應的相軌線應保持在第一象限中。求（