动物群落的稳定发展

《动物群落的稳定发展》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

动物群落的稳定发展摘要：本文通过对某公园近两年内被运出的某种动物的年龄和性别的数据进行统计分析，并针对题目的四个问题分别建立了符合实际的数学模型，在模型的求解过程中，应用C语言进行编程调试，通过统计学软件SAS,数学软件MATLAB等计算工具，编写相应的程序，对建立的模型进行求解，得出了符合实际的结果。问题一：我们假设新生幼仔的数量为小，然后通过对各年龄阶段的存活率八被运走的动物数量吗以及该动物的总体数量的分析来建立该群落的动态变化模型驾®=£号）-£＞漂，利用该群落近两年内被运走的各年龄阶段的个体数dt（=ii=i量分布，用C语言编程计算，推测出当前该动物的年龄结构（具体结果见7页表一）。并利用MATLAB软件对得出的数据用图形表示，利用对比分析法，得到该动物群落的基本分布轨迹，最后用统计软件SAS对模型进行相关性的分析检验，求得相关系数R与P的值，验正了模型的稳定性。问题二：由于现在采用注射避孕药的方法来维持该种群的稳定，而且已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数-［（1-72纪2+。3］（该种群每年的新生幼仔的数量减去该年死亡个体的数6060量的差值）；另外从£=出（即年头的数量与该年年底的数量的差值）（=11=1当£趋于0时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转化为求单目标的最优化问题建立模型；利用MATLAB对模型进行求得，得出当不考虑不确定性因素影响时要注射药物的雌性动物数量为276头，而当考虑了双胞胎和被重复注射这两个不确定性因素影响后，得到要注射药物的雌性动物数量为352头，其中有110头是被重复注射的。问题三：其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b）,即目标函数£应考虑上被运走的数量，即只是对问题二的模型进行扩6060充建立新的目标模型；£\=-［（l-/2）c2+c3］—b和£-£成,-b；利i=l/=1用MATLAB对不同b值进行求解，从而得出相应的避孕措施。（具体结果见19页表二）问题四：我们引进了增量加速度的概念，利用c语言进行编程求解，然后用MATLAB软件对得到的数据进行线性回归分析，得到该群落在减少至M时重新壮大该动物群落能力的模型：M=3.9010+0.0047Do最后应用统计软件SAS对模型进行稳定性分析。关键字：存活率年龄结构新生幼仔数稳定性最优目标增量加速度一.问题重述与提出位于非洲某国的国家公园中栖息着近11000头某种野生动物。管理员要求有

1一个健康稳定的环境以便维持这个11000头该动物的稳定群落。过去的20年中，整个该动物群是通过一些偷猎枪杀以及转移到外地而稳定下来的。但是近年来，偷猎被禁止，而且每年要转移这些动物也比较困难，因此，要控制现在的数量就使用了一种避孕注射法。用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕。要探讨这种避孕注射法的实用性，我们需要完成以下问题：1.探讨该动物年龄在2岁到60岁之间的合理的存活率的模型，推测这个动物群落的当前的年龄结构。2.估计每年在该群落中有多少雌性动物要注射避孕药，可以式群落固定在11000头左右。这里不免有些不确定性，也要估计这种不确定性的影响。3.假如每年转移50至300头此动物到别处，那么上面的避孕措施将可以有怎样的改变？4.如果由于某种原因，突然使得注射避孕的方法不得不停止（例如由于一场灾难导致大量该动物的死亡），那时重新壮大该动物群的能力如何？二.基本假设与符号说明（-）模型假设1.该公园是非开放式的，它与外界不发生关系，从而构成独立的生物群落,该动物群落不存在与其它动物种群的竞争，或虽有竞争，但其影响只局限于该动物群落的死亡率内。2.种群是通过雌性个体的繁殖而增长的，所以用雌性个体数量的变化为主要研究对象。3.为了讨论的必要，我们把新生的幼儿的存活率定为75%,而其后的存活率为95%,直到60岁为止。各年龄组的该动物经过一年后即进入高一级的年龄组，而龄超过60即认为全部死亡，退出该系统。4.由于该公园加强了对该动物群落的保护，我们认为该动物没有再被偷猎射杀。而该动物群落个体数量的减少只是因为自然死亡以及被运走。5.假设同一年龄组的动物个体之间是同质的，我们只考虑其平均水平，不

2讨论个别差异。1.题设该动物在10〜12岁开始怀孕，我们这里设定为11岁开始，经过22个月(约两年)的怀孕期后生幼仔，即可认为该雌性动物在13〜60岁的时间内可以生幼仔。2.该群落的自然死亡是在生完幼仔后才发生的，产幼仔只发生在每年的年初时段，而被运走只发生在年底时段。(二)符号说明%：新生幼儿的存活率，其值为0.75；%：1〜60岁个体的存活率，其值为0.95；片：双胞胎出生的儿率，其值为0.0135；X">(r)：该动物第k年时刻的数量；x摺：该动物第k年初i龄动物的数量；瑞：该动物第k年初底i龄动物的数量第j年被运走的动物的数量；Xa>(0)：表示该动物第k年初时的总数量；4:表示每年没有注射避孕药的雌性动物生幼仔的儿率，其值为二-；3.54：表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的儿率，其值为」5.54：表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的儿率，其值q：表示从13〜60岁该动物的雌性个体的总数;J：表示从1〜59岁该动物的个数总和;G：表示60岁该动物的个体总和；y,：表示13〜60岁雌性动物没有被注射避孕药部分的数量；%：表示13〜60岁雌性动物被注射过避孕药但在两年内不再被注射部分的数量；%：表示13〜60岁雌性动物被注射过避孕药但在两年内被重复注射部分的数量;

3马：表示每年出生幼仔的数量与该年个体死亡的数量的差值；?：表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量的和的差值；表示该动物群落在年底时的总数量与年初的数量加上被运走的个体数量b的差值。三.问题分析与模型建立问题一：1.我们要研究该动物群落的稳定性问题，首先要根据存活率确定其当前的年龄结构。该动物的新生幼仔存活率较低，题设是70%到80%之间，为了讨论的需要，我们这里设定为75%。在1岁后的存活率比较高，在这里设为95%,直到60岁，而超过60岁则认为退出该系统。因此，我们先建立出该动物群落中年龄在2岁到60岁之间的合理的存活率的模型。

4模型一:J6060“八⑺-VY*)_Vy(&Ft一乙心」一乙七，0ai»=iz=i60Z瑞=x(qo)⑴(2)i=\Bjdt（3）=瑞).・・・・・0=0,1,..59)(4)瑞)=0(/=60)(5)瑞)=%瑞)……(，=。)⑹燔)=%瑞).・・・・・(«=1,2,...,59)(7)60y淄C=A(8)式（1）表示该动物第k年增长的数量；式（2）表示该动物第k年初时的总数量，可由已有的数据计算出X（*）（0）来；式（3）表示该动物被运走的数量；式（4）和（5）表示该动物第i龄到了年底全部转化为（i+1）龄；式（6）和（7）表示该动物各年龄段的变化；式（8）表示该动物新生的幼仔数量。2.通过对该公园近两年内从这个地区运出的该动物的年龄和性别的数据进行统计分析，并利用编程工具TurboC2.0对该模型进行编程计算（源程序及计算过程见附录1）,可得到当前该动物群落的年龄结构，如下表所示：

5表一该动物的年龄结构统计表前一年数量前一年运走前一年剩下前雨年数前雨年运走前两年剩下假设无运走今年数年龄（岁）（头）数量（头）数量（头）量（头）数量（头）数量（头）数量（头）量（头）079278680680674680015940594604060455960025640564573205535315693535053554421523504540450835055161350347851254824478490124784544866457745046513452431461743420414441224194094378412940341814404388415939115376397403573683941037193623771436334937411352223303582633233135512334333134013327314337133172329432214308298320143015296305272782833031528513272289328626828716270212492741425025427217256025626012248241258182432222124620226228245192301421623325208216232202185213221172042052202120713194209141951942082219610186198101881841972318601861880188174187

6242526272829303132333435363738394041424344454647484950511761316317821761651771673013716931661561681581414416041561481591501213815141471401511420142143314013214313420114135213312513512761211283125118128120311712113108112121113510811416981061141078991081395100108101128910210929410295108596108689969038791127984918577886167079868014668112697581751065761066717671165572126067726721466819496368631350641351596459104960243656605612445617395356536475316375053503475025254750476414712354447449354445-141444113284123183841381028383443638363333613233436

75234628341618323453322111321022303254301515301713283055284242813152628562613132613132426572410142412122224582232-1022319202259201462022-2182060180181820-21718牛.幼仔的雌性数量78479273478913—60岁雌性数量27352322503277230224702569总数量117146221109211876876110001100611808注1：0岁表示新生幼仔。注2：由于每个年龄段的数据均为推测值，而实际上运走的各年龄段的数量不一定全部与预测值相符，故表中“剩下数量”两组数据中出现负数可认为是独异点，不影响模型整体的准确性。601/2每年新生幼仔的数量（%）减去生幼仔的雌性的数量（:一，由于雄性与雌性的数量比接近1：1,我们可近似地认为13——60岁个体的雌雄数量相等）,其差值即为双胞胎的数量，这个差值与生幼仔的雌性数量之比即为双胞胎的儿率（）1.35%。由表中数据可得，792-784«0.0102,806-792«0.0177,2784792746-734«nn_7890.0163,丝二上。0.0139,这些比例都基本上接近题设的双胞胎的735789儿率0.0135,说明以上推测得出的数据是准确的。

8利用Matlab软件对以上四组数据用图形表示，并进行比较，得到该动物群落的基本分布情况图（源程序见附录2）,如下图所示图1分析该图，可以看出，这四组曲线的轨迹、分布情况基本相同。由于“预测当前的年龄结构情况（无运走）”一组数据没有减去被运走的个体数量，故其每个年龄层的数量都略多于前三组的数量，因此其曲线比前三组的曲线略高一点，利用SAS软件对模型进行相关性的分析检验（源程序见附录3）,得到如下结果:21:39Saturday,ApriI24,2005TheSASSystemTheCORRProcedure4Variables：xlx2x3x4SimpleStatisticsVariableNMeanStdDevSumMinimumMaximumxl60171.66667162.2352010300-10.00000594.00000x260169.90000163.7510210194-2.00000604.00000x360171.00000151.456991026017.00000559.00000x460183.46667162.157331100818.00000600.00000PearsonCorrelationCoefficients,N=60Prob>Ir|underHO：Rho=0xlx2x3x4x11.000000.997440.998960.99896

9<.0001<.0001<.0001x20.997441.000000.998460.99854<.0001<.0001<.0001x30.998960.998461.000000.99999<.0001<.0001<.0001x40.998960.998540.999991.00000<.0001<.0001<.0001图2程序的分析及统计结论：程序中的xl是前一年的该动物群落的年龄结构，x2是前两年该动物群落的年龄结构，x3是该动物群落没有被运出是的年龄结构，x4是预测的当前的该动物群落的年龄结构。过程中的PROCCORR是分析变量中两两变量之间的PEAROS简单相关的。输出结果中的结果1是一些基本的描述统计量，结果2是两两变量之间的相关矩阵，其中包括相关系数和显著性检验的概率。由结果可知前一年的该动物群落的年龄结构(xl)与前两年的该动物群落的年龄结构(x2)的相关系数R=0.99744,P=0.000K0.01,所以前一年的该动物群落的年龄结构(xl)与前两年的该动物群落的年龄结构(x2)之间存在着极显著的正相关；前一年的该动物群落的年龄结构(xl)与没有运走是的该动物群落的年龄结构(x3)的相关系数R=0.99896,P=0.000K0.01,所以前•年的该动物群落的年龄结构(xl)与没有运走时的该动物群落的年龄结构(x3)之间存在着极显著的正相关；同理可知xl与x4的相关系数R=0.99896,P=0.0001<0.01;x2与x3的相关系数R=0.99846,P=0.000K0.01；x2与x4的相关系数R=0.99854,P=0.0001<0.01；x3与x4的相关系数R=0.99999,P=0.OOOKO.01。由以上的分析可知，xl,x2,x3,x4之间的相关系数接近1,可见模型一的稳定行很强，而且由公式推出的前一两年的数据与该公园已有的数据基本相符合，可见模型是很优的。问题二：由于目前该动物已经很少被移出或移入，而且偷猎枪杀的情况微乎其微，所以暂时不列入考虑范围内。因此对该动物群落若不采用人工手段控制，则其在一定时间范围内会大幅度增加，从而破坏该种群的动态平衡。为了保持该种群的平衡，而又不必每年运走一定数量动物，现在使用-利避孕注射法，可使该动物群落的数量固定在一定范围内，用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕，但不会引起其它附加的反应。我们所要做的就是估计出每年在该群落中要注射避孕药的雌性动物的数量，并且要考虑到各种不确定性因素的影响。为了分析的方便，我们先建立初步模型，该模型暂时不考虑注射避孕药所产

10生的不确定性因素的影响，即不考虑两年内被重复注射的雌性数量及双胞胎的儿率。在这里我们只认为新生幼仔的数量由两部分组成，一部分为没注射过避孕药的雌性个体所生，另一部分为被注射过避孕药的雌性个体所生。另外，由于已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数£（该种群每年的新生幼仔的数量减去该年死亡个体的数量的差值），当£趋于0时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转化为求单目标的最优化问题。从而建立模型如下：模型二：minz=£-xo=°<月+为=G£=九％0-[（1-/2）C2+C3]其中：4=」-,表示每年没有注射避孕药的雌性动物生幼仔的几率。3.52,=—.表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔■5.5的儿率；这里设。=2596,表示从13〜60岁该动物的雌性个体的总数为2596；设G=10243,表示从1〜59岁该动物的个数总和为10243；设J=17,表示60岁该动物的个体总和为17；£：表示该种群每年的新生幼仔的数量与该年死亡个体的数量总和的差值。把已知的数据代入上述模型，从而得到以下模型：minz-£O.2857y1+0.1818y2-x0=0<%+乃=25690.75%=529+£利用Matlab软件进行编程（源程序见附录4）,求得

11%=2293,y2=276,%=705,£=2.0286x10~16.其中£=20286xl（ri6p0,说明该模型是稳定的，即该动物群落的数量被控制在一定的范围内。由丫2=276可知，每年大约有276头雌性动物要注射避孕药，才能使该群落的数量保持在11000头左右。模型三（模型二的改进）：由于每年被注射的雌性动物数量一定，所以被注射过后两年有可能又被注射，则将其归入到“被注射过避孕药但在两年内不再被注射的雌性动物”内。因此，注射避孕药所产生的不确定性因素之一即为“被注射过避孕药但在两年避孕期内被重复注射的雌性动物”。另外不确定因素之二即双胞胎的儿率问题，所以这里我们加入一个参数为。当第一年采用注射避孕药的方法时，是不会发生有雌性个体被重复注射的情况的，故有以下模型：minz=£’4y1+4y2+%(yi+y2)-x0=o<乂+%=《£=/]X0-[(1-/2)C2+C3]其中：4=」-,表示每年没有注射避孕药的雌性动物生幼仔的几率；3.5九=—,表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的〜5.5儿率；4=」-,表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的6.5儿率；/,=0.75,表示新生幼仔的存活率；%=0.95,表示1〜60岁动物的存活率；心=00135,表示新出生的幼仔中双胞胎的概率。minz=£’0.2857yl+0.1818y2+0.0135(0.2857yl-K).1818y2)-xo=O

12利用Matlab软件求解(源程序见附录5),求得月=2202,%=367,%=705,£=2.4791x10-16其中£=24791x1(116po,说明该模型是稳定的，即该动物群落的数量被控制在一定的范围内。由必=%7可知，每年大约有367头雌性动物要注射避孕药，才能使该群落的数量保持在11000头左右。由于考虑双胞胎的机率，故必须增加注射的数量。当注射避孕药一年后再次注射时，就会有某些数量的雌性个体被重复注射的情况出现，但这部分一定比前一年注射的雌性个体的数量少，故建立以下模型模型四：ninz=£电+4y2+4丫3+%（4"4丫2+勿3）-%4）V%+必+%=。£=%%一[（1一%）C2+C3]代入已知数据，得到以下目标函数模型：minz=£0.2896y]+0.1843y2+0.1698y3-xo=O丫3«公V%+%+为=2569£=0.75%-529其中。是前一年被注射避孕药的雌性个体的数量。利用Matlab软件编写程序求解，通过对。的不同取值进行调试，（源程序见附录6）求得当。=242时，模型稳定，此时%=2217,%二242,乃二110,%=705,£=1.6409x10*其中飞=705与模型三的解相同，说明新生幼仔的数量比较稳定，而每年注射避孕药的雌性动物则减少到352头，其中有110头是在注射后一年又被重复注射的数量，而£=1.6409X10-12«0

13则说明了该动物群落的个体数量是稳定的，被控制在一定范围内。由此可知当每年注射352头时，就能把该动物群落的数量控制在11000头范围内。问题三:假如每年转移50至300头此动物到别处，其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b）,即e应表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量的和的差值。相应的避孕措施将改变如下：模型五：minz=s‘4%+办丫2+4丫3+%（4y1+4y2+々丫3）-xo=°丫3（①£=^xQ-[（l-y2）c2+c3+b]其中b表示该年被运走的动物的数量（50〜300）；代入已有数据，并把方程标准化，得到以下模型：minz=£0.2896"1843y2位1698y3-x0=0y3<69〈X+%+%=△8—0.75%0—（0.05c9+

14CommandVWndow»Exiting:Oneormoreoftheresiduals,dualitygap,ortotalrelative二hasgrown100000timesgreaterthanitsminimumvaluesofar:theprimalappearstobeinfeasible(andthedualunbounded).(Thedualresidual

15式（3）和（4）表示该动物第i龄到了年底全部转化为（i+1）龄；式（5）和（6）表示该动物各年龄段的变化；式（7）表示该动物新生的幼仔数量；式（8）表示该动物在13〜60岁的雌性个体数量之和；对b以10为间距，分别从50取到300头，利用MATLAB软件进行求解（源程序见附录8）,得到以下表二：被运走的数量（头）5060708090100110120130140被注射的数量（头）370361351342332322313304294285被运走的数量（头）150160170180190200210220230240被注射的数量（头）275266256247237228218209199190被运走的数量（头）250260270280290300被注射的数量（头）180171161152142133问题四：用人工手段可以使得该动物群落的数量减少，使其数量控制在一定范围内，但是,我们并不知道用了这种人工手段后，对于一次突发事件后要增加其数量时,该种群的复壮能力如何。因此，我们建立数学模型，分析研究其复壮能力。我们假设由于某种原因（例如由于一场大灾难导致该动物的死亡），使得该动物减少至M头，此时不得不停止注射避孕的方法。在初步模型中，我们先假设灾难发生后该种群虽然数量减少，但其年龄结构没有改变。建立该动物群落重新壮大的能力的模型如下：模型七：

16dX⑹(八6060(1)ai/=iz=i60之硝=M⑵端+□=X；；).....(z=0,1,..59)(3)靖)=0(z=60)(4)琮=%瑞・・・・・・(/=0)(5)瑞=/瑞)..…a=1,2,…,59)(6)对M取不同的值，利用c语言进行编程求解，（源程序见附录9）可以得到以下表：第t年1000090008000700060005000400030002000t=l增量dl当前数量Ml6651066559495945268526456745638863883225322250425017931791132113t=2增量d2当前数量M270711372638102325629088489794511160823425664267451719133701222235t=3增量d3当前数量M3757恢复699109326029690521846644472463656029286480320535751302365t=4增量d4当前数量M4729恢复64010330557902347277183906419305510821937941402505t=5增量d5689596507418325235145

17当前数量M511019961982256837543340292650t=6增量d6741恢复6381025753887924457272348578125142801562806当前数量M6t=7增量d76821092958193734727744371615226845481662972当前数量M7t=8增量d8729恢复62399965078251401655328645481783150当前数量M8t=9增量d9665106615438794426697930748351913341当前数量M9t=10增量diO707113685819375456743532751422043545当前数量M10dx2(t]增长加速度：。=——udr46.000046.600042.685738.881035.642428.612122.806116.436410.7636该表对M进行不同的取值，分别表示某…年(t=0)灾难发生后该动物群落剩下的数量。我们用第•题的C程序即可算出该群落在下一年(t=l)的增量，把这个增量加上上一年群落的总数即为第一年(t=l)群落的总数Ml,接着我们又用程序算出这一年的增量，加上Ml即为第二年(t=2)的个体总数M2。以这样的算法推算下去，一直到群落数量恢复到11000左右或者前十年(t=10)为止。将计算出来的各年增量按时间t(1,2……10)列在座标上，用MATLAB进行曲线拟合(具体程序见附录10),得到的图形如下所示:

18从图中可以看出，每年的增量都是遵从一定的线性关系的，也就是说，该群落的增长加速度也坐是一定的，通过MATLAB的曲线拟合可以求出不同M值所对应的增长加速度（程序见附录）,数据如下:M值1000090008000700060005000400030002000D（只/46.0046.6042.6838.8835.6428.6122.8016.4310.76年八2）000057102421616436通过这些数据，我们可以知道，在灾难后群落数量减少至M时的重新壮大能力。再用MATLAB求出M与D的关系，进行曲线拟合求出斜率k后画出图形（程序见附录11）,如下所示55

19可知M与D是服从一定的线性关系的，由MATLAB算出的的数据可知M=3.9010+0.0047Do该模型的稳定性分析见本论文第四部分。四.模型的稳定性分析1.对于问题一：我们利用C语言进行编程，求得了第k年时各年龄结构的分布情况，并利用MATLAB软件对得到的数据进行画图分析，从输出的图可看出各年该动物各个年龄的分布是围绕一定的轨迹进行分布的，从而说明了该模型是稳定的，其次，用统计软件SAS对该模型得出的数据进行相关系数的分析，得到相关系数R都接近于1,而P=0.000KO.01则表示其相关性显著。故建立的模型一是很稳定的，其各年的该动物的总数基本上都在11000左右。2.对于问题二、问题三：我们从两个不同方面对问题进行分析，分别从£=新生的幼仔数量一自然死亡的数量〜0,以及£=该动物群落在年底时的总数量与年初的数量加上被运走的数量b的差值〜0两个不同方面进行模型的建立，从而把该动物群落的稳定性问题转化为求目标函数最小化问题。（这是本论文的一

20个亮点）。利用MATLAB对模型进行求解，得出的数据最后都能使得该动物群落维持稳定，可见对于问题二、问题三所建立的模型是很稳定的，结果也是很优的。1.对于问题四，我们用C语言的编程求得了不同M值在10年内各年的增量,求出该群落在不同M值所对应的增长加速度D„为了求出不同的M与对应的D的关系，我们在MATLAB中输入M与D的数组，进行稳定性分析（程序见附录11）。在回归分析及检验中，我们得到，片=3.9010,4=0.0047；4的置信区间为[-1.5746,9.3766],仇的置信区间为[0.0039,0.0055];r2=0.9616,F=175.1021,p=0.0000,p<0.05,可知回归模型M=3.9010+0.00470成立对模型进行残差分析，作残差图，得图如下：OResidualCaseOrderPlotCaseNumber从残差图可以看出，除最后一个数据外,其余数据的残差离零点均较近,且残差的置信区间均包括零点,这说明回归模型M=3.9010+0.00470能较好的符合原始数据，而最后一个数据视为异常点.（造成异常点的原因是M=10000时群落只要两年就可恢复到原来的数量11000,数据太少无法作出较接近实际的曲线拟合。

21在预测及作图后得各数据点及回归方程得图形如图:可以看出，只有最后一个数据点离回归直线距离较远是异常点。由以上分析可知M与D的关系是比较稳定的。五.模型优缺点及改进方向我们的模型有以下优点：1.模型的稳定性好，模型给出的方案使得该动物群落的各年龄段的数量保持平衡，从而满足了维持该动物群落的稳定。2.模型的适用范围广，易于推广，模型对于其它生态经济现象（如最优捕鱼策略问题）同样适用。3.基本模型对问题的描述准确、合理、推导严谨，理论性强；4.模型结合实际，具有很高的实用价值。我们的模型有以下缺点：1.由于题目给的数据不够多，所以使得该模型无法更加接近实际的情况。2.由于题目没有以前各年偷猎枪杀该动物的数据，故模型没有考虑偷猎枪杀该

22动物所产生的影响。1.对于…个动物群落，其影响的不确定性因素很多，而模型只是局限在题目中提到的，故与实际情况有•定的偏差。模型的改进方向：1.由生活常识可知，对于该动物群落超过60岁后的死亡情况应基本上服从正态分布。故可在模型-中的死亡数中加上这个因素。2.由于动物的产幼仔是随机分布的，故可利用模拟仿真的方法建立更具有真实性的模型。参考文献[1]樊欣，邵谦谦，SAS8.X经济统计，北京希望电子出版社，2003.3.[2]苏金明，张莲花，MATLAB工具箱应用，北京电子工业出版社，2004.L[3]姜启源，谢金星，数学模型（第三版），高等教育出版社，2003.8.[4]余世孝，数学生态学导论，科学科技文献出版社，1995.8.

23附录一：(模型一用c语言实现的过程与运行结果)main(){intx[7O],n;/*定义个长度为70的数组，每个单位的值表小各年龄的数量，定义变量n*/intsum,suml,sum2;/*个变量分别代收该年群落个体数量，13—59岁雌性数量，1—59岁个体数量*/clrscr();sum=0;suml=0;sum2=0;printfCzBorneveryyear」)；scanf(*%d*,&x[O]);/*输入节年出生的幼仔数量,设幼仔年龄为0岁*/printf(/z

24Age:quantityinthisage

25

26");x[l]=0.75*x[0];/*幼仔从0岁向1岁过渡时存活率为75%*/for(n=2;n〈=60;n++){x[n]=0.95*x[n-l];/*向下一年龄段过渡时的存活率为：删*/)for(n=0;n<=60;n++){sum=sum+x[n];/*对各年龄段数量累加得到种群总体数量*/printf(*%d:%d”,n,x[n]);)for(n=13;n<=60;n++){suml=suml+x[n];/*对13—60岁个体数量累加*./)for(n=l;n<=59;n++)

27{sum2=sum2+x[n];/*对1―59岁个体数量累加*/printf("

28

29sum=%d

30”,sum);

31printf(

32sum(13―59female)=%d

33,suml/7);/*13―60岁个体总数的一printf("

34sum(l―59)=%d

35”,sum2);半即为该年龄段雌性数增再除以3.5即为当年生幼仔的雌性的数量*/通过不断输入新生幼仔的初值,直到算出当年群落的个体总数与题设数据吻合,则该组数据即为当年群落年龄结构的分布(未算入当年运走的数量).运行结果：1.前雨年运走的个体数量为876头，因此前一年的个体总数为11000+876=11876(头)，通过调试程序使总数sum达到11876左右.第一次调试结果：zquantituinthisageeuei*y731674672:84203:88--92777:0134551999:345251382152:14243:9842nu2:::612888:13453:0610:137420008456426115111245866348122:58145501:4463497133:81:4562:302863260190534461::471:44••2--4457:538:41817:26026:16035:9645:5355:2816510674195292::::3--7666:8234591ur<13——60fenale>=7922.前两年年末即前一年年头所剩下的13-60岁的雌性个体数量决定了前一年年末个体的数量即2470/3.5=706头；前两年剩下的数量加上前一年新生幼仔的数量11000+706=11706,其值为前一年年末的个体数量用程序算得11714,误差为8头，其误差由于变量是整型变量，数字取整数所致。第二次调试结果：

36Bomeueryyeai*=792Age：quantityinthisage9:39118:24327:15036:9046:5056:268538265952151:::42:555----63458727:43416:27025:16734:10144:5654:30927753^0^00••••457133:81:4562:31342630000••••46122:181454276816305011:6245479153273131:::5000:74563723:05104:13742622:::51999--34521143:99:0251210428428883451:7920:37119:23028:14237:8547:4757:24fenale>=784在前一年的年末运走622头，即该群落的个体数量剩下11714-622=11092根据表得到2503/3.5=715,所以11092+715=11807为该群落今年年尾的个体总数。euery：quantitvinthisage431674072:84203:8:::8--92777:01345211:35520:22029:13538:8148:4458:22：5693:5404:5125：4866：46112:33713:32014:30315:28721:20822:19723：18724:17730:12831:12132:11433：10839:7640:7241：6842:6443：6049:4150:3851:3652：3453：3259:2060：187:43716:27225:16834:10244:5654:30872115142:55534596385859529:39418:24527:15136:9146:5056:26feRale>=789一一59)=10990程序计算得11808,误差为1,其误差由于变量是整型变量，数字取整数所致。1.对于以上数据是否符合假设，我们还要用一-组初步数据来进行对比.因此我们设每年运出的个体数量为0,个体总数为11000头,通过调试程序使总数sum达到11000左右.

37第三次调试结果:euery：quantitvinthisage82147462:8423:7::::826669134581490684185232::655587234596438045952451:::--2--4447:53456211060350534861:::61:3362:345926530•«•・:37122581:453463767631800014561:1318832273151:::999534593:95105:027425128834263911642:94243:87427:92:01777fenale>=734该结果为当每年没有个体被运走时数量停留在11000左右的数据，用于与初步数据进行对比。附录2（用MATLAB画出模型一轨迹）：clearall;clcyl=[7925945645355084824574344123913713523343173012852702562432302152071961861761671581501421341271201131071019590858075716763595653504744413836343230282624222018]';y2=[8066045735445164904654414183973773583403223052892742602462332212091981881781691601511431351281211141081029691868176726864605653504744413836343230282624222018]';y3=[746559531504478454431409388368349331314298283268254241228216205194184174165156148140132125118112106100948984797571676359565350474441383634323028262422201817]';y4=[80060056954051248646143741539437435533732030328727225824523222020819718717716815915114313512812111410810296918681767268646056535047444138363432302826242220181）;

38xl=[O00034720915922323513210221451310013301412020635812103714101621131012636913103621154131032140]);x2=[00202113121322144014261314273141220251714100234432313161310101216121012191324171625124523341316101713131232220]';ml=yl-xl;m2=y2-x2;m3=y3;m4=y4;t=0:60plot(t,ml,"b.'),holdon,plot(t,m2,'y*'),holdon,plot(t,m3,"g-')，holdon,plot(t,m4,'m')legendC前一年剩余数量的年龄结构情况'，'前两年剩余数量的年龄结构情况没有被运走时的年龄结构情况'，’预测当前的年龄结构情况(无运走)’)titleC今年时该动物群落的年龄结构情况')附录3（用SAS对模型一进行检验的程序及调试结果）:dataA;inputxlx2x3x4cards;594604564553535523505503478478450452414419403404376357362363330332331327294308296278559600531569504540478512454486431461409437388415368394349374331355314337298320283303

392722862682872492602542722562482412582212262282452162082162322132042052201941951942081861881841971861881741871631761651771371661561681441561481591381471401511421401321431141331251351211251181281171081121211089810611499951001088992941028586899687798491787079866669758165667176556067724649636850515964493656604439535647375053472547504135444735-1414428183841284363833233436281832341122303215132830241526281313242614122224-10192022

406-2182018-217187proccorr;varxlx2x3x4;run;输出结果如下：TheSASSystem21:39Saturday,April24,20051TheCORRProcedureVariables：xlx2x4SimpleStatisticsVariableNMeanStdDevSumMinimumMaximumxl60171.66667162.2352010300-10.00000594.00000x260169.90000163.7510210194-2.00000604.00000x360171.00000151.456991026017.00000559.00000x460183.46667162.157331100818.00000600.00000PearsonCorrelationCoefficients,N=60Prob>|rlunderHO：Rho=0xlx2x3x4x11.000000.997440.998960.99896<.0001<.0001<.0001x20.987441.000000.998460.89854<.0001<.0001<.0001x30.998960.998481.000000.99999<.0001<.0001<.0001x40.998960.998540.999991.00000<.0001<.0001<.0001附录4(用MATLAB求解目标最优化问题)clearall;clcf=[0;0;0;l];Aeq=[0.28570.1818-10;l100;000.75-1];beq=[02569529];lb=zeros(4,1);[x,fval,exitflag,output,lambda]=1inprog(f,[],[],Aeq,beq,lb)输出结果如下：

41>>Optimizationterminatedsuccessfully.1.0e+003*2.29340.27560.70530.0000fval=%返回目标函数值%2.0286e-016exitflag=%表示求得的解是收敛的％output=iterations:cgiterations:algorithm:lipsol%求目标最优化的方法％lambda=%所占的空间％ineqlin:[0x1double]eqlin:[3x1double]upper:[4x1double]lower:[4x1double]附录5(用MATLAB求解模型三目标最优化问题)clearall;clcf=[0;0;0;1];Aeq=[0.28960.1843-10;l100;000.75-1];beq=[02569529];lb=zeros(4,1);[x,fval,exitflag,output,lambda]=1inprog(f,[],[],Aeq,beq,lb)输出结果如下：>Optimizationterminatedsuccessfully.

421.0e+003*2.20200.36700.70530.0000fval=2.4791e-016exitflag=1output=iterations:6cgiterations:0algorithm:'lipsoV附录6:(用MATLAB求解模型四目标最优化问题)clearall;clcf=[0;0;0;0;l];A=[00100];b=[242];Aeq=[0.28960.18430.1698-10;l1100;0000.75-1];beq=[02569529];lb=zeros(5,1);[x,fval,exitflag,output,lambda]=1inprog(f,A,b,Aeq,beq,lb)输出结果如下：>>Optimizationterminatedsuccessfully.1.0e+003*2.21710.2421%由此可知当勿=242时模型稳定％

43o.1098p.7053q.0000fval=1.6409e-012exitflag二1output=iterations:6cgiterations:0algorithm:'lipsoTlambda=ineqlin:1.9145e-014eqlin:[3x1double]upper:[5x1double]lower:[5x1double]附录7(用MATLAB求解模型五目标最优化问题)clearall;clcf=[o；o；o；i]；Aeq=[0.28960.1843-10;l100;000.75-1];beq=[02569579];lb=zeros(4,1);[x,fval,exitflag,output,lambda]=1inprog(f,口，[],Aeq,beq,lb)输出结果如下：

44>>Exiting:Oneormoreoftheresiduals,dualitygap,ortotalrelativeerrorhasgrown100000timesgreaterthanitsminimumvaluesofar:theprimalappearstobeinfeasible(andthedualunbounded).(Thedualresidual

45370361351342332323313304294285275266256247228218209199190180171161152142133附录9(用C语言求不同M值时各年的增量)main(){intx[70],n,d;intsum,suml,sum2;clrscr();sum=0;sum1=0;sum2=0;printf(HBomeveryyear=K);scanf("%dM,&x[0]);printf("

46Age:quantityinthisage

47

48H);x[l]=0.75*x[0];for(n=2;n<=60;n++){x[n]=0.95*x[n-1];

49{sum=sum+x[n];printf(M%d:%dM,n,x[n]);/*打印出各年龄段的数量*/)for(n=13;nv=60;n++){suml=suml+x[n];)printf("

50

51sum=%d

52",sum);/*用总和来捽制，使之等于M在第t年数量*/printf("

53sum(increase)=%d

54",sum1/7);/*每年的增鼠约等于13—16岁雌性个体的数量*/)附录10(求不同M值增量d的关系)clearall;clct=[l2345678910];tl=[l2];t2=[l23];t5=[l234]t3=[l23456];t4=[l2345678];Ml=[665707757];M2=[594638699729];M3=[526562602640689741];M4=[456489521557596638682729];M5=[388414444472507538581623665707];M6=[322342365390418445472507543581];M7=[250267286305325348371401426456];M8=[179191205219235251268286307327];M9=[113124134141152161173185199212];

55f=inline(，k(l)*t+k(2)','k','t');kl=[O.530];ml=lsqcurvefit(f,kl,t2,Ml)zl=f(ml,t2)k2=[0.530];m2=lsqcurvefit(f,k2,t5,M2)z2=f(m2,t5)k3=[0.530];m3=lsqcurvefit(f,k3,t3,M3)z3=f(m3,t3)k4=[0.530];m4=lsqcurvefit(f,k4,t4,M4)z4=f(m4,t4)k5=[0.530];m5=lsqcurvefit(f,k5,t,M5)z5=f(m5,t)k6=[0.530];m6=lsqcurvefit(f,k6,t,M6)z6=f(m6,t)k7=[0.530];m7=lsqcurvefit(f,k7,t,M7)z7=f(m7,t)k8=[0.530];m8=lsqcurvefit(f,k8,t,M8)z8=f(m8,t)k9=[0.530];m9=lsqcurvefit(f,k9,t,M9)z9=f(m9,t)plot(t2,Ml,'b.'),holdon,plot(t2,zl,'r」),plot(t5,M2,'b.'),holdon,plot(t5,z2,'r--'),plot(t3,M3,'b.'),holdon,plot(t3,z3,'y，),plot(t4,M4,'b.'),hold

56on,plot(t4,z4,'y-'),plot(t,M5,'b.'),holdon,plot(t,z5,'g」),holdon,plot(t,M6,'b.'),holdon,plot(t,z6,'g一，),plot(t,M7,'b.'),holdon,plot(t,z7,',plot(t,M8,'b.'),holdon,plot(t,z8,'m-->),holdon,plot(t,M9,‘b.'),holdon,plot(t,z9,'k」)titleCM取不同的值时每年增量的分布情况')其运行结果如下：Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunml=46.0000617.6667zl=663.6667709.6667755.6667Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm2=46.6000548.5000z2=595.1000641.7000688.3000734.9000Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm3=42.6857477.2667z3=519.9524562.6381605.3238648.0095690.6952733.3810Optimizationterminatedsuccessfully:NormofthecurrentstepislessthanOPTIONS.TolXm4=38.8810408.5357z4=Columns1through7447.4167486.2976525.1786564.0595602.9405641.8214680.7024Column8719.5833Optimizationterminatedsuccessfully:

57NormofthecurrentstepislessthanOPTIONS.TolXm5=35.6424337.8667z5=Columns1through7373.5091409.1515444.7939480.4364516.0788551.7212587.3636Columns8through10623.0061658.6485694.2909Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm6=28.6121281.1333z6=Columns1through7309.7455338.3576366.9697395.5818424.1939452.8061481.4182Columns8through10510.0303538.6424567.2545Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm7=22.8061218.0667z7=Columns1through7240.8727263.6788286.4848309.2909332.0970354.9030377.7091Columns8through10400.5152423.3212446.1273Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm8=16.4364156.4000z8=Columns1through7172.8364189.2727205.7091222.1455238.5818255.0182

58271.4545Columns8through10287.8909304.3273320.7636Optimizationterminatedsuccessfully:First-orderoptimalitylessthanOPTIONS.TolFun,andnonegative/zerocurvaturedetectedm9=10.7636100.2000z9=Columns1through7110.9636121.7273132.4909143.2545154.0182164.7818175.5455Columns8through10186.3091197.0727207.8364附录11(M、D关系的稳定性检验)x=[2000300040005000600070008000900010000];y=[10.763616.436422.806128.612135.642437.535740.400046.600046.0000];f=inline('k(l)*x+k(2)','k','x');k=[0.530];m=lsqcurvefit(f,k,x,y)z=f(m,x)plot(x,y,'.'),holdon,plot(x,z,')其输出结果如下：»Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunm=0.00464.0893z=Columns1through7

5913.274217.866622.459127.051631.644036.236540.8290Columns8through945.421450.0139附录12(问题四稳定性分析)»x=[2000300040005000600070008000900010000]';>>X二[ones⑼1)x];Warning:FutureversionsofMATLABwillrequirethatwhitespace,acomma,orasemicolonseparateelementsofamatrix.Pleasetype“helpmatrix_element_separators/zattheMATLABpromptformoreinformation.»Y=[10.763616.436422.806128.612135.642438.881042.685746.600046.0000]);>>[b,bint,r,rint,stats]=regress(Y,X);»b,bint,statsb=3.90100.0047bint=-1.57469.37660.00390.0055stats=0.9616175.10210.0000>>rcoplot(r,rint)»z=b⑴+b(2)*xz=13.283217.974322.665327.356432.047536.7386

6041.429646.120750.8118>>plot(x,Y〉k+',x,z,'r')