数学模型 4.2 常微分方程组模型.pdf

安徽大学剧等科学学院o 4.2 常微分方程组模型：由几个微分方程联立而成的方程组称为微分方程 组。本节建立传染病模型、男生追女生模型和种 群增长模型，并判断方程组的稳定性。4.2.1传染病模型我国目前的法定传染病有甲、乙、丙三类，共39种 0传染病的特点是有病原体，有传染性和流行性，感染后常有免疫性。有些传染病还有季节性或地方 性。传染病的分类尚未统一，有人按病原体分类，有人按传播途径分类。传染病的预防应采取以切断 主要传播环节为主导的综合措施。传染病的传播和0：流行必须具备3个环节，即传染源（能排出病原标的 人或动物）、传播途径（病原体传染他人的途径）及易感者（对该种传染病无免疫力者）。若能完全切断其中的一个环节，即可防止该种传染 病的发生和流行。各种传染病的薄弱环节各不相同 o在预防中应充分利用。除主导环节外对其他环节 也应采取措施，只有这样才能更好地预防各种传染 病。不同类型的传染病，其传播过程有着各自不同的特点，了解这些具体的传染病的传播过程需要了解其病理知识，这里不可能从医学角度一一进行分析，而主要按照一般的传播机理建立几类一般的传 染病模型分析受感染人数的变化规律，讨论终止传 染病蔓延的方法和手段。模型L用表示!时刻的病人数量，假设病人一 旦与健康人群接触就会使健康人群患病，且单位时 间内每个病人能够使健康人患病的人数为2。初始 时刻的病人数为40)=/。利用微元法，考虑依什旬 内病人数的变化，贝|有M%)=疝上式 两边同时除以4,并令4.0,得到dxj dt-x(0)=x0(4.2.1)式(4.2.1)的求解程序及结果如下：dsolveC D x-la mbda*x=0,x(0)=x0,t)a ns=xO exp（la mbda t）即式（4.2.1）的解为x）=（4.2.2）式（4.2.2）为指数形式，故称模型（4.2.1）为指 数增长模型。根据式（4.2.2）,当/一长 时,8,即所有人都会患病。很显然，这与事实不符。模型2.（SI模型）考虑以下假设:（1）不考虑人口的出生、死亡和迁移等种群动力因 素，在疾病传染期所考察地区内总人数始终保持常 数不变；（2）人群分健康人群（易感染者，Susc eptible）和病人（已感染者，Inf ec tive）,t 时刻这两类人在总人数中所占比例分别为s。）和,即有s+近）=1,并设初始时刻患病人数比例为/（0）=i0.(3)病人一旦与健康人群接触，就必然具有一定的 传染性。假设方时刻单位时间内每个病人能传染的 易感染者数量与该地区易感染者比例5(。成正比，比例系数为即，时刻单位时间内每个病人可使 a s个易感染者患病，而病人总数为z()N,因而 t时刻单位时间内共有改N个易感染者患病。利用微元法，内患病人数的变化量为：Ni(t+Ar)-z(O =as(t)i(t)N kt上式两段同时除以,并令加一。，有didt=hsi o结合假设（2）,有di/dt=a z(l-z);z(0)=z0(4.2.3)先观察山/的图像，建M文件myf unl.m。命令如下:f unc tion y=myf unl(i);y=0.01*i*(l-i);再输入命令:f plot C rnyf unT,0,1);f plot函数的命令格式为f plot(f un,lims),表示绘制字符串f un制定的函数在lims=xmin,xma x或 lims=xmin,xma x,ymin,yma x上的图形。但是f un必 须是M文件的函数名或是独立变量为x的字符串。图形如下:图4.2.1 图像由图4.2.1可以看出，当，=0.5时，成/应达到最大值,此时病人数增长得最快，意味着传染病高潮的到 来。模型(4.2.3)的求解程序如下：syms i a lph a t;dsolve(D i-a lph a*i*(l-i),i(0)=i0,t)结果为：a ns=1/(-exp(-a lph a*t)*(T+iO)/i0)这里syms是M a tla b中定义多个变量的。根据运彳亍结 果，有式（4.2.3）的解：止和-e）止阿私（4.2.4）由式（4.2.4）,当 Z=0.5 时=1口（1儿 一 1）/。,这是传染病最高峰时刻，此时才与。成反比。因 为。反映了单位时间每个病人的传染数量，所以也 被称为单位时间接触率，它直接反映了当地的卫生水平可以推迟传染病高潮的到来。函数 的作图命令及图形如下(这里取。=O.5o=O.O l)：f plot c l/(l-exp(-0.5*t)*(-1+0.01)/0.01)0,30);xla bel(t);yla bel(i);图4.2.2患病人数比例 的图像由图4.2.1可以看出,，单调递增，且在，0.5时下 凸，而在，0.5时上凸，且有渐近线，=1,即病人 数量一直递增，在病人比例不到时增长速度非常快,而在病人比例超过时增长速度放缓，最终病人的 比例将接近于100%。显然这与事实不符，究其原因，是因为模型（4.2.3）中没有考虑病人被治愈的情 况。模型3.（SIS模型）增加假设：假设被治愈的人不 具有免疫能力，将成为易感染者，且，时刻单位时 间内治愈者变化率与病人数量成正比，比例系数为4,利用微元法，印+加内患病人数的变化量为Ni+X)ig=as)iQ)NMX,故有dijdt-ods-pi-ai(l-z)-/3 i i=(4.2.5)显然，当尸=0,即没有病人治愈时，模型(4.2.5)变为模型(4.2.3)。令a/尸，则有di/dt(J)(4.2.6)因为。反映了单位时间每个病人的传染数量，而尸 反映了单位时间内每个病人的治愈率，贝1)1/夕可表示 为平均每个病人的治愈时间，/为传染期内平 均每个病人的有效传染数量。图4.2.3(a)时的图像图4.2.3(b)时的图像由图4.2.3(a)和(b)可以看出，当。1时,山/力 随着，的增大呈现先增后减的趋势，且在1-1/。处为 0,即此时，)达到最高峰；而当时，力/力小 于0且随着1的增大一直递减，即海)单调递减且递 减的速度越来越快。模型(4.2.6)的程序为：dsolveCD i=-a lph a*i*(i-(l/sig ma),，i(0)=i0,，t);运行结果为：a ns=(sig ma-1)/(sig ma-exp(-a lph a*(sig ma-1)*t/sig ma)*(-sig ma+l+iO sig ma)/iO/(sig ma-1)*sig ma+exp(-a lph a*(sig ma-1)*t/sig ma)*(-s i g ma+1+i 0 s i g ma)/iO/(sig ma-1)图4.2.4(b)时患病人数比例i的图图4.2.4(a Wwi时患病人数比例1的图像图4.2.4(c)a=Q5，。=2，i0=0.6时患病人数比例，的图形。图4.2.4(a)为 a=0.5，c r=0.2?=0.01 时的图形,图4.2.4(b)为a=0.5,。=2%=0.01 时的图 形，图 4.2.4(c)为 a=0.5,0=2,办=0.611/。时 的图形。可以看出，当。（1时,i从。点出发，单调递减趋于3即此时病人数会越来越小最终趋于0；当01时，近）的单调性取决于。的大小，但是 不论初始值办21-1/。还是io 0,s(0)=s。,r(0)=0；(2)单位时间内每个病人的接触率为a,单位时间 内每个病人的治愈率为夕，(r=a/为传染期内平均 每个病人的有效传染数量；利用微元法，甲+的内移出者的变化量为r(t+A t)-r(t)N=/3 i(t)NA t)等式两边同除以，并令470,则有drjdt=/3 i(4.2.7)故有SIR模型：ds/dt=-ocsi,5(0)=s0；闻应=而_戊7(0)=%;如应=戊40)=石=0.(4.2.9)且S+9)+厂=1.因模型(4.2.9)无法得到解析解,因此考虑作数值分析，先建立M文件：f unc tion y=ill(t,x);a lph a=l,beta=0.3;y=a lph a x(1)*x(2)-beta*x(1),-a lph a*x(l)*x(2);这里假设a=l,尸=。3,然后输入命令:ts=0:50;xO=O.02,0.98;t,x=ode45(ill,ts,xO);plot(t,x(:,1),t,x(:,2),g rid,pa use;plot(x(:,2),x(:,1),g rid得到s)、i以及心)图形如下:0.350.30.250.20.150.10.0500 0.1 0.2 0.3 0.4 0.5 0.6 07 0.8 0.9 1图4.2.5(a)s)、，的图像图4.2.5(b)，(s)的图像下面分析模型(4.2.9)中的s、，和r的变化 情况：(1)由主/d%=-a si知恒有ds/d1(。，即s单调 递减，且对任意时刻有sWs0；(2)力/dt=ais-优可表K为di/dt=ai(s-U(J),(y-alp 若s0l/b,由sM s。知恒有山/a 0,t t。di/dt-ai(s-1/c r)=0,t=t。tG.(4.2.10)由式(4.2.10)知，在，。之前，患病人数在增加，而在之后患病人数在减少，即传染病传播基本得 到控制。(3)由(4.2.9)知，若存在匕使得，(%)=0,而/0,故有网力L=成/山曰二0,进 而有以/说曰=0,即传染病传播完全终止，此时所 有的患病者全部被移出，故传染病传播的最终结果是ii)=O,s(%)+r(%)=l.(4)将(4.2.9)中的第一个方程和第三个方程相 除得到ds/dr=-(js;s(O)=sQ;r(0)=为=。.(4.2.11)求解得到s=“一,令t t,得到2soe-0,该式表明传染病终止后并非所有的人都要患病，然后治愈成为移出者，总有未曾患病的人存在。(5)由 s(G+“G=i得到厂(4)=1 5储)(4.2.12)式(4.2.12)反映了传染病传播终止后，被移出人 数占群体总人数的比率，可作为衡量传染病传染强 度的一种指标，该比值越大传染强度也越大。(6)将式(4.2.9)中的第二个方程和第一个力程 相除得到成/所=1/(5)-1;6)=%;(j=a/3(4.2.13)对式(4.2.13)进行求解，得到即i=I+防 _s+ln(s/so)/b(4.2.14)z(Z)+5(Z)=z0+%+ln(s/so)/(r(4.2.15)令t T t、,且有ii)=O,则s(%)满足%+%一58）+111（汽）/防）/。=0.即:。二（1115（%）-111%）/（霏）_“_10）（4.2.16）而由式(4.2.13)知效i=10)-l(Xsl/b，易感人群此例 随着看的增加而减少，患病人数比例随1的增加而 增加，易感人数比例从l/b递减至患 病人数比例从，。增加至而在之后至传染病结 束时刻匕，易感人数比例继续减少到S/),患病人 数比例也从最多的开始减少，直到传染病终止=。，患病人群完全痊愈。也就是说，设法提高模型的阀值1/。，使得,就可以控制传染病的蔓延。而提高1/。，就是降 低。=。/夕的值，即降低。值或者提高夕值，。反映 了该地区的卫生水平，而反映了该地区的医疗水 平，因此提高卫生水平和医疗水平就能降低提 高夕，就能降低。，从而控制传染病蔓延。SIR案例分析：某高校相对独立，总人口 1万人，某 年冬天该校最初有20人患流感，流行时间持续1个月,累计患病人数大约为3000人。考虑用上述SIR模型 来进行模拟。9980/10000=0.998,s&)=7000/10000=Q7,由式(4.2.16)得至Ug(ns)lnM)/(sG)心幻=(lnQ7-ln0.998)/(Q7-Q998-0.002)=1.1822因为防l/b=0.8459,所以流感会蔓延，且患病人数最 高比例U=力+跖-l/b-ln5o/b=0.0142 o即患病人数从 初始时刻20人不断增加，当易感人群比例从初始时刻的0.998降低到0.8459时，患病人数最多时达到0.0142x10000=142人，之后患病人数开始减少，直至 流感终止。模型5.(SIRP模型)对于天花、乙肝、流感等传染 病来说，注射疫苗可以起到提高免疫力，防止患病 的效果，因此需要对模型4进行改进，增加因为注射 疫苗而具有免疫力的人群(Preventive),即在模 型4的基础上增加假设：易感人群因防疫而减少的变 化率与易感人数比例成正比，比例系数为九/称为 预防系数。依据该假设，可将模型(4.2.9)修改为ds/dt=-asi-ys,s(。)=didt=ods/3 iO)=ia；drl出=阮丫=2.18)将式(4.2.18)中的第二个方程与第一个方程相除，得至Udlds=(l/(j+&/z)I(s+&/z)-1，yoc!0,8 y!a.而由0,sl/(y闻肩=0,s=l/(Jl/a可类似模型4中讨论，且当s=l/。时，患病人数比 例达到最大值，ma x。但是与模型4相比，当s1/(T时,+t Ct T乂。模型5中阀值1/。提高，患病高峰期推迟，患病人数比例最大值会降低。其他定性讨论可类似模型4进行。4.2.2种群增长模型种群生态学的研究起源于人口统计学、渔业资源学 和应用动物学，它以人类、昆虫和动物为主要研究 对象，其理论和方法来源于M.O dum、M.B eg on M ortimer和Prec e,并成为生态学中最为活跃的一个 领域。微分方程是研究种群生态学最重要且最常用 的理论和工具，主要研究两个方面的问题：（1）种 群随时间的演变规律；（2）如何实施人工干扰对种 群进行保护、开发和利用。我们这里主要讨论第一 个问题，考虑随着时间的推移，种群是持续生存还 是走向灭绝？考虑种群的规模是否具有一个或者多 个平衡状态？这种平衡态是静平衡还是动平衡？这 种平衡态是否稳定？在一定的生态环境下有多个物种的生物群体，每一 物种的群体中生物数量的变化既受到本群体自限规 律的制约，同时又受到其他群体的影响。有的群体之间为争夺赖以生存的同一资源和生活空间而相互 竞争，可能会导致竞争能力较弱的群体灭绝，竞争 力较强的群体达到环境容许的最大数量。有的群体 之间相互依存：或两个群体都不能独立生存，但共 处时可相互提供食物，或两个群体都能独立生存，共处时又能相互提供食物，或一个群体可以独立生 存另一个则不能独立生存，共处时可以相互提供食 物等等。还有的群体之间是弱肉强食：一个群体靠丰富的天然资源生长，另一个则靠捕食前一个种群 为生，例如海洋中的食用鱼和鲨鱼等。本节讨论两 种群的相互竞争模型、相互依存模型和食饵一捕食 者模型。模型1.两种群相互竞争模型假设有甲乙两个种群，/时刻的数量分别为尤和加）,固有增长率分别为。和 它们共同生存于同一 系统中，最大容量分别为2和M。对于种群甲而 言，它独自生存时，数量的变化应服从于Log istic 规律，即dx/dt=rxx(1-x!N)(4.2.19)式(4.2.19)中1-M N1反映了种群甲的增长对自身 的阻滞作用，X/V表示这种阻滞作用是线性的，即 相对于甲的最大容量而言，单位数量的种群甲消 耗掉的供养种群甲的食物量为x/N(总量设为1)。当种群乙与种群甲生存于同一系统中时，种群 乙会消耗甲的资源，从而对甲的增长产生进一步的 阻滞作用，且假设相对于种群乙的最大容量生而言乙的存在对种群甲的阻滞作用与种群乙的数量加)成 正比，比例系数为明,即相对于M而言，单位数 量的种群乙消耗掉的供养种群甲的食物量为相对于9而言单位数量甲消耗掉的供养自身的食物量的由 倍。故式(4.2.19)可改写为dx/dt=1x(1-x/Nr-y/N(4.2.20)可以看出，*反映了种群乙占用种群甲的能力，即 种群乙相对于种群甲的竞争能力。若61,则说明 在消耗供养甲的资源时，乙的消耗多于甲自身的消 耗，故而对甲增长的阻滞作用大于甲本身，乙的竞 争力强于甲。同理，得到种群乙的演变规律：dyldt=r2y(-ylN2-(J2xlN(4.2.21)式(4.2.21)中。2反映了种群甲相对于种群乙的竞 争能力。综合(4.2.20)和(4.2.21)得到种群的 相互竞争模型如下：助dt=nx(l-x/NqylN);cfy/dt=ry(l-y/Nex/(4.2.22)由于模型(4.2.22)为非线性模型，无法得到其解 析解，下面通过稳定性及相轨线进行分析。先介绍 二阶微分方程的平衡点与稳定性理论。二阶微分方 程的一般形式可用两个一阶微分方程表示为dxjdt=/(x,y);dy/dt=g(x,y)(4.2.23)代数方程组y（x,y）=0；g（羽 y）=0（4.2.24）的实根（%，%）称为方程（4.2.23）的平衡点,记为 P（x，y。）。若存在某个邻域，使得方程（4.2.23）的解。,y）从这个邻域内的某个点出发满足lim x=x0,lim y=Jo（4.2.25）则称平衡点p。是稳定的（渐近稳定），否则称P。是 不稳定的（不渐近稳定）。用式(4.2.25)判别平衡点的稳定性需要先对微分 方程求出其解析解，但是大多数二阶微分方程都难 以找到其解析解，因此无法根据(4.2.25)来判别,下面考虑一种直接的判别方法。对于二阶线性常系数微分方程dx dt=ax+by;dyj di=cx+dy(4.2.26)(a b令方程右端变量的系数矩阵为人=,且A标逆(C Ct),因此此是(4.2.26)的唯一的平衡点，其稳定性 由特征方程|花T=0的特征根2决定，这里为与A同阶的单位阵。该特征方程又表述为如下形式：牙p=-(a+d);归川(4.2.27)则方程(4.2.26)的一般解的形式为 i”种(4 W%2)或俨如+。24即(4=丸2)，这里Q，Q 为任意常数，4,2=(-P77F)/2为特征方程(4.2.27)的特征根。显然4，%均不为0,否则根据(4.2.27),q=O,这与A可逆矛盾。根据式(4.2.25)知，若，和心均为负数或者均具有负的实 部时，则心为稳定的平衡点；而当41和也有一个为正数或者有正实部，则P。为不稳定的平衡点。实际 上，平衡点又可细分为结点、焦点、鞍点、中心等 类型，其具体的判别方法可依据4，%或P，4，表 4.2.1给出了具体的判别结果。表4.2.1二阶微分方程平衡点类型及稳定性判别准则A Ap,q平衡点类型方程稳定性4 4 0,g 0方-4q 0稳定结点稳定4 o4q 0p 05p:-4g 0不稳定结点不稳定4=办 0,g 0s p:-4g=0稳定退化结点稳定4=%op Qtp:-4q=0不稳定退化结点不稳定九二夕士曲，a 0：g 0 炉-4g 0中心At:=afii,a0p0,p:-4g1（田1（。：i是平衡点e（Nd）稳 定的条件，即在对供养乙的资源竞争中甲的竞争能 力强于乙，因而种群乙将会在竞争中处于劣势，最 后灭绝，同时种群甲会增长到最大容量2。另外，51代表在消耗乙的资源时，甲的消耗要高于乙的,即乙的消耗要小于甲的，因此在同一环境系统中,消耗甲的资源时，乙的消耗也会小于甲的，即有?I。也就是说CT 1是平衡点尸2（。加2）稳定的条件，同时。21，表 明在该竞争系统中，种群甲会灭绝，种群乙达到最 大容量乂。（3）01,生1是平衡点稳定的条件，因为在消耗甲的资源竞争中，乙的竞争力比甲要弱,而在消耗乙的资源竞争中，甲的竞争力也比乙弱,所以竞争双方可以达到共存。（4）舄（0,0）恒为不稳定，因为P1且。2 1，即在 消耗甲的资源时，甲的消耗高于乙，同时在消耗乙 的资源时，乙的消耗高于甲，这种情况长期来讲无 法达到平衡，但是可能有短期的或局部的平衡性。为了更加直观的说明平衡点的稳定性，下面用相轨 线来进行分析，即基于式(4.2.22)分析M t)和y 的变动。在式(4.2.22)中令M人 y)=-R Nfyl N2 M人 y)=l-y1凡-印1 N(4.2.30)下面就上述4个平衡点具体进行分析：(1)%1(6 0,y(?)(Z:x(?)0,y(?)0;Z:x(00,y(?)0 y t图 4.2.6(a)才图(1)yS 4.2.6(c)才图(10 L%1)对于区域2，因为a 0且v0,由式（4.2.22）知f Q）0,y）0 o若相轨线从区域2中的某点出发 移动，随着珀勺增加，。）会逐渐增大,yQ）也会逐渐 增大，即点的横坐标和纵坐标都会变大，相轨线会 向右上方移动。对于区域。2，若相轨线从区域。2中某点出发，随 着，的增加，%）逐渐增大，逐渐减小，即点的 横坐标变大，纵坐标变小，故相轨线会向右下方移 动。对于区域外若相轨线从区域心中某点出发，随着方的增力口，1）和y均逐渐减小，即点的横坐标和 纵坐标都变小，故相轨线会向左下方移动。也就是 说，无论相轨线从哪个区域的任一点出发，当 时都将趋于点4（N,0）o(2)01(%1)o如图4.2.6(b),在相平面。上m=0和v=0无交点，且m=0在v=0下方，这两 条直线将相平面分成三个区域：4:a a：%(oo,y(?)qy）a a：%（oo,y）a 2：%（oo,y（o0,)Q a：f Q 4：f G)0,yG)Q Q：f G)0,yG)0.当看逐渐增加时，若相轨线从区域A的某点出发，则横坐标和纵坐标会变大，即相轨线往右上方移动；若相轨线从区域。2的某点出发，则横坐标变小，纵坐标变大，即相轨线往左上方移动；若相轨线从 区域a的某点出发，则横坐标变大，纵坐标变小，即相轨线往右下方移动；若相轨线从区域2的某点 出发，则横坐标和纵坐标均变小，即相轨线往左下 方移动。这样，最终相轨线将趋于点或因此,这只是一种局部的稳定性，而不是全局的稳定性 O要想巴达到全局稳定，需要附加条件6 1；而若 想鸟达到全局稳定，则需要附加条件巴1。下面考虑两种初值情况下的方程的数值解，仅考虑 第三种情形下的数值解，即令r=0.5,1=Q6,4=2.5,5=L8，N=1.6，M=1。两种初始条件：(1)X0)=X0)=0.1；(2)X0)=lX0)=2o 先画出x)和V的图像，然后画出相轨线。建立M文件：f unc tion dy=xt(t,y);pl=0.5;p2=0.6;rl=2.5;r2=l.8;N1=1.6;N2=1;dy二zeros(2,1);dy(1)=rl*y(1).*(1-y(1)./Nl-pl*y(2)./N2);dy(2)=r2*y(2).*(1-y(2)./N2-p2*y(1)./Nl);输入命令：t0=0;tf=10;t,y=ode45(，xtJ,tO tf,0.1 0.1);%取初始条件x(0)二y(0)=0.1时，求微分 方程数值解subplot(1,2,1);plot(t,y(:,1),t,y(:,2),r);%画出x(t),y(t)曲线图 xla belC t);yla bel(种群数量)；g text(x(t);g text)y(t);%作标记 title(初值分别为0.1,0 1时两种群密度与时间关系);g rid on;t2,y2二ode45(xt,t0 tf,1 2);%取初始条件x(0)=1,y(0)=2时，求微分方 程数值解subplot(1,2,2);plot(t2,y2(:,1),t2,y2(:,2),r);%画出x(t),y(t)曲线图 xla bel C t);yla bel C 种群数量)；g text(x(t);g text(y(t);%作标记 title(初值分别为1,2时两种群密度与时间关系)；图形如下:蛔 翔 整 悬初值分别为1.2时两种群数量与时间关系1.O.图4.2.7种群数量变化图该模型的相轨线画法命令为：tOR；tf=10;t,y=ode45(5 xt,tO tf,0.1 0.1);plot(y(:,1),y(:,2),5 b，);hold on;plot(y(l,1),y(l,2),，r+?);xlabel(，甲种群x);ylabel(乙种群y);%作标记t2,y2=ode45(，xt,tO tf,1 2);hold on;plot(y2(:,1),y2(:,2),b);hold on;plot(y2(1,1),y2(1,2),r*);grid on;titlef甲乙种群相轨线);图形如下:图4.2.8两种群相轨线图形由图4.2.8可以看出，两种群最后区域稳定且数量均 大于0.模型2.两种群相互依存模型 假设种群甲能够独立生存，种群乙无法独立生存，两个种群相互提供食物。因为甲可以独立生存，数 量的变化应服从于Log istic规律，即dx/dt=rxx(1-x/Nj当种群乙与种群甲生存于同一系统中时，种群乙会 种群乙为甲提供食物，有助于甲的增长。假设相对于种群乙的最大容量“而言，乙的存在对种群甲的 促进作用与种群乙的数量y成正比，比例系数为 6，即相对于N?而言，单位数量的种群乙提供给种 群甲的食物量为相对于V而言单位数量甲的供养自 身的食物量的6倍。故上式可改写为dx dt=中c _ N N6N(4.2.31)可以看出，*反映了种群乙对种群甲的供养能力。种群乙不能独自生存，独自生存时会灭亡，设其死 亡率为r2,则种群乙独自生存时满足方程效/或=-与丁 o而甲的存在为乙提供食物，对乙的增长有促进作用，即有dy/dt=4丁(一1+。2%/乂)(4.2.32)同时，乙的增长又会受到自身的阻滞作用，因此式(4.2.32)需要修正为对dt=r2y电+81 N(4.2.33)式(4.2.33)中反映了种群甲对种群乙的供养能 力,且根据式(4.2.32),应有。2田乂1。故有两 种群的相互依存模型：dxfdt=g_xj风+3N);创成二社(一1y/M+0 x/N)(4.2.34)类似于模型（4.2.22）,计算模型（4.2.34）的平衡点，并判断其稳定性，结果如表4.2.3所示。表4.2.3两种群相互依存模型的平衡点及稳定性平衡点Pq稳定性条件祠彳 F(l-6)-仍(If)”仍卬叫（1-城（1巧）+的1）1-印：侦（1-琐。1）1-研。1即0)-仍不稔定由表4.2.3可知,平衡点月(0,0)不稳定；当01就生1 时，平衡点不乂,0)稳定，此时能够独立生存的甲种 群达到最大容量而不能独立生存的乙种群则会 灭亡;0 1,0出1说明种群甲要为种群乙提供足够 的食物才能维持种群乙增长，而此时34/乂+4乂.如图4.2.9,=0和=0将相平面(第一象限)分成四个部分:1:x(r)Q y(0 Q y(t)Q|:直。口:的)0,的)0(4.2.36)当y)存在时，即男生注重学习学业成绩时，会导 致疏远度增长率减小，设单位时间内加)的减少量与y)值成正比，比例系数为(80),故(4.2.36)可修改为dx/dt-(a-by)x=ax-bxy(4.2.37)(2)当男生A发起对女生B的追求后，立即使追求攻 势转化为B对A的好感，设转化系数为仪而随着A发 起对B的追求，A的学业成绩的自然下降率与学业成绩成正比，比例系数为小则有 dy/dt-abxy-ry(4.2.38)%因此得到无外界干扰情况下的男生追女生模型：dxdi=ax-bxy;dy/dt=-ry+cxy(4.2.39)这里C=M。该方程为一非线性自治系统，它是由意 大利数学家V olterra提出，因此也称模型(4.2.39)为V olterra模型。令依-g=0;-zy+oc y=0,得平 衡点P(GQ/b)和。(。，。)o且有人,一 二:j,故(0-br/cyA p=ac/b 0 J,p=0,q=ar0 9 此时为临界状fa 0 态，无法判别尸点的稳定性；A p=_J,p=r-a,q-arQ）的交线在相平面M y上的投影方（%）=E是环绕点尸的闭曲线簇，也就是说学业成绩和疏远度 的指数成周期性变化。为了更进一步说明这个性质,下面画出X和了以及相轨线的图形。先建立M 文件：f unc tion xdot=zh uiqiu(t,x);a=0.5;b=0.1;r=0.3;c=0.05;xdot=(a-b*x(2).*x(1);(-r+b*x(1).*x(2);这里取 a=0.5,/?=0.1/=0.3,a=0.5。再输入命令：ts=0:0.1:50;x0=5,2;t,x=ode45(zh uiqiu,ts,xO);t,x;plot(t,x),g rid,g text(，x(t)，),g text C y(t)，),pa use,plot(x(:,1),x(:,2),g rid得到和图形如下:图4.2.10周期性变化图再画出相轨线，命令如下：tO=O;tf=50;t,y=ode45(?zh uiqiu5,tO tf,5 2);plot(y(:,1),y(:,2),b);h oldonjplot(y(l,1),y(1,2),r+);xla bel(x);yla bel(y);%作标记 g rid on;相轨线如图4.2.n所示Q 111 111_012345678x图4.2.11相轨线二%C由图4.2.10可以看出户和y。）均为周期函数，L两 者有一部分时间的滞后。而相轨线（图4.2.13）则 说明M%）或/（%,y）=k为环绕点P的闭曲线簇。也就是 说，当A的学习成绩y。）下降时，B会疏远A,从而导 致疏远度无上升；于是A就开始奋发图强，学习成 绩丁。）又会上升，故而B看到这一点，就会继续和A来 往，疏远度X下降。与B交往多了，分散了学习时 间，A的学习成绩y又会下降，如此反复。%如图4.2.n,还可以看出，曲线簇中每一个封闭相 轨线略有不同，但是在一周期内X。)和y的平均值元 和y均为常数，且恰为平衡点尸点两个坐标。事实 上，可以将模型(4.2.42)中的第二个方程变形为 dy/ydt=-r+cx 9对该式两边在一个周期时间了内积 分，得到力(dy/ydt)dt=-rT+cxdt(4.2.41)当看经过一个周期T时，点(x,y)绕闭轨线运行一 圈又回到初始点，故也(力/yM d=也力/y=0,由式(4.2.41)得到 tT+c f xdt=0,故十=xdt=中 T c同理可得.=:巾=/O若进一步考虑追求攻势对上述模型的影响，则可设追求攻势与该时刻的疏远度成正比，比例系数为h(),则。反映了追求攻势的作用力。此时模型(4.2.39)可改写为d)c/dt=ax-bxy-hx=(a-h)-bxydy/dt=-ry+cxy-hy=r+h)y+ocy(4.2.42)比较式(4.2.42)与(4.2.39),可以看出，两者 形式完全相同，前者把后者中和y的系数分别换 成了q-/z和厂+/z,因此对式(4.2.42)来说，有 x=(r+h)/c,y=(a-h)/bo也就是说攻势作用力/z的增 大会使下增加，y减小。这说明，强大的爱情攻势有 时不一定能达到满意的效果，反而不利于学业的增 长；有时通过慢慢接触，慢慢了解，再加上适当的 追求行动，女生的疏远度就会慢慢降低，而学习成 绩也不会降低。