一类微分方程模型稳定的数值模拟毕业论文

1、一类微分方程模型稳定性的数值模拟摘要： 数值模拟也叫计算机模拟。当前对于数值模拟的研究成果已经应用于诸多领域，如一类微分方程模 型稳定性的数值模拟更是因为它能直观和准确地反映和解决问题而受到广泛的注意。而本文通过搜集数据 和种群模型，一方面借助 MATLAB编写程序从而实现对种群关系问题的直观表述，验证理论结果，另一方 面尽可能地分析模型之间的联系，理解和分析问题的实质。关键词： 数值模拟 稳定性 平衡点 微分方程1 前言20 世纪以来，随着大量的边缘科学诸如电磁流体力学、化学流体力学、动力气象学、海 洋动力学、地下水动力学等等的产生和发展，微分方程稳定性的应用也越来越广泛。尤其当 我们描述实

2、际对象的某些特性随时间或空间而演变的过程 , 分析它的变化规律，预测它的未来 形态时，就要研究微分方程模型的稳定性。但是在理论上研究微分方程的稳定性并不求解微 分方程，而是用微分方程稳定性理论研究平衡状态的稳定性。由于理论本身的抽象性致使这 种动态变化的过程不能形象地表述。而微分方程稳定性的数值模拟则弥补了这样的不足，它 可以借助 matlab 等数学软件对方程解的集合进行几何描述， 同时根据几何图形间接判断微分 方程的稳定性。这是一个从抽象的稳定解的概念到形象的图形表示的过程，它能够直观地通 过方程的解模拟图形，从而很好地帮助分析方程的稳定性。鉴于微分方程稳定性 数值模拟的特点 ，本文主要从

3、以下几个方面准备：首先，对本文的 研究背景和主要工作进行简单的论述；其次，对本文需要的知识点再进行简单的介绍；接着， 从一类关于种群的微分方程模型入手，对这些方程进行数值模拟，同时尝试利用数值模拟的 结果说明在一定范围内图形所具有的稳定性。利用微分方程平衡点稳定性结合图形，对两种 生物种群的模型进行讨论分析。着重讨论微分方程的稳定性，突出数值模拟的优越性和微分 方程平衡点稳定性在实际问题中的重要应用；最后对这些工作进行简单的总结。1.1 研究背景 数值模拟是对具体对象抽取数学模型，然后用数值分析方法，通过计算机求解。经过几 十年的发展，开发了许多不同的科学方法，其中有： （1） 差分法法；（2

4、）有限元法；（3）数 值积分法；（4）蒙特卡洛法等。而数值模拟分析方法则是从结构化矩阵分析发展而来, 逐步推广到板、壳和实体等连续体固体力学分析。近年来数值模拟方法已发展到流体力学、温度场、电传导、磁场、渗流等求解计算 , 最 近又发展到求解几个交叉学科的问题。究其原因主要是数值模拟存在以下特点：首先体现在 它的理论意义 上。它可以通过对复杂或不可观察的现象进行定量分析和对极端情况下尚不知 的规则的推测和预测，实现对复杂现象的模拟，以助于认清现象的本质，弄清整个过程内含 的规律。其次则是体现在 数值模拟的现实意义 上。它可以根据对研究现象和过程的数值模拟， 优化结构设计或者工艺设计， 从而减少

5、试验工作量， 提高产品或研究成果的质量。 可以说 , 继 理论分析和科学试验之后 , 数值模拟已成为科学技术发展的主要手段之一。20世纪八十年代以来，国际数值模拟研究就取得显著发展，它已不仅是一种现代化的实 验手段，而且已发展为具有独立特征的学科分支。而随着计算数学理论和方法的迅速发展和 各种高级计算机语言的出现，使国内数值模拟技术更是得到更好的发展和应用。国内与国外的研究存在一个共同点 , 即对微分方程模型稳定性的数值模拟及它在各个领 域中相关研究和应用的关注。 然而, 在这一现象之下却存在着差异。 国外是建立在微分方程稳 定性研究发展比较充分的基础上具有较大的开拓性和灵活性 , 而国内数值

6、模拟的应用领域在 很大程度上受国外研究方向和方式的影响。国内现在的研究一方面介绍微分方程稳定性理论 在数值模拟方面的扩展应用及国外研究近况 , 另一方面是探讨微分方程稳定性数值模拟在实 践中的运用。基于这些事实可知在以计算机为基础的其他技术带动下，数值模拟技术必将发生更大的 变化，并对微分方程理论证明方向的研究起到相当大的推动作用。而在众多关于微分方程稳 定性模型中，以种 群模型 为代表的一类微分方程稳定性 模型作为数值模拟研究的典 型之 一，它必将会得到更多的 关注。1.2 本文的主要工作本文主要以种群模型为例对一类微 分方程稳定 性模型进行数值模拟的研究。首 先针对相互竞争和相互依存两种种

7、群模 型，对种 群的演变过程进行 数值模拟，验证 理论 结果；接着再对食饵-捕食 者模型进 行分析并 进行数值 模拟，最后尝试改进模型 并运用数值模拟验证结论的可行性。对模型进行改进时，主要考虑的内容如下：因为种群 模 型有相互竞争模型、相互依存及 Volterra 食饵捕食模型三种，所以可以从以下几 个方 面阐 述。（ 1）首先，从两种群的竞争关系入手建立微分方程模型，并从理论上对模型 的动力学性质进行分析，揭示两种竞争种群之间的数量变化关系。同时，结合生物学对模型 作出相应的解释，并运用数值模拟验证结论的可行性； （2）其次，基于前 两种 种群模型 是 目前 研究 和引 用较 广 的模 型

8、， 改进 和完 善 的空 间不 是很 大， 而 Volterra 食 饵捕 食 者 模 型相 对而 言具 有 比较 多 的 局限 性 ， 所以可 以改 进的 空间 比 较大 。（ 3）最 后， 对微分方程结构以及部分参数进行修改，达到 改进和完善“食饵和捕食者”模 型的 目的，使之 可以与前两种模型比较并运用极限环把周期性变 化的结 构体现出来，同 时通过数值模拟直观和准确地揭示问题 的实质，验 证理论讨论的正确性。2 预备知识2.1 二阶微分方程的平衡点和稳定性在这里，我们以二阶微分方程为例，简要叙述微分方程的稳定性理论。下面讨论以下形式的微分方程：dxd1t(t ) f (x1,x2)dt

9、 1 2(1)dxd2t(t) g(x1,x2)dt右端不显含 t ，代数方程组2)f (x1,x2) 0 g(x1,x2) 0的实根(x10,x20)称为方程( 1)的平衡点，记为 P0(x10,x20)。如果从所有可能的初始条件出发，方程( 1)的解 x1(t), x2(t)都满足lim x1(t) tlim x2(t) x0t23)则称平衡点 P0(x10,x20) 是稳定的(渐近稳定)；否则，称 P0 是不稳定的(不渐近稳定) 。为了用直接法讨论方法方程1)的平衡点的稳定性，先看线性常系数方程：dx1(t)1a1x1b1x2dt1 11 2dx2 (t ) a2x1 b2x2dt2 1

10、2 24)系数矩阵记作：b1b2并假定 A的行列式 detA 0。于是原点 P0 (0,0) 是方程( 4)的唯一平衡点，它的稳定性由 的特征方程： det(A I) 0的根 (特征根)决定，上方程可以写成更加明确的形式：2 p q 0p (a1 b2)(5)q detA将特征根记作 1, 2，则： 1, 2 1( p p2 4q)(6)方程( 4)的解一般有形式 c1e 1t c2e 2t( 1 2)或(c1 c2t)e t( 1 2)， c1, c2为任意实数。由定义(3)，当 1, 2 全为负数或有负的实部时 P0(0,0) 是稳定的平衡点， 反之，当 1, 2 有 一个为正数或有正的实

11、部时 P0(0,0) 是不稳定的平衡点。微分方程稳定性理论将平衡点分为结 点、焦点、鞍点、中心等类型，完全由特征根1, 2 或相应的 p, q取值决定，下表简明地给出了这些结果，表中最后一列指按照定义( 3)式得下列关于稳定性的结论。表1 ：由特征方程决定的平衡点的类型和稳定性1, 2p,q平衡点类型稳定性12p 0,q 0, p2 q稳定结点稳定122p 0,q 0, p q不稳定结点不稳定10q0鞍点不稳定12p 0,q 0, p2 q稳定退化结点稳定3122p 0,q 0, p q不稳定退化结点不稳定1, 2p 0,q 0, p2 q稳定焦点稳定1, 22p 0,q 0, p q不稳定焦

12、点不稳定1, 2p 0,q 0中心不稳定由上表可以看出，根据特征方程的系数 p,q 的正负很容易判断平衡点的稳定性，准则如 下 ： 若 A1 ： p 0,q 0 则 平 衡 点 稳 定 。 若 A2 ： p 0 或q 0 则平衡点不稳定。以上是对线性方程( 4)的平衡点 P0 (0,0) 稳定性的结论，对于一般的非线 性方程(1)，可以用近似线性方法判断其平衡点 P0(x10,x20)的稳定性，在 P0(x10,x20)点将 f(x1,x2) 和 g(x1,x2) 作泰勒展开，只取一次项，得( 1)的近似线性方程(7)dx1(t) f x1 ( x10 , x20 )( x1 x10) fx2

13、 ( x10 , x20 )( x2 x20) dt12系数矩阵记作：fx1fx2AP (x0,x0 )gx1gx2 P0(x10,x20 )特征方程系数：p (fx1 gx2)P0(x10,x20), q det A 。显然， P0 (x10 , x20 )点对于方程( 7)的稳dx2(t) gx1 (x10,x20)(x1 x10) g x2 ( x10 , x20 )( x2 x20) dt12定性由表 1 或准则 A1、A2决定，而且已经证明了如下结论 :若方程( 7)的特征根不为零或实部不为零，则 P0 (x10, x20 )点对于方程( 1)的稳定性与对 于近似方程( 7)的稳定性

14、相同。这样， P0(x10,x20) 点对于方程( 1)的稳定性也由准则 A1、A2 决定。2.2 微分方程与极限环及其稳定性设系统8)dddxyt P(x,y) ddyt Q(x,y)具有闭轨线 C。假如在 C 充分小邻域中，除 C之外，轨线全不是闭轨线，且这些非 闭轨线当 t或 t 时趋近于闭轨线 C，则说闭轨线 C是孤立的，并称之为 (15) 的一个极限环 . 极限环 C 将相平面分成两个区域：内域和外域。如果极限环 C的内域的靠近 C的轨线当 t +(- ) 时盘旋地趋近于 C，则称 C是内 稳定 ( 内不稳定的 ) ；如果在极限环 C的外域的靠近 C的轨线当 t +( ) 时盘旋地趋

15、近于 C，侧称 C是外稳定的 （外不稳定的 ）；如果当 t （ ）时， C的内部及外部靠 近 C 的轨线都盘旋地趋近于 C，则称 C 是稳定的 （ 不稳定的 ）, 如果当 t （ ） 时， C 的内外部的稳定性相反，则称C 为半稳定的。3 一类微分方程模型稳定性的数值模拟3.1 模型 : 相互竞争的种群模型1、模型假设及说明 :（1）、假设有甲乙两个种群，当它们独自在一个自然环境中生存时，数量的演变均遵从 Logistic 规律，记 x1（t）, x2（t）是两个种群的数量， r1 ， r2是它们的固有增长率， N1, N2是它 们的最大容量。（2）、当甲、乙两个种群在同一自然环境中生存时，乙

16、对甲增长的阻滞作用与乙的数量成 正比；甲对乙有同样的作用。（3）、 1 ：单位数量乙（相对于 N2 而言）消耗的供养甲的食物量为单位数量甲 （相对于 N1 而言）消耗的供养甲的食物量的 1倍。 2 表示的意义与之相反。4）、x1 可解释为相对于 N1 而言单位数量的甲消耗的供养甲的食物量 （设食物总量为 1） N1Nx21可做类似解释2、模型的建立x1(t) r1x1(1 N1 1 N2)x2(t) r2x2(1 2 Nx1 Nx2)3、稳定性分析： 利用平衡点的稳定性分析，讨论时间足够长以后两种群的变化趋势。分 析结果列入表 2。表 2：种群竞争模型的平衡点及稳定性平衡点pq稳定条件p1 (

17、N1 ,0)r1 r2(1 2)r1r2(1 2 )1 1 ， 2 1p2 (0, N2 )r1(1 1) r2r1r2(1 1)1 1 ， 2 1p (N1(1 1), N2(1 2) p3( 1 1 2 , 1 1 2 ) 1 2 1 2r1(1 1) r2(1 2)r1r2(1 1)(1 2)1 1 ， 2 11 1 21 1 2P4(0,0)(r1 r2)r1r2不稳定根据相轨迹的性质， 相关文献中对于 1 ， 2的不同取值范围的情况已经进行了详尽的分析，概括如下：（1）、 1 1， 2 1。由表 2 知对于 p1（N1,0） 有 p0, q0。 p1（N1,0） 稳定； （2） 1

18、1， 2 1 。类似的分析可知 p2（0, N2）稳定；（3） 1 1， 2 1。在 p3点 p0, q0，故 p3稳定；(4) 1 1， 2 1。由表 2可知在 p3点，q0,故 p3不稳定(鞍点)。4、数值模拟根据表 2 中的稳定条件，探究 1 ， 2 在不同的取值范围内的情况。通过数值模拟验证 其理论结果。探究 1：设已知 1 =0.5 ， 2 =1.6 ，r1 =2.5 ，r2 =1.8 ，N1 =1.6 ，N 2 =1，初始条件 1：x1(0) 0.1， x2(0) 0.1；初始条件 2： x1(0) 1， x2(0) 2 。(1) 首先，建立 m-文件 model1.m 如下：fu

19、nction dx=model1(t,x) a1=0.5;a2=1.6;r1=2.5;r2=1.8;N1=1.6;N2=1;dx=r1*x(1).*(1-x(1)./N1-a1*x(2)./N2);r2*x(2).*(1-x(2)./N2-a2*x(1)./N1) 2)其次，建立主程序 Untitled1.m 如下：t,x=ode45(shier1,0 30,0.1 0.1);subplot(2,2,1)plot(t,x(:,1),t,x(:,2),r);xlabel(时间 t);ylabel(种群密度 x);gtext(x1(t);gtext(x2(t);grid on;title(当初值分

20、别为 0.1,0.1 时);subplot(2,2,2);plot(x(:,1),x(:,2),r); grid on;title(相轨线(x1(t),x2(t);t,x=ode45(model1,0 30,1 2);subplot(2,2,3);plot(t,x(:,1),t,x(:,2),r);xlabel(时间 t);ylabel(种群密度 x);gtext(x1(t);gtext(x2(t);grid on;title( 当初值分别为 1,2 时); subplot(2,2,4);plot(x(:,1),x(:,2),r); grid on;title(相轨 (x1(t),x2(t);

21、时间t相 轨 线 (x1(t),x2(t)(图形： 1)探究 2：设已知 1 =1.5 ， 2 =1.6 ，r1 =2.5 ，r2 =1.8 ，N1 =1.6 ，N 2 =1，初始条件 1：x1(0) 0.1， x2(0) 0.1；初始条件 2： x1(0) 1， x2(0) 2 。2x度1x1(t)度密密群种x2(t)0种-10102030时间t当 初 值 分别为1,2 时21x度 度密x2(t)x1(t)群种0-10种102030当 初 值 分 别 为 0.1,0.1 时时间t相 轨 线 (x1(t),x2(t)图形： 2)探究 3：设已知 1=1.5 ， 2 =0.6r1 =2.5 ，r

22、2 =1.8 ，N1 =1.6 ，N 2 =1，初始条件 1：x1(0)0.1，x2(0) 0.1；初始条件2： x1(0) 1，x2(0) 21.5-0.5x密度 1群种 0-1x2(t)x1(t)值 分 别 为 0.1,0.1时10.50当初10 20 30 时间t当 初 值 分 别 为 1,2 时x2(t)x1(t)301020相 轨 线 (x1(t),x2(t)时间t探究 4：设已知 1=0.52 =0.6(图形： 3)r1 =2.5 ，r2 =1.8 ，N1 =1.6 ，N 2 =1，初始条件 1：x1(0)0.1，x2(0) 0.1；初始条件2： x1(0) 1，x2(0) 2值

23、分 别 为 0.1,0.1时间t当初相 轨 线 (x1(t),x2(t)时图形：4)5、结果解释根据数值模拟的结果。 （1）当 1 1， 2 1时，如图 1 所示。 1 1 意味着在供养的资 源的竞争中乙弱于甲， 2 1意味着在对供养乙的资源的竞争中甲强于乙，于是种群乙终将 灭绝，种群甲趋于最大容量，即 x1（t）、x2（t）趋向平衡点 p1（N1,0） 。（2）当 1 1， 2 1时， 如图 2 所示。在竞争甲的资源中乙较强，而在竞争乙的资源中甲较强。随着时间的推移，有 可能甲占优势，乙趋于灭亡。但也有可能乙占优势，而甲灭亡。究竟趋于哪个平衡点，由初 始时刻两群生物的总数决定。 （3）当 1

24、 1 ， 2 1 时，如图 3 所示。这种情况和（ 1）刚好相反。（4）当 1 1， 2 1时，如图 4 所示。因为在竞争甲的资源中乙较弱，而在竞争乙的 资源中甲较弱，于是可达到一个双方共存的稳定的平衡状态p3，这是种群竞争中很少出现的情况。显然数值模拟的结果可以很好地解释和说明理论结果。3.2 一类相互依存的种群模型自然界中两种群相互依存有三种形式： （ 1）甲可以独自生存，乙不能独自生存；甲乙一 起生存时相互提供食物、促进增长； （2） 甲乙均可以独自生存；甲乙一起生存 时相互提供食 物、促进增长；（ 3） 甲乙均不能独自生存；甲乙一起生存时相互提供食物、促进增长。下面 分别对着三种情形进

25、行讨论。3.2 模型 1：种群的相互依存（情形 1：）1、模型假设及说明：（1）以x1（t）、x2（t）表示甲、乙二种群在时刻 t的数量， r1表示甲种 群的固有增长率 , Ni（i 1,2） 分别表示甲、乙二种群在单种群情况下自然资源所能承受的最大 种群数量；（2）甲乙一起生存时，乙为甲提供食物、促进增长； （ 3）乙种群没有甲的存在会 灭亡,死亡率为 r2 ,甲乙一起生存时甲为乙提供食物、促进增长；乙的增长又受到本身的阻滞 作用；（4）乙为甲提供食物是甲消耗的 1 倍, 甲为乙提供食物是乙消耗的 2倍。2、模型建立x1(t) r1x2(1 x11 x2N1N2x2(t) r2x1( 1 2

26、 x1 )N13、稳定性分析： 利用平衡点的稳定性分析，讨论时间足够长以后两种群的变化趋势。分 析结果列入表 3。表 3 种群依存模型的平衡点及稳定性平衡点pq稳定条件P1 N1,0r1 r2 ( 2 1)r2r2( 2 1)2 1,1 2 1N1(1 1) N2( 2 1)r1 1 1 r2( 2 1)r1r2(1 1) 2 11 1, 2 1,p2(, )1 1 2 1 1 21 1 21 1 21 2 1P3(0,0),r1 r2r1r2不稳定8显然， p2是甲乙相互依存而共生的平衡点，下面我们着重分析 p2 稳定的条件。由 p2的 表达式容易看出，要使平衡点 p2 有实际意义，即位于相

27、平面第一象限，必须满足下面两个条 件中的一个：A1 : 1 1, 2 1, 1 2 1A2 : 1 1, 2 1, 1 2 1 由上面的分析知：仅在条 A1件下 p2 才是稳定的，而在 A2条件下 p2是不稳，而是鞍点 关于这一问题，在相关文献中已经做了详尽的描述。4、数值模拟根据表 3中的稳定条件，探究 1 ， 2 在不同的取值范围内的情况。通过数值模拟验证其 理论结果。探 究 1 ： 设 1 0.1, 2 3,r1 1.8,r2 1.5,N1 1.6,N 2 1 ， 初 始 值 分 别 取 ：x1(0) 0.1,x2(0) 0.1 x1(0) 1,x2 (0) 2。(1) 首先，建立 m-

28、文件 model1.m 如下： function dx=model1(t,x)r(1)=1.8;r(2)=1.5;a=0.1;b=3;N(1)=1.6;N(2)=1;dx=r(1).*x(1).*(1-x(1)/N(1)+a.*x(2)/N(2);r(2).*x(2).*(-1+b.*x(1)/N(1) -x(2)/N(2);(2) 其次，建立主程序 Untitled1.m 如下：t,x=ode45(model1,0 8,0.1 0.1);subplot(2,2,1) ; plot(t,x);gtext(x1(t);gtext(x2(t);grid on;title( 初值为 0.1 0.1

29、时);subplot(2,2,2);plot(x(:,1),x(:,2),grid on; title(初值为 0.1 0.1 时种群相轨线 );初值为 1 2 时初值为 1 2 时种群相轨t,x=ode45(model1,0 8,1 2);subplot(2,2,3); plot(t,x);gtext(x1(t);gtext(x2(t);grid on;title( );subplot(2,2,4);plot(x(:,1),x(:,2),grid on; title( 线);初 值 为 0.1 0.1 时1.52.5图形 5)r2 1.5,N1 1.6,N初始值分别取：1 1.1, 2 1.

30、 3, r1x1(0) 0.1,x2(0)(图形：1.8,r26)1.5,N11.6,N 2 1始值分别取：0.1初值 为 0.1 0.1时x1(0) 1,x2 (0) 2。0.40.30.20.101.521.510.5012x 10314初值为 0.1 0.1相x 1初0 值14 为 1 2图形：7)1 3, 2 0.1,r1 1.8,r2 1.5,N1 1.6,N2值分别取：x1(0) 0.1,x2(0)0.1初值 为 0.1 0.1时x1(0) 1,x2 (0) 2。12345图形：8)105、结果解释根据数值模拟的结果。 （1）当 11, 1 2 1时，如图 5 所示。甲为乙提供足够

31、的食 物，甲本身食物不足。 随着时间的增长， 甲、乙于是达到一个双方共存的稳定的平衡状态 p2 ； （2）当 11, 21, 21 时，如图 7 所示。由于甲、乙都能给对 方提供足够的食物，因此甲和乙有可能共存也有可能乙灭绝，究竟趋于哪个平衡点，由初始 时刻两群生物的总数决定； （4）当 11, 21， 1 2 1时，如图 8 所示。乙为甲提供足够的 食物，乙本身食物不足。种群乙灭绝，甲达到最大种群数量。显然数值模拟的结果可以很好 地解释和说明理论结果。3.2 模型 2：相互依存的种群模型（情形 2）1、模型假设及说明： 与第一种情况的假设一样，只需要将 r2 修改为固有增长率。x1x2x1(

32、t) r1x1(1 1 1 2 )N1N2x1 x2 2 N12、模型建立x2（t） r2x2（ 1 2 N11 N22）3、稳定性分析：析结果列入表 4。利用平衡点的稳定性分析，讨论时间足够长以后两种群的变化趋势。分表 4 独立种群相互依存模型的平衡点和稳定性Pipq稳定条件p1(N1,0)r1 r2(1 2 )r2 r1 (1 2 )不稳定p2(0,N2)r2 r1(1 1)r2r1(1 1)不稳定p ( 1 1)N1 ,( 2 1)N2) p3( 1 2 1 , 1 2 1 )r1(1 1) r2(1 2)r1r2(1 1)(1 2 )1 2 11 1 21 1 2p4(0,0)r1 r

33、2r1r2不稳定由 P3点的表达式容易看出，要使平衡点 P3有实际意义，即位于相平面第一象限（ x1,2 0）， 必须满足下面两个条件中的一个：A1： 11， 1 21， 2 1， 1 21由表 4 可知，仅在条件 A1下 P3才是稳定的。4、数值模拟根据表 4中的稳定条件，探究 1 ， 2 在不同的取值范围内的情况。通过数值模拟验证其 理论结果。探 究 1 ： 设 1 0.1, 2 1.6, r1 1.8,r2 1.5,N1 1.6,N2 1 ， 初 始 值 分 别 取 ：11x1(0) 0.1,x2(0) 0.1 x1(0) 1,x2 (0) 2。(1) 首先，建立 m-文件 model5

34、.m 如下：function dx=model5(t,x) r(1)=2.5;r(2)=1.8;a=0.1;b=1.6;N(1)=1.6;N(2)=1; dx=r(1).*x(1).*(1-x(1)/N(1)-a.*x(2)/N(2);r(2).*x(2).*(1-b.*x(1)/N(1)- x(2)/N(2);(2) 其次，建立主程序 Untitled3.m 如下： t,x=ode45(model5,0 8,0.1 0.1);subplot(2,2,1) plot(t,x);gtext(x1(t);gtext(x2(t);grid on;title(初值为 0.1 0.1 时);subplo

35、t(2,2,2)plot(x(:,1),x(:,2),grid on;title( 初值为 0.1 0.1 两种群相轨线 ); t,x=ode45(model5,0 8,1 2);subplot(2,2,3) plot(t,x);gtext(x1(t);gtext(x2(t);grid on;title( 初值为 1 2 时);subplot(2,2,4) plot(x(:,1),x(:,2),grid on; title(初值为 1 2 两种群相轨线 );初 值 为 0.1 0.1 时0.40.30.20.10初值为 0.1 0.1两种群相 轨线0 0.5 1 1.5 2初 值 为 1 2

36、时初值为 1 2两 种群相轨线(图形： 9)探究 2：设 1 1.6, 2 0.1,r1 1.8,r2 1.5,N1 1.6,N2 1，初始值分别取： x1(0) 0.1,x2(0) 0.1 x1(0) 1,x2 (0) 2。120.1 0.1 时10)探 究 3：1 1.1, 2 1.6, r1 1.8,r2图形：1.5,N1 1.6,N初始值分别取：x1(0) 0.1,x2(0)0.12初 值 为 0.1 0.1 时x1(0) 1,x2 (0) 2。11)探 究 4：1 0.1, 2 0.1, r1 1.8图形：1.5,N1 1.6,N2初始值分别取：x1(0) 0.1,x2(0)0.1初

37、 值 为 0.1 0.1时x1(0) 1,x2 (0) 2。1图形：12)5、结果解释 根据数值模拟的结果：1)当 1 1,1 2 1, 21 ， 1 21时，如图 11 所示.种群乙灭绝，种群甲达到最大种群数量；（4）当 1 1, 2 1， 1 21, 21 ， （3）当 1 1, 2 1， 1 21 2 1时，如图 13 所示。种群甲、乙都种群甲灭绝， 种群乙达到最大种群数量；（ 1）当 1 1,1 21时，如图 15 所示。种群甲、乙都灭绝；（4） 当 1 1, 2 1，2 1 ，14 所示。1 21时, 如图 16 所示。甲、乙两种群将分别趋于非零的有限值，达到平衡状态。显然数 值模拟

38、的结果可以很好地解释和说明理论结果。3.3 模型 1： Volterra 食饵捕食者模型（不考虑人工捕捞）1、模型假设及说明： 假设 x1（t）表示食饵在时刻 t 的数量， x2（t） 表示捕食者在时刻 t 的 数量，即， r1表示食饵独立存在是以指数规律增长，与捕食者相对的增长率，r2 表示捕食者独立存在时的死亡率， 1表示单位数量的乙（相对于 N1 ）捕食单位数量甲（相对于 N2）的 能力， 2表示单位数量的甲（相对于 N 2）供养单位数量乙（相对于 N1）的能力， N1表示生 存环境允许食饵的最大生存量， N2 表示生存环境允许捕食者的最大生存量。2、模型建立：x1(t) r1x1(1

39、x11 x2 )N1N2x2(t) r2x2( 1 2 x1 x2 )N1 N23、稳定性分析：利用平衡点的稳定性分析，讨论时间足够长以后两种群的变化趋势。分析结果列入表 616表 6 独立种群“食饵捕食者”模型的平衡点和稳定性平衡点pq稳定条件p1(0,0)r1 r2r1r2不稳定p2 (N1,0)r1 r2( 2 1)r1r2 ( 2 1)21N1( 1 1) N2( 2 1)r1( 1 1) r2( 2 1)r1r2( 1 1)( 2 1)21p3,1 1 2 1 1 21 1 21 1 2表中可以看出 p1(0,0) 是不稳定的点；当 2 1时， p2(N1,0) 点是稳定点；当 2

40、1时，p3 N1( 1 1),N2( 2 1) 是稳定点。1 1 2 1 1 24、数值模拟根据表 5中的稳定条件，探究 1 ， 2 在不同的取值范围内的情况。通过数值模拟验证其理论结果。探究 1：设 r1 1,r2 0.5,a(1) 1 1.6,a(2)2 1.6, N 1 100 , N 2 35 ，x0 =25，y0 =2则用 MATLAB软件编程结果如下：90x1(t)x2(t)80706050403020100102030初 值 为 25 2120(x1(t),x2(t)5.554.543.532.521.5406080100初 值 为 25 2 的 相 轨 线探究2：设 r11,r

41、2 0.5,a(1) 1图形：0.6,a(2)17)2 1.6, N1 100 ， N2 35， x0 =25，y0=2。探究y0=2。10090807060504030201003：设 r1x1(t)x2(t)102030初 值 为 25 2(x1(t),x2(t)4060801001110987654321初 值 为 25 2 的 相 轨 线20图形：18)1,r2 0.5,a(1) 1 1.6,a(2)2 0.6, N 1 100 ， N2 35 ， x0 =25，17100初 值 为 25 221.81.6(x1(t),x2(t)1.41.210.80.60.40.2020406080

42、100初 值 为 25 2 的 相 轨 线19)图形：探究 4：设 r1 1,r2 0.5, a(1) 1 0.6,a(2) 2 0.6, N1 100 ， N2 35， x0 =25，y0=2。10090x1(t)8070605040302010x2(t)00102030初 值 为 25 22初值为 25 2的相轨线1.81.61.4t),x2(t)1.210.80.60.40.2020406080100(图形： 20)5、结果解释根据数值模拟的结果： (1)当 2 1时，如图 17 和图 18 所示， p3点稳定。食饵和捕食 者共存，无论是食饵和捕食者的数量如何， 只要同时存在二者就能分别

43、趋向于非零的有限值， 而且始终都不能达到自己环境所允许的最大值； (2) 2 1时，如图 19和图 20 所示， p2点 稳定。此时因为食饵供养捕食者的能力低于捕食者自身的基本生存要求，所以捕食者就慢慢 灭绝了，而食饵则趋向环境允许的最大容量。显然数值模拟的结果可以很好地解释和说明理 论结果。3.3 模型 2:( 考虑人工捕捞 )1、模型假设及说明： 与模型 1 的假设一样，只需添加 m1捕食者掠取食饵的能力， m2表 示食饵对捕食者的供养能力， e 表示捕获能力系数。2、模型建立x1(t) x1r1 (1 x1 ) ex1 m1x1x2xN1x2(t) x2r2 x2 ex2 m2x1x22

44、 2 2 N2 2 2 1 23、数值模拟设r1 1,r2 0.5, (1) 0.1, (2) 0.02,N1 100,N2 35 ， x0=25， y0 =2则用MATLA软B 件编程如下：(1) 首先，建立 m-文件 :7.m 如下：18function dx=model8(t,x)r1=1;r2=0.5;m1=0.1;m2=0.02;e1=0.1;N1=100;N2=35;dx=x(1)*r1*(1-(x(1)/N1)-e1*x(1)-m1*x(1)*x(2);-r2*x(2)+(x(2)/N2)-e1*x(2)+m2 *x(1)*x(2);function dx=model9 (t,x

45、) r1=1;r2=0.5;m1=0.1;m2=0.02;e2=0.3;N1=100;N2=35; dx=x(1)*r1*(1-(x(1)/N1)-e2*x(1)-m1*x(1)*x(2);-r2*x(2)+(x(2)/N2)-e2*x(2)+m2 *x(1)*x(2);(2) 其次，建立主程序 Untitled7.m 如下：t,x=ode45( model8 ,0 50,25 2);subplot(2,2,1) plot(t,x(:,1),t,x(:,2) ,gtext( x1(t) ),gtext(x2(t),gridon;subplot(2,2,2)plot(x(:,1),x(:,2), ,gridt,x=ode45( model9 ,0 50,25 2);subplot(2,2,3) plot(t,x(:,1)