第4章_種群數(shù)量的狀態(tài)轉(zhuǎn)移

1、數(shù)學(xué)實驗數(shù)學(xué)實驗n第第4章章 種群數(shù)量的狀態(tài)轉(zhuǎn)移種群數(shù)量的狀態(tài)轉(zhuǎn)移2022-5-3 刻畫世界千變?nèi)f化的規(guī)律，微分方程是最有力刻畫世界千變?nèi)f化的規(guī)律，微分方程是最有力的工具；求解微分方程特征值問題，數(shù)值迭代為最的工具；求解微分方程特征值問題，數(shù)值迭代為最直接的方法；表示數(shù)值迭代結(jié)果與定性分析結(jié)論，直接的方法；表示數(shù)值迭代結(jié)果與定性分析結(jié)論，幾何圖形乃最直觀的方式，它們的結(jié)合在哪里？幾何圖形乃最直觀的方式，它們的結(jié)合在哪里？ -作者作者2022-5-3n狀態(tài)轉(zhuǎn)移與很多因素有關(guān)，如歷史狀態(tài)及時間間隔等系統(tǒng)內(nèi)部的因素、控制變量等外部因素，集中表現(xiàn)為變化率n微分方程是研究變化規(guī)律的有利工具，在科技、工程

2、、經(jīng)濟(jì)管理、生態(tài)、環(huán)境與人口等各個領(lǐng)域有廣泛的應(yīng)用。n如何求解微分方程？高等數(shù)學(xué)介紹了很多求解方法，但是很多實際問題根本求不出解析解，其他方法：數(shù)值解法、數(shù)值解法、圖形方法、微分方程的定性分析方法。圖形方法、微分方程的定性分析方法。n本章將對人口變化、動物種群變遷、網(wǎng)絡(luò)系統(tǒng)的可靠性分析，介紹微分方程（組）的模型建立、數(shù)值解和圖形解等方法，并用MATLAB幾何直觀地展示各種求解方法的求解結(jié)果。引言引言2022-5-34.1 人口問題人口問題n18世紀(jì)末，英國人馬爾薩斯馬爾薩斯（Malthus, 1766-1834）在一本專著中對人口數(shù)量的增長趨勢進(jìn)行了模擬，指出人口的指數(shù)增長模式，提出了人口爆炸

3、問題。n馬爾薩斯人口模型及其改進(jìn)：n單位時間內(nèi)，人口的出生數(shù)量和死亡數(shù)量與人口總單位時間內(nèi)，人口的出生數(shù)量和死亡數(shù)量與人口總數(shù)成比例，即人口出生率和死亡率都是常數(shù)，因此數(shù)成比例，即人口出生率和死亡率都是常數(shù)，因此人口的凈增長率為常數(shù)人口的凈增長率為常數(shù)n設(shè)時刻設(shè)時刻t的人口數(shù)量為的人口數(shù)量為p(t)，人口出生率為，人口出生率為b，死亡率，死亡率為為d，則有：，則有： ()( )( )( )p ttp tbp ttdp tt 2022-5-34.1 人口問題人口問題 從而 其中：k=b-d稱為凈增長率（常數(shù)）。因此馬爾薩斯人口模型如下： 利用高數(shù)知識可得解析解：( )( )( )( )dp tb

4、p tdp tkp tdt00( )( )( )dp tkp tdtp tp0()0( )k t tp tp e2022-5-34.1 人口問題人口問題 利用下表中某國家的人口歷史數(shù)據(jù)對模型進(jìn)行驗證2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法 確定參數(shù)k：因為 ，所以 利用1790年和1840年的數(shù)據(jù)計算得出： 由此預(yù)測，1850年的人口數(shù)量為22898000，誤差為1%，1900年的人口數(shù)量為99476000，誤差為31%，2000年的人口數(shù)量為1877463000，誤差為567%，2050年將達(dá)到80多個億？則該模型對長期預(yù)測不合理。00( )ln()p tk ttp00l

5、n( ( )/)p tpkttln(17069/3929)0.02941840 1790k 2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法 解析解 y=f(x) 數(shù)值解 (xi,yi) 圖形解l歐拉方法l梯形法l改進(jìn)歐拉法l龍格-庫塔法簡單的微分方程復(fù)雜、大型的 微分方程2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法 一階微分方程一階微分方程： 獲取解析解的方法歸類獲取解析解的方法歸類：n分離變量法分離變量法；如dy/dx=x*y, sin(x*y)n齊次方程的變換法齊次方程的變換法；如dy/dx=f(y/x)n線性方程的常數(shù)變異法或公式法線性方程的常數(shù)變異法或公式

6、法。 ( , )dyf x ydx2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法nMATLAB軟件實現(xiàn)軟件實現(xiàn) dsolve(eqn1, eqn2, , c1, , var1, ) 注意注意：y Dy，y D2y 自變量名可以省略，默認(rèn)變量名 t微分方程組初值條件變量組2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n例1 輸入： y=dsolve (Dy=1+y2) y1=dsolve(Dy=1+y2, y(0)=1, x) 輸出：y=tan(t-C1) （通解，一簇曲線） y1=tan(x+1/4*pi) （特解，一條曲線）21,(0)1dyyydx 2022-5

7、-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n例例2 常系數(shù)的二階微分方程 輸入：y=dsolve(D2y-2*Dy-3*y=0, x) y=dsolve(D2y-2*Dy-3*y=0, y(0)=1, Dy(0)=0, x)例例3 非常系數(shù)的二階微分方程 x=dsolve(D2x-(1-x2)*Dx+x=0, x(0)=3,Dx(0)=0)上述兩例的計算結(jié)果怎樣？由此得出什么結(jié)論？ 2 30,(0)1,(0)0yyyyy2( )(1( ) ( )( )0,(0)3,(0)0 x tx tx tx txx2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n例4 非線性微分方程 輸

8、入：x=dsolve(Dx)2+x2=1, x(0)=0) 輸出： x= sin(t) -sin(t) 若欲求解的某個數(shù)值解，如何求解？ t=pi/2; eval(x)22( )( )1,(0)0 x tx tx輸出結(jié)果：ans= -1 12022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法例5 輸入：x,y=dsolve(Dx=3*x+4*y, Dy=-4*x+3*y) x,y=dsolve(Dx=3*x+4*y, Dy=-4*x+3*y, x(0)=0, y(0)=1)34(0)0(0)143dxxyxdtdyyxydt 2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法

9、n數(shù)值解n歐拉法n龍格-庫塔法n研究常微分方程的數(shù)值解法是十分必要的n數(shù)值求解的思想: (變量離散化)n引入自變量點(diǎn)列xnyn 在x0 x1x2xn上求y(xn)的近似值yn 。通常取等步長h，即xn=x0+nh，或xn=xn-1+h，(n=1,2,.)。00( , ),()dyf x yy xydx2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n歐拉方法歐拉方法n在小區(qū)間xn,xn+1上用差商代替微商（近似），n向前歐拉公式： （迭代式） （近似式）特點(diǎn)：f(x,y)取值于區(qū)間xn,xn+1的左左端點(diǎn)1()()nny xy xyh( ( , )yf x y1()()(, ()n

10、nnny xy xhf xy x1(,)nnnnyyhf xy2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n向后歐拉公式 特點(diǎn)：f(x,y)取值于區(qū)間xn,xn+1的右端點(diǎn) 非線性方程，稱“隱式公式” 方法：迭代（y= f(x,y) ）111(,)nnnnyyhf xyyn+1=yn+hf(xn,yn)x=; y=; x(1)=x0; y(1)=y0;for n=1:k x(n+1)=x(n)+n*h; y(n+1)=y(n)+h*f(x(n),y(n); 向前end2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n歐拉方法例1y=- y+x+1, y(0)=1,

11、h=0.1其中f(x,y)=-y+x+1 觀察向前歐拉、向后歐拉計算法計算情況，與精確解進(jìn)行比較，誤差有多大？解：1) 解析解：y=x+e-x y=dsolve(Dy=-y+x+1,y(0)=1,x) y=f(x,y)=-y+x+1; 2)向前歐拉法： yn+1=yn+h(-yn+xn+1) =(1-h)yn+hxn+h2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法3) 向后歐拉法： yn+1=yn+h(-yn+1+xn+1+1) =(yn+hxn+1+h)/(1+h) x1(1)=0; y1(1)=1; y2(1)=1; h=0.1; (died.m) for k=1:10 x

12、1(k+1)=x1(k)+h; y1(k+1)=(1-h)*y1(k)+h*x1(k)+h; y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h); end x1, y1, y2, %(y1-向前歐拉法解，y2-向后歐拉法解) x=0:0.1:1; y=x+exp(-x) %（解析解） plot(x,y,x1,y1,k:,x1,y2,r-)2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法00.20.40.60.8111.11.21.31.42022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法x精確解向前歐拉向后歐拉01110.11.004811.00910.2

13、1.01871.011.02640.31.04081.0291.05130.41.07031.05611.0830.51.10651.09051.12090.61.14881.13141.16450.71.19661.17831.21320.81.24931.23051.26650.91.30661.28741.324111.36791.34871.38552022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n(2)步長h=0.0100.20.40.60.8111.11.21.31.42022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法x精確解向前歐拉向后歐拉01110.11.

14、00481.00441.00530.21.01871.01791.01950.31.04081.03971.04190.41.07031.0691.07170.51.10651.1051.1080.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471.308411.36791.3661.3697結(jié)論結(jié)論：顯然迭代步長h的選取對精度有影響！2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n使用數(shù)值積分 對方程y=f(x,y)，兩邊由xi到xi+1積分，并利用梯形公式，有：111111

15、1()( )( , ( ) ( , ( )(, ()21 ( , ( )(, ()2iixiixiiiiiiiiiiy xy xf t y t dtxxf x y xf xy xhf x y xf xy x11100 ( ,)(,)2()iiiiiihyyf x yf xyyy x即梯形法：2022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n梯形公式n改進(jìn)歐拉公式111 (,)(,)20,1,2,.nnnnnnhyyf xyf xyn1121211()2(,)(,)nnnnnnhyykkkf xykf xyhk(,)nnnyhf xy2022-5-34.2 微分方程的數(shù)值解法微分方程

16、的數(shù)值解法n使用改進(jìn)歐拉公式的例子 以例1為例，用改進(jìn)歐拉公式編程計算，再與精確解比較 yn+1=yn+(h/2)*(-yn+xn+1)+(-yn+1+xn+1+1) =yn+(h/2)*(-yn+xn+1)-(yn+h*(-yn+xn+1) +xn+1+1 =yn+(h/2)*(1-h)*xn+xn+1+2-h+(h-2)*yn Matlab程序： x1(1)=0; y1(1)=1; y2(1)=1; h=0.1; for k=1:10 x1(k+1)=x1(k)+h; y1(k+1)=y1(k)+(h/2)*(1-h)*x1(k)+x1(k+1)+2-h+(h-2)*y1(k); end2

17、022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法x精確解向前歐拉向后歐拉改進(jìn)歐拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.07031.05611.0831.07080.51.10651.09051.12091.10710.61.14881.13141.16451.14940.71.19661.17831.21321.19720.81.24931.23051.26651.250.91.30661.28741.32411.307211.36791.34871.38551

18、.36852022-5-34.2 微分方程的數(shù)值解法微分方程的數(shù)值解法n使用泰勒公式n以此方法為基礎(chǔ)，有龍格-庫塔法、線性多步法等方法n數(shù)值公式的精度n當(dāng)一個數(shù)值公式的截斷誤差可表示為O(hk+1)時（k為正整數(shù)，h為步長），稱它是一個k階公式。k越大，則數(shù)值公式的精度越高；歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆公式和內(nèi)插公式。2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解 t, x=solver(f, ts, x0, options)自變量值函數(shù)值ode23ode45ode113ode15sode23sode23：

19、組合的2/3階龍格-庫塔算法ode45：運(yùn)用組合的4/5階龍格-庫塔算法由待解方程寫出的m-函數(shù)文件ts=t0,tf，t0、tf為自變量的初值和終值函數(shù)的初值用于設(shè)定誤差限（缺省時設(shè)定相對誤差10-3，絕對誤差10-6），命令為：options=odeset(reltol,rt,abstol,at)，rt，at：分別為設(shè)定的相對誤差和絕對誤差。2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n范例 例例1 y=-y+x+1, y(0)=1 標(biāo)準(zhǔn)形式： y=f(x,y) 1) 首先建立M-文件 （weif.m） function f=weif(x,y) f=-y+x+1; 2)

20、求解：x,y=ode23(weif,0,1,1) 3) 作圖形：plot(x,y,r); 4) 與精確解進(jìn)行比較 hold on ezplot(x+exp(-x),0,1)2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n注意注意n在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-函數(shù)文件中的待解方程組應(yīng)以x的分量形式寫成。n使用Matlab軟件求數(shù)值解時，高階微分方程必須等價的變換成一階微分方程組。2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n注意1：n建立M文件函數(shù)function xdot=fun(t,x,y)xdot=f1(t,x(t),y(t)

21、; f2(t,x(t),y(t) ;n數(shù)值計算（執(zhí)行以下命令）t,x,y=ode23(fun,t0,tf,x0,y0)12( )( , ( ), ( )( )( , ( ), ( )dx tf t x ty tdtdy tf t x ty tdtxd(1)=f1(t,x(t),y(t)xd(2)=f2(t,x(t),y(t)xdot=xd; %列向量注意：執(zhí)行命令不能寫在注意：執(zhí)行命令不能寫在M函數(shù)文件中函數(shù)文件中2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n注意2：n例如：y=f(y,y,t) 令 y=x y=xM文件函數(shù)如何寫呢？注意：注意：y(t)是原方程是原方程的解

22、。的解。x(t)只是中間變量只是中間變量( , , )xf x y tyxfunction xdot=fun1(t,x,y) (fun1.m) xdot=f(t,x(t),y(t);x(t)t,x,y=ode23(fun1,t0,tf,x0,y0)如果方程形式是：z=f(t,z,z)？2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n范例n例2 Van der pol 方程2( )(1( ) ) ( )( )0(0)3,(0)0 x tx tx tx txx122212112;(1);(0)3,(0)0;yyyyyyyy令y1=x(t)，y2=x(t)；該方程是否有解？2022

23、-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n編寫M文件（文件名為vdpol.m）： function yp=vdpol(t,y); yp=y(2);(1-y(1)2)*y(2)-y(1);n編寫程序如下：（vdj.m） t,y=ode23(vdpol,0,20,3,0); y1=y(:,1); %原方程的解 y2=y(:,2); plot(t,y1,t,y2,-) pause plot(y1,y2), grid, % 相軌跡圖，即y2(y1)曲線2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解05101520-3-2-10123Van Der Pol Soluti

24、on y1y22022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解-3-2-10123-3-2-10123y1y22022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n范例n例3 考慮Lorenz模型： 其中參數(shù) 解：1）編寫M函數(shù)文件（lorenz.m） 2) 數(shù)值求解并畫三維空間的相平面軌線（ltest.m）112322332123( )( )( )( )( )( )( )( )( )( )( )( )x tx tx t x tx tx tx tx tx t x tx tx t 8,10,2832022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n lo

25、renz.m function xdot=lorenz(t,x) xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x;n ltest.m x0=0 0 0.1; t,x=ode45(lorenz,0,10,x0); plot(t,x(:,1),-,t,x(:,2),*,t,x(:,3),+) pause plot3(x(:,1),x(:,2),x(:,3), grid on2022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解0246810-20-1001020304050 x1x2x32022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解0204060-20020-20020402022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n自變量區(qū)間取0,20，結(jié)果如下：05101520-40-200204060 x1x2x30204060-20020-40-20020402022-5-34.3 Matlab軟件計算數(shù)值解軟件計算數(shù)值解n自變量區(qū)間取0,40，結(jié)果如下：010203040-40-200204060 x1x2x30204060-200