食饵—捕食者模型稳定性分析

食饵—捕食者模型稳定性分析【摘要】自然界中不同种群之间还存在着一种非常有趣的既有相互依存、又有相互制约的生活方式：种群甲靠丰富的天然资源生存，种群乙靠捕食甲为生，形成食饵-捕食者系统，如食用鱼和鲨鱼，美洲兔和山猫，害虫和益虫等。本文是基于食饵—捕食者之间的有关规律，建立具有自身阻滞作用的两种群食饵—捕食者模型，分析平衡点的稳定性，进行相轨线分析，并用数值模拟方法验证理论分析的正确性。【关键词】食饵—捕食者模型相轨线平衡点稳定性一、问题重述在自然界中，存在这种食饵—捕食者关系模型的物种很多。下面讨论具有自身阻滞作用的两种群食饵-捕食者模型，首先根据该两种群的相互关系建立模型，解释参数的意义，然后进行稳定性分析，解释平衡点稳定的实际意义，对模型进行相轨线分析来验证理论分析的正确性。二、问题分析本文选择渔场中的食饵(食用鱼)和捕食者(鲨鱼)为研究对象，建立微分方程，并利用数学软件MATLAB求出微分方程的数值解，通过对数值结果和图形的观察，猜测出它的解析解构造。然后，从理论上研究其平衡点及相轨线的形状，验证前面的猜测。三、模型假设1.假设捕食者（鲨鱼）离开食饵无法生存；2.假设大海中资源丰富，食饵独立生存时以指数规律增长；四、符号说明)(tx/)(1tx——食饵(食用鱼)在时刻t的数量；)(ty/)(2tx——捕食者(鲨鱼)在时刻t的数量；1r——食饵(食用鱼)的相对增长率；2r——捕食者(鲨鱼)的相对增长率；1N——大海中能容纳的食饵(食用鱼)的最大容量；2N——大海中能容纳的捕食者(鲨鱼)的罪的容量；1——单位数量捕食者（相对于2N）提供的供养食饵的实物量为单位数量捕食者(相对于1N)消耗的供养甲实物量的1倍；2——单位数量食饵（相对于1N）提供的供养捕食者的实物量为单位数量捕食者(相对于2N)消耗的供养食饵实物量的2倍；d——捕食者离开食饵独立生存时的死亡率。五、模型建立食饵独立生存时以指数规律增长，且食饵(食用鱼)的相对增长率为1r，即rxx，而捕食者的存在使食饵的增长率减小，设减小的程度与捕食者数量成正比，于是)(tx满足方程axyrxayrxtx)()((1)比例系数a反映捕食者掠取食饵的能力。由于捕食者离开食饵无法生存，且它独立生存时死亡率为d，即dyy，而食饵的存在为捕食者提供了食物，相当于使捕食者的死亡率降低，且促使其增长。设这种作用与食饵数量成正比，于是)(ty满足bxydybxdyty)()((2)比例系数b反映食饵对捕食者的供养能力。方程(1)、(2)是在自然环境中食饵和捕食者之间依存和制约的关系，这里没有考虑种群自身的阻滞作用，是Volterra提出的最简单的模型。下面，我们加入种群自身的阻滞作用，在上两式中加入Logistic项，即建立以下数学模型：22111111)(1NxNxxrtx(3)22112122)(2NxNxxrtx(4)六、模型求解在此，我们采用MATLAB软件求解此微分方程组中的)(1tx、)(2tx的图形及相轨线图形。设5.11，42，11r，4.02r，35001N,5002N,使用MATLAB软件求解，程序代码如下：1)建立M文件functiony=fun(t,x)y=[x(1).*(1-x(1)./3500-1.5*x(2)./500),0.4.*x(2).*(-1+4.*x(1)./3500-x(2)./500)]';2)在命令窗口输入如下命令：[t,x]=ode45('fun1',[0,40],[2000,35])得到数值解如下：t(x(1),x(2))1.0e+003*(单位：千克)00.10330.20660.30990.41320.80791.20262.00000.03502.06540.03692.12760.03892.18630.04122.24120.04382.41130.05602.51110.07321.59731.99192.38062.76933.15793.54663.93534.32394.71265.10125.41675.73226.04776.36316.74467.12617.50767.88918.29678.70429.11179.519310.017310.515311.013311.511312.045312.579313.113313.647314.223114.798915.374715.950516.552317.154117.755918.357719.040319.722920.405421.088021.857422.62682.53580.09612.48700.12482.37410.15772.21270.19222.02280.22461.82470.25141.63600.26971.47270.27931.34310.28131.24270.27751.17800.27171.12940.26421.09510.25611.07280.24761.05990.23771.05930.22861.06840.22061.08510.21391.10940.20811.13760.20381.16760.20091.19750.19931.23130.19901.25990.20001.28150.20201.29540.20471.30220.20791.30170.21111.29580.21381.28640.21591.27450.21741.26270.21811.25290.21811.24550.21771.24020.21691.23770.21601.23750.21511.23910.21441.24200.21381.24540.21341.24840.21341.25060.21351.25240.21371.25300.214023.396224.165625.065625.965726.865727.765728.765729.765730.765731.765732.765733.765734.765735.765736.824337.882838.941440.00001.25270.21421.25200.21441.25090.21451.24990.21451.24950.21441.24950.21431.24960.21431.24980.21421.25000.21421.25010.21431.25010.21431.25010.21431.25000.21431.25000.21431.25000.21431.25000.21431.25000.21431.25000.2143plot(t,x),grid,gtext('x(t)'),gtext('y(t)')图1.数值解)(1tx，)(2tx的图形plot(x(:,1),x(:,2)),grid,图2.相轨线图形从数值解及)(1tx，)(2tx的图形可以看出他们的数量变化情况，随着时间的推移，都趋于一个稳定的值，从数值解中可以近似的得到稳定值为：(1250,214)。下面对其平衡点进行稳定性分析：由微分方程(3)、(4)2211212222111111),(),(2121NxNxxrNxNxxrxxfxxf得到如下平衡点:)0,(11NP,)1)1(,1)1((212221112NNP,)0,0(3P因为仅当平衡点位于平面坐标系的第一象限时(0,21xx)才有意义，所以，对2P而言要求20。按照判断平衡点稳定性的方法计算：)21()21(221122122221112211112121NxNxrNxrNxrNxNxrggffAxxxx根据p等于主对角线元素之和的相反数，而q为其行列式的值，我们得到下表：平衡点pq稳定条件)0,(11NP)1(221rr)1(221rr21)1)1(,1)1((212221112NNP2122111)1()1(rr2121211)1)(1(rr21)0,0(3P21rr21rr不稳定七、模型分析与检验1.平衡点稳定性的分析及其实际意义：1)对)0,(11NP而言，有p=)1(221rr，q=)1(221rr，故当21时，平衡点)0,(11NP是稳定的。意义：如果)0,(11NP稳定，则种群乙灭绝，没有种群的共存。2)对)1)1(,1)1((212221112NNP而言，有p=2122111)1()1(rr，q=2121211)1)(1(rr，故当21时，平衡点)1)1(,1)1((212221112NNP是稳定的。意义：如果)1)1(,1)1((212221112NNP稳定，则两物种恒稳发展，会互相依存生长下去。3)对)0,0(3P而言，由于21rrp，21rrq，又有题知1r0,2r0,故q0,即)0,(11NP是不稳定的。2.平衡点的检验：对于平衡点)1)1(,1)1((212221112NNP，把前面给出的初始值带入，在这使用MATLAB软件进行简单的求解，在命令窗口输入如下代码：x(1)=(3500.*(1+1.5))./(1+1.5.*4);x(2)=(500.*(4-1))./(1+1.5.*4);[x(1);x(2)]ans=1.0e+003*1.25000.2143把此处求解出的解和前面得出的数值解进行比较可知，平衡点)1)1(,1)1((212221112NNP是稳定的。八、模型的评价与推广1.模型的评价自然界中，任何物种即使是捕食者也有自身的阻滞作用，该模型从原始的没带自身阻滞作用模型中加入了阻滞项，使得此模型更接近于生态平衡系统。从此模型中，我们知道两物种同时灭绝是不稳定的，也就是不太可能的，但两种群有一种灭绝一种生存是完全有可能的，两种群共存的可能也是可能的。2.模型的推广本文只考虑两物种模型，我们完全可以把此模型推广到三物种的情形。自然界里长期存在的呈周期变化的生态平衡系统应该是结构稳定的，即系统受到不可避免的干扰而偏离原来的周期轨道后，其内部制约作用会使系统自动回复原状，如恢复原有的周期和振幅，而Volterra模型描述的周期变化状态却不是结构稳定的。要得到能反映周期变化的结构模型，要用到极限环的概念参考文献[1]姜启源，谢金星，叶俊．数学模型,高等教育出版社.2003年[2]冯杰，黄力伟，王勤.《数学建模原理与案例》科学出版社，2007年1月