第九講常微分方程模型

第九講常微分方程模型第1頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型分析人口將按指數(shù)規(guī)律無(wú)限增長(zhǎng)！人口將始終保持不變！人口將按指數(shù)規(guī)律減少直至絕滅！人口倍增時(shí)間模型求解第2頁(yè)，共75頁(yè)，2023年，2月20日，星期一Malthus模型預(yù)測(cè)美國(guó)人口第3頁(yè)，共75頁(yè)，2023年，2月20日，星期一Malthus模型預(yù)測(cè)美國(guó)人口第4頁(yè)，共75頁(yè)，2023年，2月20日，星期一Malthus模型預(yù)測(cè)的優(yōu)缺點(diǎn)優(yōu)點(diǎn)短期預(yù)報(bào)比較準(zhǔn)確缺點(diǎn)不適合中長(zhǎng)期預(yù)報(bào)原因預(yù)報(bào)時(shí)假設(shè)人口增長(zhǎng)率r

為常數(shù)。沒(méi)有考慮環(huán)境對(duì)人口增長(zhǎng)的制約作用。第5頁(yè)，共75頁(yè)，2023年，2月20日，星期一2.阻滯增長(zhǎng)模型假設(shè)人口增長(zhǎng)率r(t)是t時(shí)刻人口x(t)的減函數(shù)：其中，xm

為考慮到受自然資源和環(huán)境條件限制所能容納的最大人口數(shù)量（稱最大人口容量）

模型假設(shè)模型建立第6頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型分析（定性分析）人口將遞減并趨向于xm！人口將始終保持xm不變！人口將遞增并趨向于xm！

無(wú)論在哪種情況下，人口最終將趨向于最大人口容量！模型求解第7頁(yè)，共75頁(yè)，2023年，2月20日，星期一人口增長(zhǎng)率達(dá)到最大值第8頁(yè)，共75頁(yè)，2023年，2月20日，星期一阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口第9頁(yè)，共75頁(yè)，2023年，2月20日，星期一阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口第10頁(yè)，共75頁(yè)，2023年，2月20日，星期一阻滯增長(zhǎng)模型預(yù)測(cè)的優(yōu)缺點(diǎn)優(yōu)點(diǎn)中期預(yù)報(bào)比較準(zhǔn)確缺點(diǎn)理論上很好，實(shí)用性不強(qiáng)原因預(yù)報(bào)時(shí)假設(shè)固有人口增長(zhǎng)率r以及最大人口容量xm為定值。實(shí)際上這兩個(gè)參數(shù)（特別是xm

）很難確定，而且會(huì)隨著社會(huì)發(fā)展情況變化而變化。前面圖中曲線末端分叉就是由于這個(gè)原因。第11頁(yè)，共75頁(yè)，2023年，2月20日，星期一6.2藥物在體內(nèi)的分布與排除藥物進(jìn)入機(jī)體形成血藥濃度(單位體積血液的藥物量)血藥濃度需保持在一定范圍內(nèi)——給藥方案設(shè)計(jì)

藥物在體內(nèi)吸收、分布和排除過(guò)程——藥物動(dòng)力學(xué)

建立房室模型——藥物動(dòng)力學(xué)的基本步驟房室——機(jī)體的一部分，藥物在一個(gè)房室內(nèi)均勻分布(血藥濃度為常數(shù))，在房室間按一定規(guī)律轉(zhuǎn)移本節(jié)討論二室模型——中心室(心、肺、腎等)和周邊室(四肢、肌肉等)第12頁(yè)，共75頁(yè)，2023年，2月20日，星期一中心室周邊室給藥排出模型假設(shè)中心室(1)和周邊室(2),容積不變藥物在房室間轉(zhuǎn)移速率及向體外排除

速率，與該室血藥濃度成正比藥物從體外進(jìn)入中心室，在二室間

相互轉(zhuǎn)移,從中心室排出體外模型建立第13頁(yè)，共75頁(yè)，2023年，2月20日，星期一線性常系數(shù)非齊次方程對(duì)應(yīng)齊次方程通解模型建立第14頁(yè)，共75頁(yè)，2023年，2月20日，星期一幾種常見的給藥方式1.快速靜脈注射t=0瞬時(shí)注射劑量d的藥物進(jìn)入中心室,血藥濃度立即為d/V1給藥速率f

(t)和初始條件第15頁(yè)，共75頁(yè)，2023年，2月20日，星期一2.恒速靜脈滴注t>T時(shí),c1(t)和c2(t)按指數(shù)規(guī)律衰減趨于零藥物以恒定速率k進(jìn)入中心室0Tt￡￡第16頁(yè)，共75頁(yè)，2023年，2月20日，星期一吸收室中心室3.口服或肌肉注射相當(dāng)于藥物(劑量d)先進(jìn)入吸收室，吸收后再進(jìn)入中心室吸收室藥量x0(t)第17頁(yè)，共75頁(yè)，2023年，2月20日，星期一參數(shù)估計(jì)各種給藥方式下的c1(t),c2(t)取決于參數(shù)k12,k21,k13,V1,V2以快速靜脈注射為例,在ti(i=1,2,n)測(cè)得c1(ti)由較大的用最小二乘法定A,由較小的用最小二乘法定B,第18頁(yè)，共75頁(yè)，2023年，2月20日，星期一參數(shù)估計(jì)法一進(jìn)入中心室的藥物全部排除第19頁(yè)，共75頁(yè)，2023年，2月20日，星期一參數(shù)估計(jì)法二%構(gòu)造非線性擬合函數(shù)[TWOEXPS.M]functionE=twoexps(a,x,y)x=x(:);y=y(:);Y=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x);E=sum((y-Y).^2)第20頁(yè)，共75頁(yè)，2023年，2月20日，星期一a0=[10011];options=optimset('fminsearch');options.TolX=0.01;options.Display='off';a=fminsearch(@ps,a0,options,x,y)

a=[112.23780.18232.1773]第21頁(yè)，共75頁(yè)，2023年，2月20日，星期一6.3傳染病模型問(wèn)題描述傳染病的傳播過(guò)程分析受感染人數(shù)的變化規(guī)律預(yù)報(bào)傳染病高潮到來(lái)的時(shí)刻預(yù)防傳染病蔓延的手段按照傳播過(guò)程的一般規(guī)律，用機(jī)理分析方法建立模型第22頁(yè)，共75頁(yè)，2023年，2月20日，星期一已感染人數(shù)(病人)i(t)每個(gè)病人每天有效接觸(足以使人致病)人數(shù)為模型1假設(shè)若有效接觸的是病人，則不能使病人數(shù)增加必須區(qū)分已感染者(病人)和未感染者(健康人)建模?第23頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型2區(qū)分已感染者(病人)和未感染者(健康人)假設(shè)1）總?cè)藬?shù)N不變，病人和健康人的比例分別為2）每個(gè)病人每天有效接觸人數(shù)為,且使接觸的健康人致病建模~日接觸率SI模型第24頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型21/2tmii010ttm~傳染病高潮到來(lái)時(shí)刻(日接觸率)tmLogistic模型病人可以治愈！?t=tm,di/dt最大第25頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型3傳染病無(wú)免疫性——病人治愈成為健康人，健康人可再次被感染增加假設(shè)SIS模型3）病人每天治愈的比例為~日治愈率建模~日接觸率1/~感染期~一個(gè)感染期內(nèi)每個(gè)病人的有效接觸人數(shù)，稱為接觸數(shù)。第26頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型3i0i0接觸數(shù)=1~閾值感染期內(nèi)有效接觸感染的健康者人數(shù)不超過(guò)病人數(shù)1-1/i0模型2(SI模型)如何看作模型3(SIS模型)的特例idi/dt01>10ti>11-1/i0t1di/dt<0第27頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型4傳染病有免疫性——病人治愈后即移出感染系統(tǒng)，稱移出者SIR模型假設(shè)1）總?cè)藬?shù)N不變，病人、健康人和移出者的比例分別為2）病人的日接觸率,日治愈率,

接觸數(shù)=/建模需建立的兩個(gè)方程第28頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型4SIR模型無(wú)法求出的解析解在相平面上研究解的性質(zhì)第29頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型4消去dtSIR模型相軌線的定義域相軌線11si0D在D內(nèi)作相軌線的圖形，進(jìn)行分析第30頁(yè)，共75頁(yè)，2023年，2月20日，星期一si101D模型4SIR模型相軌線及其分析傳染病蔓延傳染病不蔓延s(t)單調(diào)減相軌線的方向P1s0imP1:s0>1/

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

i(t)單調(diào)降至01/

~閾值P3P4P2S0第31頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型4SIR模型預(yù)防傳染病蔓延的手段(日接觸率)衛(wèi)生水平(日治愈率)醫(yī)療水平傳染病不蔓延的條件——s0<1/的估計(jì)降低s0提高r0提高閾值1/降低(=/),群體免疫第32頁(yè)，共75頁(yè)，2023年，2月20日，星期一SIR模型被傳染人數(shù)的估計(jì)法一記被傳染人數(shù)比例x<<s0i0P1i00,s01小,s01提高閾值

1/

降低被傳染人數(shù)比例xs0-1/=第33頁(yè)，共75頁(yè)，2023年，2月20日，星期一被傳染人數(shù)的估計(jì)法二X=fzero(‘x-1.2*log(x/0.96)-0.99’,0.5)

X=0.8651第34頁(yè)，共75頁(yè)，2023年，2月20日，星期一6.4多種群生態(tài)數(shù)學(xué)模型意大利生物學(xué)家D’Ancona曾致力于魚類種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚等的比例有明顯增加（見下表），而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭(zhēng)使捕魚量下降，食用魚增加，鯊魚等也隨之增加，但為何鯊魚的比例大幅增加呢？他無(wú)法解釋這個(gè)現(xiàn)象，于是求助于其岳父，著名的意大利數(shù)學(xué)家V.Volterra，希望建立一個(gè)食餌—捕食系統(tǒng)的數(shù)學(xué)模型，定性或定量地回答這個(gè)問(wèn)題.第35頁(yè)，共75頁(yè)，2023年，2月20日，星期一

該模型反映了在沒(méi)有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒(méi)有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡(jiǎn)單的模型.模型（一）不考慮捕獲第36頁(yè)，共75頁(yè)，2023年，2月20日，星期一定理

Volterra微分方程組對(duì)應(yīng)初值問(wèn)題的解是周期函數(shù)，且解的周期平均值為

第37頁(yè)，共75頁(yè)，2023年，2月20日，星期一首先，建立m-文件shier.m如下：functiondx=shier(t,x)dx=zeros(2,1);dx(1)=x(1)*(1-0.1*x(2));dx(2)=x(2)*(-0.5+0.02*x(1));其次，建立主程序shark.m如下：[t,x]=ode45('shier',[015],[252]);plot(t,x(:,1),'-',t,x(:,2),'*')plot(x(:,1),x(:,2))ToMatlab(shark)第38頁(yè)，共75頁(yè)，2023年，2月20日，星期一求解結(jié)果：左圖反映了x1（t）與x2（t）的關(guān)系?？梢圆聹y(cè)：x1（t）與x2（t）都是周期函數(shù)。第39頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型（二）考慮人工捕獲設(shè)表示捕獲能力的系數(shù)為e，相當(dāng)于食餌的自然增長(zhǎng)率由a降為a-e，捕食者的死亡率由c增為c+e設(shè)戰(zhàn)前捕獲能力系數(shù)e=0.3,戰(zhàn)爭(zhēng)中降為e=0.1,則戰(zhàn)前與戰(zhàn)爭(zhēng)中的模型分別為:Volterra原理第40頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型求解:1、分別用m-文件shier1.m和shier2.m定義上述兩個(gè)方程2、建立主程序shark1.m,求解兩個(gè)方程，并畫出兩種情況下鯊魚數(shù)在魚類總數(shù)中所占比例x2(t)/[x1(t)+x2(t)]ToMatlab(shark1)實(shí)線為戰(zhàn)前的鯊魚比例，“*”線為戰(zhàn)爭(zhēng)中的鯊魚比例結(jié)論：戰(zhàn)爭(zhēng)中鯊魚的比例比戰(zhàn)前高！第41頁(yè)，共75頁(yè)，2023年，2月20日，星期一functiony=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;y=diag([r-a*x(2),-d+b*x(1)])*x;shier.mts=0::0.1:35;x0=[25,2];[t,x]=ode45('shier',ts,x0);[t,x],plot(t,x),grid,gtext('x1(t)'),gtext('x2(t)'),pause,plot(x(:,1),x(:,2)),grid,xlabel('x1'),ylabel('x2')shiyan42注：ts中終值(=15)和步長(zhǎng)=(0.1)的確定第42頁(yè)，共75頁(yè)，2023年，2月20日，星期一第43頁(yè)，共75頁(yè)，2023年，2月20日，星期一計(jì)算結(jié)果（數(shù)值，圖形）x(t),y(t)是周期函數(shù)，相圖(x,y)是封閉曲線；觀察，猜測(cè)x(t),y(t)的周期約為10.7;xmax=99.3,xmin=2.0,ymax=28.4,ymin=2.0.用數(shù)值積分可算出x(t)一周期的平均值為25,y(t)一周期的平均值為10.第44頁(yè)，共75頁(yè)，2023年，2月20日，星期一6.5其它生態(tài)數(shù)學(xué)模型存在一大類生態(tài)模型源于對(duì)Volterra模型的改造模型1考慮食餌種群與外界有遷入或遷出

G.R.Walsh（1978）外界有食餌遷入外界有食餌遷出也可以表示人工干預(yù)，如投放或捕獲模型討論食餌－捕食者模型第45頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型2考慮食餌種群內(nèi)部存在生存競(jìng)爭(zhēng)G.Bojadziev表示當(dāng)沒(méi)有捕食者存在時(shí)食餌種群的環(huán)境容納量模型3考慮食餌和捕食者種群內(nèi)部都存在生存競(jìng)爭(zhēng)張錦炎（1979）表示當(dāng)沒(méi)有捕食者存在時(shí)食餌種群的環(huán)境容納量第46頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型4考慮雙方內(nèi)部都存在生存競(jìng)爭(zhēng)，且捕食者另有食物來(lái)源E.C.Pielou平衡點(diǎn)：第47頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型兩種群競(jìng)爭(zhēng)模型競(jìng)爭(zhēng)模型競(jìng)爭(zhēng)排斥原理（CompetitionExclutionLaw）

多個(gè)種群依靠同一個(gè)生存資源而生活，如果生活在同一個(gè)地理空間，

獵取相同食物或營(yíng)養(yǎng)物。在有限的相同生存資源條件下，如果存在競(jìng)爭(zhēng)關(guān)系，

它們必然相互排斥，展開激烈的生存競(jìng)爭(zhēng)。結(jié)局是競(jìng)爭(zhēng)力較弱的種群滅絕，

競(jìng)爭(zhēng)力最強(qiáng)的種群達(dá)到其環(huán)境容納量。

第48頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型討論種群y最終將被滅絕,種群x最終趨于最大容量種群x最終將被滅絕,種群y最終趨于最大容量存在過(guò)正平衡點(diǎn)的一條分界線，將第一象限分成

種群x和種群y的兩個(gè)吸引域。種群x和種群y最終達(dá)到穩(wěn)定的正平衡態(tài)競(jìng)爭(zhēng)排斥原理

是針對(duì)前三種情形得出的結(jié)論，

第四種情況極為罕見。第49頁(yè)，共75頁(yè)，2023年，2月20日，星期一%M函數(shù)functiondy=cwf1(t,y)dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)-0.0008*y(2));dy(2)=0.2*y(2)*(1-0.0012*y(1)-0.001*y(2));第50頁(yè)，共75頁(yè)，2023年，2月20日，星期一

主程序[T,X]=ode45(‘cwf1’,[0200],[200200]);[T,Y]=ode45(‘cwf1’,[0200],[500200]);[T,Z]=ode45(‘cwf1’,[0200],[1200500]);Plot(X(:,1),X(:,2),(Y(:,1),Y(:,2),(Z(:,1),Z(:,2))第51頁(yè)，共75頁(yè)，2023年，2月20日，星期一

第52頁(yè)，共75頁(yè)，2023年，2月20日，星期一模型兩種群互惠模型互惠模型研究多個(gè)種群之間相互依賴、共生現(xiàn)象。模型討論模型有三個(gè)平衡點(diǎn)，分別為第53頁(yè)，共75頁(yè)，2023年，2月20日，星期一兩種群最終達(dá)到穩(wěn)定平衡態(tài)兩種群共生P3為正平衡點(diǎn)P3穩(wěn)定兩種群最終達(dá)不到穩(wěn)定平衡態(tài)P3不穩(wěn)定第54頁(yè)，共75頁(yè)，2023年，2月20日，星期一%m程序functiondy=cwf2(t,y)

dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)+0.0005*y(2));dy(2)=0.2*y(2)*(-1+0.0015*y(1)-0.001*y(2));第55頁(yè)，共75頁(yè)，2023年，2月20日，星期一[T,X]=ode45(‘cwf2’,[0,100],[1600,2800]);[T,Y]=ode45(‘cwf2’,[0,100],[1200,2500]);[T,Z]=ode45(‘cwf2’,[0,100],[2500,2200]);Plot(X(:,1),X(:,2),Y(:,1),Y(:,2),Z(:,1),Z(:,2))Text(2000,2600,2000,2600,'{\sigma_1<1,\sigma_2>1,\sigma_1\sigma_2<1}')第56頁(yè)，共75頁(yè)，2023年，2月20日，星期一第57頁(yè)，共75頁(yè)，2023年，2月20日，星期一6.6常微分方程的數(shù)值解及實(shí)驗(yàn)（一）常微分方程數(shù)值解的定義在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問(wèn)題，一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的。第58頁(yè)，共75頁(yè)，2023年，2月20日，星期一歐拉公式1、用差商代替導(dǎo)數(shù)若步長(zhǎng)h較小，則有故有公式：此即歐拉法（向前歐拉法）。第59頁(yè)，共75頁(yè)，2023年，2月20日，星期一2、使用數(shù)值積分對(duì)方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實(shí)際應(yīng)用時(shí)，與歐拉公式結(jié)合使用：此即改進(jìn)的歐拉法。故有公式：第60頁(yè)，共75頁(yè)，2023年，2月20日，星期一例1求解初值問(wèn)題第61頁(yè)，共75頁(yè)，2023年，2月20日，星期一第62頁(yè)，共75頁(yè)，2023年，2月20日，星期一龍格—庫(kù)特方法

考慮微分中值定理

第63頁(yè)，共75頁(yè)，2023年，2月20日，星期一2階龍格－庫(kù)特公式

第64頁(yè)，共75頁(yè)，2023年，2月20日，星期一利用泰勒展式得到第65頁(yè)，共75頁(yè)，2023年，2月20日，星期一4階龍格－庫(kù)特公式

第66頁(yè)，共75頁(yè)，2023年，2月20日，星期一

第67頁(yè)，共75頁(yè)，2023年，2月20日，星期一一階方程組與高階方程第68頁(yè)，共75頁(yè)，2023年，2月20日，星期一

第69頁(yè)，共75頁(yè)，2023年，2月20日，星期一第70頁(yè)，共75頁(yè)，2023年，2月20日，星期一

Volterra方程解曲線與Volterra方程相軌線第71頁(yè)，共75頁(yè)，2023年，2月20日，星期一3、使用泰勒公式以此方法為基礎(chǔ)，有