动物群落的稳定发展

论文评阅要点

一、主要标准： 1、 假设的合理性； 2、 建模的创造性； 3、 文字表达的清晰性； 4、 结果的正确性。 二、论文组成概要： 1、 题目 2、 摘要 3、 问题重述

4、 模型假设与符号 5、 分析建立模型 6、 模型求解

7、 模型检验与推广 8、 参考文献与附录

三、参考给分步骤（10分制）

1、 摘要部分（论文的方法、结果、表达饿清晰度）。。。。。。。。。。。。。。3分 2、 假设部分（合理性与创造性）。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。1分 3、 数学模型（创造性与完整性）。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。3分 4、 解题方法与结果（创造性与正确性）。。。。。。。。。。。。。。。。。。。。。。。。2分 5、 模型的优缺点与推广（合理性）。。。。。。。。。。。。。。。。。。。。。。。。。。。。1分 四、评阅方法

1、 每位教师把卷号、分数及主要理由记录在白纸上，以便专人统计；

2、 每份论文至少要三位教师评阅过，选出获奖论文的2倍数量，对分歧大的试卷讨论给分； 3、 对入选论文至少要六位教师评阅过。按分数高低排序；

4、 对一、二等奖的论文要求写出30字左右的评语，与论文一起在网上发表。 五、评阅时间：5月21日（星期六）

A题：动物群落的稳定发展

摘要：本文通过对某公园近两年内被运出的某种动物的年龄和性别的数据进行统计分析，并针对题目的四个问题分别建立了符合实际的数学模型，在模型的求解过程中，应用C语言进行编程调试，通过统计学软件SAS，数学软件MATLAB等计算工具，编写相应的程序，对建立的模型进行求解，得出了符合实际的结果。

问题一：我们假设新生幼仔的数量为x0，然后通过对各年龄阶段的存活率?、被运走的动物数量Bj以及该动物的总体数量的分析来建立该群落的动态变化模型

d?(k)?t?60(k)60(k)??xi,1??xi,0，利用该群落近两年内被运走的各年龄阶段的个体数量分布，dti?1i?1用C语言编程计算，推测出当前该动物的年龄结构（具体结果见 7 页表一）。并利用MATLAB软件对得出的数据用图形表示，利用对比分析法，得到该动物群落的基本分布轨迹，最后用统计软件SAS对模型进行相关性的分析检验，求得相关系数R与P的值，验正了模型的稳定性。

问题二：由于现在采用注射避孕药的方法来维持该种群的稳定，而且已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数

?1??1x0?[(1??2)c2?c3]（该种群每年的新生幼仔的数量减去该年死亡个体的数量的差

值）；另外从???xi?160(k)i,1k)（即年头的数量与该年年底的数量的差值）当?趋于0??xi(,0i?160时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转化为求单目标的最

优化问题建立模型；利用MATLAB对模型进行求得，得出当不考虑不确定性因素影响时要注射药物的雌性动物数量为276头，而当考虑了双胞胎和被重复注射这两个不确定性因素影响后，得到要注射药物的雌性动物数量为352头，其中有110头是被重复注射的。

问题三：其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b），即目标函数?应考虑上被运走的数量，即只是对问题二的模型进行扩充建立新的目标模型；?1??1x0?[(1??2)c2?c3]－b和???xi?160(k)i,1k)－b；利用MATLAB对不同b??xi(,0i?160值进行求解，从而得出相应的避孕措施。（具体结果见 19 页表二）

问题四：我们引进了增量加速度的概念，利用c语言进行编程求解，然后用MATLAB软件对得到的数据进行线性回归分析，得到该群落在减少至M时重新壮大该动物群落能力的模型：M＝3.9010+0.0047D。最后应用统计软件SAS对模型进行稳定性分析。

关键字：存活率 年龄结构 新生幼仔数 稳定性 最优目标 增量加速度

一． 问题重述与提出

位于非洲某国的国家公园中栖息着近11000头某种野生动物。管理员要求有一个健康稳定的环境以便维持这个11000头该动物的稳定群落。过去的20年中，整个该动物群是通过一些偷猎枪杀以及转移到外地而稳定下来的。但是近年来，偷猎被禁止，而且每年要转移这些动物也比较困难，因此，要控制现在的数量就使用了一种避孕注射法。用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕。

要探讨这种避孕注射法的实用性，我们需要完成以下问题：

1．探讨该动物年龄在2岁到60岁之间的合理的存活率的模型，推测这个动物群落的当前的年龄结构。

2．估计每年在该群落中有多少雌性动物要注射避孕药，可以式群落固定在11000头左右。这里不免有些不确定性，也要估计这种不确定性的影响。

3．假如每年转移50至300头此动物到别处，那么上面的避孕措施将可以有怎样的改变？

4．如果由于某种原因，突然使得注射避孕的方法不得不停止（例如由于一场灾难导致大量该动物的死亡），那时重新壮大该动物群的能力如何？

二． 基本假设与符号说明 （一）模型假设

1．该公园是非开放式的，它与外界不发生关系，从而构成独立的生物群落，该动物群落不存在与其它动物种群的竞争，或虽有竞争，但其影响只局限于该动物群落的死亡率内。

2．种群是通过雌性个体的繁殖而增长的，所以用雌性个体数量的变化为主要研究对象。

3．为了讨论的必要，我们把新生的幼儿的存活率定为75％，而其后的存活率为95％，直到60岁为止。各年龄组的该动物经过一年后即进入高一级的年龄组，而龄超过60即认为全部死亡，退出该系统。

4．由于该公园加强了对该动物群落的保护，我们认为该动物没有再被偷猎射杀。而该动物群落个体数量的减少只是因为自然死亡以及被运走。

5．假设同一年龄组的动物个体之间是同质的，我们只考虑其平均水平，不讨论个别

差异。

6．题设该动物在10～12岁开始怀孕，我们这里设定为11岁开始，经过22个月（约两年）的怀孕期后生幼仔，即可认为该雌性动物在13～60岁的时间内可以生幼仔。

7．该群落的自然死亡是在生完幼仔后才发生的，产幼仔只发生在每年的年初时段，而被运走只发生在年底时段。

（二）符号说明

?1：新生幼儿的存活率，其值为0.75；

?2：1～60岁个体的存活率，其值为0.95； ?3：双胞胎出生的几率，其值为0.0135；

?(k)(t)：该动物第k年时刻的数量；

k)

：该动物第k年初i龄动物的数量； xi(,0

k)

：该动物第k年初底i龄动物的数量 xi(,1

Bj：第j年被运走的动物的数量；

?(k)?0?：表示该动物第k年初时的总数量；

?1:表示每年没有注射避孕药的雌性动物生幼仔的几率，其值为

1； 3.51； 5.51表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的几率，其值为； ?3：

6.5表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的几率，其值为?2：

c1：表示从13～60岁该动物的雌性个体的总数；

c2：表示从1～59岁该动物的个数总和；

c3：表示60岁该动物的个体总和；

y1：表示13～60岁雌性动物没有被注射避孕药部分的数量；

y2： 表示13～60岁雌性动物被注射过避孕药但在两年内不再被注射部分的数量； y3：表示13～60岁雌性动物被注射过避孕药但在两年内被重复注射部分的数量；

?1：表示每年出生幼仔的数量与该年个体死亡的数量的差值；

表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量的和的?2：差值；

?3：表示该动物群落在年底时的总数量与年初的数量加上被运走的个体数量b的差值。

三． 问题分析与模型建立 问题一：

1．我们要研究该动物群落的稳定性问题，首先要根据存活率确定其当前的年龄结构。该动物的新生幼仔存活率较低，题设是70％到80％之间，为了讨论的需要，我们这里设定为75％。在1岁后的存活率比较高，在这里设为95％，直到60岁，而超过60岁则认为退出该系统。因此，我们先建立出该动物群落中年龄在2岁到60岁之间的合理的存活率的模型。

?d?(k)?t?60(k)60(k)??xi,1??xi,0????????????(1)?dti?1i?1??60(k)(k)x??(0)??????????????????(2)??i,0?i?1?d?(k)?t??Bj???????????????????(3)dt??x(k?1)?x(k)......(i?0,1,..59)???????????(4)i,1?i,0?x(k)?0............(i?60)??????????????(5)?i,1模型一：?x(k)??x(k)......(i?0)??????????????(6)

i,11i,0?k)k)?xi(,1??2xi(,0......(i?1,2,...,59)??????????(7)?60?(k)x?i,0?(k)?????????????????(8)?x0,0??1i?13?2

式（1）表示该动物第k年增长的数量；

式（2）表示该动物第k年初时的总数量，可由已有的数据计算出?(k)(0)来； 式（3）表示该动物被运走的数量；

式（4）和（5）表示该动物第i龄到了年底全部转化为（i＋1）龄； 式（6）和（7）表示该动物各年龄段的变化； 式（8）表示该动物新生的幼仔数量。

2．通过对该公园近两年内从这个地区运出的该动物的年龄和性别的数据进行统计分析，并利用编程工具Turbo C 2.0对该模型进行编程计算（源程序及计算过程见附录1），可得到当前该动物群落的年龄结构，如下表所示：

表一 该动物的年龄结构统计表

前一年数量前一年运走前一年剩下前兩年数前兩年运走前两年剩下假设无运走今年数年龄（岁） （头） 数量（头） 数量（头） 量（头） 数量（头） 数量（头） 数量（头） 量（头） 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 792 594 564 535 508 482 457 434 412 391 371 352 334 317 301 285 0 0 0 3 4 7 20 9 15 9 22 3 23 5 13 786 594 564 535 505 478 450 414 403 376 362 330 331 294 296 272 806 604 573 544 516 490 465 441 418 397 377 358 340 322 305 289 0 20 21 13 12 13 22 14 40 14 26 13 14 27 3 806 604 553 523 503 478 452 419 404 357 363 332 327 308 278 286 746 559 531 504 478 454 431 409 388 368 349 331 314 298 283 268 800 600 569 540 512 486 461 437 415 394 374 355 337 320 303 287 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

270 21 256 0 243 22 230 14 218 5 207 13 196 10 186 0 176 13 167 30 158 14 150 12 142 0 134 20 127 6 120 3 113 5 107 8 101 12 95 10 90 3 85 7 80 14 75 10 71 16 67 21 63 13 59

10

249 256 221 216 213 194 186 186 163 137 144 138 142 114 121 117 108 99 89 85 87 78 66 65 55 46 50 49

274 14 260 12 246 20 233 25 221 17 209 14 198 10 188 0 178 2 169 3 160 4 151 4 143 3 135 2 128 3 121 13 114 16 108 13 102 10 96 10 91 12 86 16 81 12 76 10 72 12 68 19 64 13 60

24

250 248 226 208 204 195 188 188 176 166 156 147 140 133 125 108 98 95 92 86 79 70 69 66 60 49 51 36

254 241 228 216 205 194 184 174 165 156 148 140 132 125 118 112 106 100 94 89 84 79 75 71 67 63 59 56

272 258 245 232 220 208 197 187 177 168 159 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 生幼仔的雌性数量 13—60岁雌性数量 总数量 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 12 6 3 6 9 13 10 3 6 21 15 4 13 10 32 14 0 44 47 47 41 35 28 28 33 28 11 15 24 13 14 -10 6 18 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 17 16 25 12 45 23 34 13 16 10 17 13 13 12 3 22 20 39 37 25 35 -1 18 4 23 18 22 13 15 13 12 19 -2 -2 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 17 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 784 792 734 789 2735 11714 232 622 2503 11092 2772 11876 302 876 2470 11000 2569 11006 11808 注1：0岁表示新生幼仔。

注2：由于每个年龄段的数据均为推测值，而实际上运走的各年龄段的数量不一定全部与预测值相符，故表中“剩下数量”两组数据中出现负数可认为是独异点，不影响模型整体的准确性。

每年新生幼仔的数量（x0）减去生幼仔的雌性的数量（

?x/2i13603.5，由于雄性与雌性

的数量比接近1：1，我们可近似地认为13——60岁个体的雌雄数量相等），其差值即为双胞胎的数量，这个差值与生幼仔的雌性数量之比即为双胞胎的几率（?3）1.35％。由表中数据可得，

792?784806?792746?734?0.0102，?0.0177，?0.0163，784792735800?789?0.0139，这些比例都基本上接近题设的双胞胎的几率0.0135，说明以上推测789得出的数据是准确的。

利用Matlab软件对以上四组数据用图形表示，并进行比较，得到该动物群落的基本分布情况图（源程序见附录2），如下图所示

图1

分析该图，可以看出，这四组曲线的轨迹、分布情况基本相同。由于“预测当前的年龄结构情况（无运走）”一组数据没有减去被运走的个体数量，故其每个年龄层的数量都略多于前三组的数量，因此其曲线比前三组的曲线略高一点，

利用SAS软件对模型进行相关性的分析检验（源程序见附录3），得到如下结果：

图2

程序的分析及统计结论：

程序中的x1是前一年的该动物群落的年龄结构，x2是前两年该动物群落的年龄结构，x3是该动物群落没有被运出是的年龄结构，x4是预测的当前的该动物群落的年龄结构。过程中的PROC CORR是分析变量中两两变量之间的PEAROS简单相关的。

输出结果中的结果1是一些基本的描述统计量，结果2是两两变量之间的相关矩阵，其中包括相关系数和显著性检验的概率。由结果可知前一年的该动物群落的年龄结构（x1）与前两年的该动物群落的年龄结构（x2）的相关系数R=0.99744,P=0.0001<0.01，所以前一年的该动物群落的年龄结构（x1）与前两年的该动物群落的年龄结构（x2）之间存在着极显著的正相关；前一年的该动物群落的年龄结构（x1）与没有运走是的该动物群落的年龄结构（x3）的相关系数R=0.99896,P=0.0001<0.01，所以前一年的该动物群落的年龄结构（x1）与没有运走时的该动物群落的年龄结构（x3）之间存在着极显著的正相关；同理可知x1与x4的相关系数R=0.99896,P=0.0001<0.01;x2与x3的相关系数R=0.99846,P=0.0001<0.01；x2与x4的相关系数R=0.99854,P=0.0001<0.01；x3与x4的相关系数R=0.99999,P=0.0001<0.01 。

由以上的分析可知，x1,x2,x3,x4之间的相关系数接近1，可见模型一的稳定行很

强，而且由公式推出的前一两年的数据与该公园已有的数据基本相符合，可见模型是很优的。 问题二：

由于目前该动物已经很少被移出或移入，而且偷猎枪杀的情况微乎其微，所以暂时不列入考虑范围内。因此对该动物群落若不采用人工手段控制，则其在一定时间范围内会大幅度增加，从而破坏该种群的动态平衡。为了保持该种群的平衡，而又不必每年运走一定数量动物，现在使用一种避孕注射法，可使该动物群落的数量固定在一定范围内，用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕，但不会引起其它附加的反应。我们所要做的就是估计出每年在该群落中要注射避孕药的雌性动物的数量，并且要考虑到各种不确定性因素的影响。

为了分析的方便，我们先建立初步模型，该模型暂时不考虑注射避孕药所产生的不确定性因素的影响，即不考虑两年内被重复注射的雌性数量及双胞胎的几率。在这里我们只认为新生幼仔的数量由两部分组成，一部分为没注射过避孕药的雌性个体所生，另一部分为被注射过避孕药的雌性个体所生。另外，由于已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数?（该种群每年的新生幼仔的数量减去该年死亡个体的数量的差值），当?趋于0时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转化为求单目标的最优化问题。从而建立模型如下: 模型二：

minz???1y1+?2y2?x0=0?? ? y1?y2?c1????x?[(1??)c?c]10223?其中：?1＝

1，表示每年没有注射避孕药的雌性动物生幼仔的几率。 3.51?2?，表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的几率；

5.5这里设c1＝2596，表示从13～60岁该动物的雌性个体的总数为2596； 设c2＝10243，表示从1～59岁该动物的个数总和为10243； 设c3＝17，表示60岁该动物的个体总和为17；

?：表示该种群每年的新生幼仔的数量与该年死亡个体的数量总和的差值。 把已知的数据代入上述模型，从而得到以下模型：

minz???0.2857y1+0.1818y2?x0=0?y1?y2?2569 ??0.75x0?529???利用Matlab软件进行编程（源程序见附录4），求得

y1?2293,y2?276,x0?705,??2.0286?10?16.

其中??2.0286?10?16?0，说明该模型是稳定的，即该动物群落的数量被控制在

一定的范围内。由y2?276可知，每年大约有276头雌性动物要注射避孕药，才能使该群落的数量保持在11000头左右。

模型三（模型二的改进）：

由于每年被注射的雌性动物数量一定，所以被注射过后两年有可能又被注射，则将其归入到“被注射过避孕药但在两年内不再被注射的雌性动物”内。因此，注射避孕药所产生的不确定性因素之一即为“被注射过避孕药但在两年避孕期内被重复注射的雌性动物”。另外不确定因素之二即双胞胎的几率问题，所以这里我们加入一个参数?3。

当第一年采用注射避孕药的方法时，是不会发生有雌性个体被重复注射的情况的，故有以下模型：

minz????1y1+?2y2+?3(y1+y2)?x0=0?y1?y2?c1 ?????x?[(1??)c?c]10223?其中：

1，表示每年没有注射避孕药的雌性动物生幼仔的几率； 3.51?2?，表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的几率；

5.5?1＝

?3?1，表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的几率； 6.5?1?0.75，表示新生幼仔的存活率； ?2＝0.95, 表示1～60岁动物的存活率；

?3?0.0135, 表示新出生的幼仔中双胞胎的概率。

minz???0.2857y1+0.1818y2+0.0135(0.2857y1+0.1818y2)?x0=0?y1?y2?2569 ????0.75x0?[(1?0.95)10243?17]?

利用Matlab软件求解（源程序见附录5），求得

y1?2202,y2?367,x0?705,??2.4791?10?16

其中??2.4791?10?16?0，说明该模型是稳定的，即该动物群落的数量被控制在

?367可知，每年大约有367头雌性动物要注射避孕药，才能使该

一定的范围内。由y2群落的数量保持在11000头左右。由于考虑双胞胎的机率，故必须增加注射的数量。

当注射避孕药一年后再次注射时，就会有某些数量的雌性个体被重复注射的情况出现，但这部分一定比前一年注射的雌性个体的数量少，故建立以下模型

模型四：

minz????1y1+?2y2+?3y3+?3(?1y1+?2y2+?3y3)?x0=0?y3???? y1?y2?y3?c1?????1x0?[(1??2)c2?c3]?代入已知数据，得到以下目标函数模型：

minz???0.2896y1+0.1843y2+0.1698y3?x0=0?y3????

y?y?y?2569123????0.75x0?529?其中?是前一年被注射避孕药的雌性个体的数量。利用Matlab软件编写程序求解，通过对?的不同取值进行调试，（源程序见附录6）求得当?＝242时，模型稳定，此时

y1?2217,y2?242,y3?110,x0?705,??1.6409?10?12

其中x0?705与模型三的解相同，说明新生幼仔的数量比较稳定，而每年注射避孕

药的雌性动物则减少到352头，其中有110头是在注射后一年又被重复注射的数量，而

??1.6409?10?12?0则说明了该动物群落的个体数量是稳定的，被控制在一定

范围内。由此可知当每年注射352头时，就能把该动物群落的数量控制在11000头范围内。 问题三:

假如每年转移50至300头此动物到别处，其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b），即?应表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量的和的差值。相应的避孕措施将改变如下： 模型五：

minz????1y1+?2y2+?3y3+?3(?1y1+?2y2+?3y3)?x0=0?y3????

y?y?y?c1231?????1x0?[(1??2)c2?c3?b]?其中b表示该年被运走的动物的数量（50～300）； 代入已有数据，并把方程标准化，得到以下模型：

minz???0.2896y1+0.1843y2+0.1698y3?x0=0?y3????y1?y2?y3?c1

????0.75x0?(0.05c2?c3)?b?其中b的值在50～300之间变化，我们以50头为间距（即b的值分为50，100，150，200，250，300六种情况），利用软件MATLAB6.0计算得到每年转移从50至300头到别处时所应采取的避孕措施（源程序见附录7），当取b＝50时，得到以下输出结果： ：

由输出结果的提示可知，在100000次迭代后仍然无法找到最优值；由exitflag=-1可知，以上模型得到的数据是发散的，从而说明上述模型本身并不存在问题，只是计算机无法在有限的迭代次数内找到最优值。为此解决该问题，我们对以上模型进行修改，把目标函数?3定义为该动物群落在年底时的总数量与年初的数量加上被运走的数量b的差值，从而得到以下模型： 模型六：

minz??36060?(k)(k)??3??xi,1??xi,0?b－－－－－－－－（1）i?1i?1??60(k)??xi,0?11000－－－－－－－－－－－（2）?i?1k?1)k)?xi(,0?xi(,1......(i?0,1,..59)－－－－－－（3）?(k)??xi,1?0............(i?60)－－－－－－－－（4） ?(k)(k)x??x?i,11i,0......(i?0)－－－－－－－－（5）?(k)(k)x??x,2,...,59)－－－－－（6）2i,0......(i?1?i,1?x(k)??y(k)??y(k)－－－－－－－－－（7）1122?0,060?k)xi(,0??(k)(k)?y1?y2?i?13－－－－－－－－－（8）??2其中

式（1）表示目标函数，其值趋于0；

式（2）表示该动物在第k年初的数量，其值取11000；

式（3）和（4）表示该动物第i龄到了年底全部转化为（i＋1）龄； 式（5）和（6）表示该动物各年龄段的变化； 式（7）表示该动物新生的幼仔数量；

式（8）表示该动物在13～60岁的雌性个体数量之和；

对b以10为间距，分别从50取到300头，利用MATLAB软件进行求解（源程序见附录8），得到以下表二：

被运走的数量（头） 被注射的数量（头） 50 60 70 80 90 100 110 120 130 140 370 361 351 342 332 322 313 304 294 285

被运走的数量（头） 被注射的数量（头） 150 160 170 180 190 200 210 220 230 240 275 266 256 247 237 228 218 209 199 190

被运走的数量（头） 被注射的数量（头） 250 260 270 280 290 300 180 171 161 152 142 133

问题四：

用人工手段可以使得该动物群落的数量减少，使其数量控制在一定范围内，但是，我们并不知道用了这种人工手段后，对于一次突发事件后要增加其数量时，该种群的复壮能力如何。因此，我们建立数学模型，分析研究其复壮能力。我们假设由于某种原因（例如由于一场大灾难导致该动物的死亡），使得该动物减少至M头，此时不得不停止注射避孕的方法。在初步模型中，我们先假设灾难发生后该种群虽然数量减少，但其年龄结构没有改变。建立该动物群落重新壮大的能力的模型如下： 模型七：

?d?(k)?t?60(k)60(k)??xi,1??xi,0－－－－－－－(1)?i?1i?1?dt?60(1)??xi,0?M－－－－－－－－－－－－－(2)??i?1k?1)k)?xi(,0?xi(,1.....(i?0,1,..59)－－－－－－－(3)?k)?xi(,1?0...........(i?60)－－－－－－－－－(4) ?(k)(k)?xi,1??1xi,0......(i?0)－－－－－－－－－(5)?(k)(k)x??x?2i,0.....(i?1,2,...,59)－－－－－－(6)?i,1对M取不同的值，利用c语言进行编程求解，（源程序见附录9）可以得到以下表： M值 第t年 t=1 增量d1 665 594 526 456 388 322 250 179 113 10000 9000 8000 7000 6000 5000 4000 3000 2000 当前数量M1 t=2 增量d2 当前数量M2 t=3 10665 9594 707 638 8526 562 7456 489 7945 521 8466 557 6388 414 6082 444 7246 472 7718 507 8225 538 5322 342 5664 365 6029 390 6419 418 6837 445 7272 472 7744 507 8251 543 4250 267 4517 286 4803 305 5108 325 5433 348 5781 371 6152 401 6553 426 6979 456 7435 22.8061 3179 191 3370 205 3575 219 3794 235 4029 251 4280 268 4548 286 4548 307 4835 327 5142 16.4364 2113 122 2235 130 2365 140 2505 145 2650 156 2806 166 2972 178 3150 191 3341 204 3545 10.7636 11372 10232 9088 757 恢复 699 602 增量d3 当前数量M3 10932 9690 729 恢复 640 t=4 增量d4 当前数量M4 10330 9023 689 596 t=5 增量d5 当前数量M5 11019 9619 741 恢复 638 t=6 增量d6 当前数量M6 10257 8792 682 581 t=7 增量d7 当前数量M7 10929 9373 729 恢复 623 9996 665 t=8 增量d8 当前数量M8 t=9 增量d9 当前数量M9 10661 8794 707 581 t=10 增量d10 当前数量M10 11368 9375 46.6000 42.6857 38.8810 35.6424 28.6121 d?2?t?46.00增长加速度：D? 2dt00

该表对M进行不同的取值，分别表示某一年（t＝0）灾难发生后该动物群落剩下的数量。我们用第一题的C程序即可算出该群落在下一年（t＝1）的增量，把这个增量加上上一年群落的总数即为第一年（t＝1）群落的总数M1,接着我们又用程序算出这一年的增量，加上M1即为第二年（t＝2）的个体总数M2。以这样的算法推算下去，一直到群落数量恢复到11000左右或者前十年（t＝10）为止。

将计算出来的各年增量按时间t（1，2??10）列在座标上，用MATLAB进行曲线拟合（具体程序见附录10），得到的图形如下所示：

从图中可以看出，每年的增量都是遵从一定的线性关系的，也就是说，该群落的增

d?2?t?长加速度是一定的，通过MATLAB的曲线拟合可以求出不同M值所对应的增2dt长加速度(程序见附录 )，数据如下： M值 10000 9000 46.6000 8000 42.6857 7000 38.8810 6000 35.6424 5000 28.6121 4000 22.8061 3000 16.4364 2000 10.7636 D（只/46.00年^2） 00

通过这些数据，我们可以知道，在灾难后群落数量减少至M时的重新壮大能力。

再用MATLAB求出M与D的关系，进行曲线拟合求出斜率k后画出图形（程序见附录11），如下所示

可知M与D是服从一定的线性关系的，由MATLAB算出的的数据可知

M＝3.9010+0.0047D。

该模型的稳定性分析见本论文第四部分。

四． 模型的稳定性分析

1．对于问题一：我们利用C语言进行编程，求得了第k年时各年龄结构的分布情况，并利用MATLAB软件对得到的数据进行画图分析，从输出的图可看出各年该动物各个年龄的分布是围绕一定的轨迹进行分布的，从而说明了该模型是稳定的，其次，用统计软件SAS对该模型得出的数据进行相关系数的分析，得到相关系数R都接近于1，而P＝0.0001<0.01则表示其相关性显著。故建立的模型一是很稳定的，其各年的该动物的总数基本上都在11000左右。

2．对于问题二、问题三：我们从两个不同方面对问题进行分析，分别从?＝新生的幼仔数量－自然死亡的数量～0，以及?＝该动物群落在年底时的总数量与年初的数量加上被运走的数量b的差值～0两个不同方面进行模型的建立，从而把该动物群落的稳定性问题转化为求目标函数最小化问题。（这是本论文的一个亮点）。利用MATLAB对模型进行求解，得出的数据最后都能使得该动物群落维持稳定，可见对于问题二、问题三所建立的模型是很稳定的，结果也是很优的。

3．对于问题四，我们用C语言的编程求得了不同M值在10年内各年的增量，求出该群落在不同M值所对应的增长加速度D。为了求出不同的M与对应的D的关系，我们在MATLAB中输入M与D的数组，进行稳定性分析（程序见附录11）。在回归分析及检验中，我们得到，?0＝3.9010，?1＝0.0047；?0的置信区间为[-1.5746 , 9.3766], ?1的置信区间为[0.0039, 0.0055]; r2=0.9616,F=175.1021,p=0.0000,p<0.05,可知回归模型 M＝3.9010+0.0047D 成立 对模型进行残差分析,作残差图,得图如下:

从残差图可以看出,除最后一个数据外,其余数据的残差离零点均较近,且残差的置信区间均包括零点,这说明回归模型M＝3.9010+0.0047D能较好的符合原始数据,而最后一个数据视为异常点.(造成异常点的原因是M=10000时群落只要两年就可恢复到原来的数量11000,数据太少无法作出较接近实际的曲线拟合。 在预测及作图后得各数据点及回归方程得图形如图：

可以看出,只有最后一个数据点离回归直线距离较远是异常点。 由以上分析可知M与D的关系是比较稳定的。 五． 模型优缺点及改进方向 我们的模型有以下优点：

1．模型的稳定性好，模型给出的方案使得该动物群落的各年龄段的数量保持平衡，从而满足了维持该动物群落的稳定。

2．模型的适用范围广，易于推广，模型对于其它生态经济现象（如最优捕鱼策略问题）同样适用。

3．基本模型对问题的描述准确、合理、推导严谨，理论性强； 4．模型结合实际，具有很高的实用价值。 我们的模型有以下缺点：

1. 由于题目给的数据不够多，所以使得该模型无法更加接近实际的情况。

2. 由于题目没有以前各年偷猎枪杀该动物的数据，故模型没有考虑偷猎枪杀该动物所产生的影响。

3. 对于一个动物群落，其影响的不确定性因素很多，而模型只是局限在题目中提到的，故与实际情况有一定的偏差。

模型的改进方向：

1．由生活常识可知，对于该动物群落超过60岁后的死亡情况应基本上服从正态分布。故可在模型一中的死亡数中加上这个因素。

2．由于动物的产幼仔是随机分布的，故可利用模拟仿真的方法建立更具有真实性的模型。

参考文献

【1】 樊欣,邵谦谦，SAS 8.X经济统计，北京希望电子出版社，2003.3. 【2】 苏金明，张莲花，MATLAB工具箱应用，北京电子工业出版社，2004.1. 【3】 姜启源，谢金星，数学模型（第三版），高等教育出版社，2003.8. 【4】 余世孝，数学生态学导论，科学科技文献出版社，1995.8. 附录一:（模型一用C语言实现的过程与运行结果） main()

{int x[70],n; /* 定义一个长度为70的数组,每个单位的值表示各年龄的数

量,定义变量n */

int sum,sum1,sum2; /* 三个变量分别代表该年群落个体数量,13—59岁

雌性数量,1—59岁个体数量 */

clrscr(); sum=0; sum1=0; sum2=0;

printf(\

scanf(\/* 输入当年出生的幼仔数量,设幼仔年龄为0岁 */ printf(\

x[1]=0.75*x[0]; /* 幼仔从0岁向1岁过渡时存活率为75% */ for(n=2;n<=60;n++)

{x[n]=0.95*x[n-1]; /* 向下一年龄段过渡时的存活率为95% */ }

for(n=0;n<=60;n++)

{sum=sum+x[n]; /* 对各年龄段数量累加得到种群总体数量 */ printf(\ }

for(n=13;n<=60;n++)

{sum1=sum1+x[n]; /* 对13—60岁个体数量累加 */ }

for(n=1;n<=59;n++)

{sum2=sum2+x[n]; /* 对1--59岁个体数量累加 */ }

printf(\

printf(\/* 13—60岁个体总数的一 printf(\半即为该年龄段雌性数量

再除以3.5即为当年生幼仔的雌性的数量 */

通过不断输入新生幼仔的初值,直到算出当年群落的个体总数与题设数据吻合,则该组数据即为当年群落年龄结构的分布(未算入当年运走的数量). 运行结果:

1.前兩年运走的个体数量为876头,因此前一年的个体总数为11000+876=11876(头),通过调试程序使总数sum达到11876左右. 第一次调试结果:

2．前两年年末即前一年年头所剩下的13—60岁的雌性个体数量决定了前一年年末个体的数量即2470/3.5＝706头；

前两年剩下的数量加上前一年新生幼仔的数量 11000+706=11706，其值为前一年年末的个体数量

用程序算得11714，误差为8头，其误差由于变量是整型变量，数字取整数所致。 第二次调试结果:

在前一年的年末运走622头，即该群落的个体数量剩下11714-622=11092 根据表得到2503/3.5=715，所以11092+715=11807为该群落今年年尾的个体总数。

程序计算得11808，误差为1，其误差由于变量是整型变量，数字取整数所致。

3.对于以上数据是否符合假设,我们还要用一组初步数据來进行对比.因此我们设每年运出的个体数量为0, 个体总数为11000 头,通过调试程序使总数sum达到11000左右.

第三次调试结果:

该结果为当每年没有个体被运走时数量停留在11000左右的数据，用于与初步数据进行对比。

附录2（用MATLAB画出模型一轨迹）： clear all;clc

y1=[792 594 564 535 508 482 457 434 412 391 371 352 334 317 301 285 270 256 243 230 215 207 196 186 176 167 158 150 142 134 127 120 113 107 101 95 90 85 80 75 71 67 63 59 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

y2=[806 604 573 544 516 490 465 441 418 397 377 358 340 322 305 289 274 260 246 233 221 209 198 188 178 169 160 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

y3=[746 559 531 504 478 454 431 409 388 368 349 331 314 298 283 268 254 241 228 216 205 194 184 174 165 156 148 140 132 125 118 112 106 100 94 89 84 79 75 71 67 63 59 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 17]';

y4=[800 600 569 540 512 486 461 437 415 394 374 355 337 320 303 287 272 258 245 232 220 208 197 187 177 168 159 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

x1=[0 0 0 0 3 4 7 20 9 15 9 22 3 23 5 13 21 0 22 14 5 13 10 0 13 30 14 12 0 20 6 3 5 8 12 10 3 7 14 10 16 21 13 10 12 6 3 6 9 13 10 3 6 21 15 4 13 10 32 14 0]';

x2=[0 0 20 21 13 12 13 22 14 40 14 26 13 14 27 3 14 12 20 25 17 14 10 0 2 3 4 4 3 2 3 13 16 13 10 10 12 16 12 10 12 19 13 24 17 16 25 12 45 23 34 13 16 10 17 13 13 12 3 22 20]';

m1=y1-x1; m2=y2-x2; m3=y3; m4=y4; t=0:60

plot(t,m1,'b.'),hold on,plot(t,m4,'m-')

legend('前一年剩余数量的年龄结构情况','前两年剩余数量的年龄结构情况','没有被运走时的年龄结构情况','预测当前的年龄结构情况（无运走）') title('今年时该动物群落的年龄结构情况')

附录3（用SAS对模型一进行检验的程序及调试结果）：

data A;

input x1 x2 x3 x4 @@; cards;

594 604 559 600 564 553 531 569 535 523 504 540 505 503 478 512 478 478 454 486 450 452 431 461 414 419 409 437 403 404 388 415 376 357 368 394 362 363 349 374 330 332 331 355 331 327 314 337 294 308 298 320 296 278 283 303 272 286 268 287 249 260 254 272 256 248 241 258 221 226 228 245 216 208 216 232 213 204 205 220 194 195 194 208 186 188 184 197 186 188 174 187 163 176 165 177

on,plot(t,m2,'y*'),hold on,plot(t,m3,'g--'),hold

137 166 156 168 144 156 148 159 138 147 140 151 142 140 132 143 114 133 125 135 121 125 118 128 117 108 112 121 108 98 106 114 99 95 100 108 89 92 94 102 85 86 89 96 87 79 84 91 78 70 79 86 66 69 75 81 65 66 71 76 55 60 67 72 46 49 63 68 50 51 59 64 49 36 56 60 44 39 53 56 47 37 50 53 47 25 47 50 41 35 44 47 35 -1 41 44 28 18 38 41 28 4 36 38 33 23 34 36 28 18 32 34 11 22 30 32 15 13 28 30 24 15 26 28 13 13 24 26 14 12 22 24 -10 19 20 22 6 -2 18 20 18 -2 17 18 ;

proc corr; var x1 x2 x3 x4; run;输出结果如下：

附录4（用MATLAB求解目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2857 0.1818 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 529 ]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

>> Optimization terminated successfully. x =

1.0e+003 * 2.2934 0.2756 0.7053 0.0000

fval = ％ 返回目标函数值 ％ 2.0286e-016

exitflag = ％ 表示求得的解是收敛的 ％ 1 output =

iterations: 6 cgiterations: 0

algorithm: 'lipsol' ％ 求目标最优化的方法 ％ lambda = ％ 所占的空间 ％ ineqlin: [0x1 double] eqlin: [3x1 double] upper: [4x1 double] lower: [4x1 double]

附录5（用MATLAB求解模型三目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2896 0.1843 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 529 ]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

> Optimization terminated successfully. x =

1.0e+003 * 2.2020 0.3670 0.7053 0.0000 fval = 2.4791e-016

exitflag = 1 output =

iterations: 6 cgiterations: 0 algorithm: 'lipsol'

附录6: （用MATLAB求解模型四目标最优化问题） clear all;clc f=[0;0;0;0;1]; A=[0 0 1 0 0]; b=[242];

Aeq=[0.2896 0.1843 0.1698 -1 0;1 1 1 0 0;0 0 0 0.75 -1]; beq=[0 2569 529]; lb=zeros(5,1);

[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) 输出结果如下：

>> Optimization terminated successfully. x =

1.0e+003 * 2.2171

0.2421 ％由此可知当?＝242时模型稳定％ 0.1098 0.7053 0.0000 fval = 1.6409e-012 exitflag = 1 output =

iterations: 6

cgiterations: 0 algorithm: 'lipsol' lambda =

ineqlin: 1.9145e-014 eqlin: [3x1 double] upper: [5x1 double] lower: [5x1 double]

附录7（用MATLAB求解模型五目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2896 0.1843 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 579]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

>> Exiting: One or more of the residuals, duality gap, or total relative error has grown 100000 times greater than its minimum value so far: the primal appears to be infeasible (and the dual unbounded). (The dual residual < TolFun=1.00e-008.) x =

1.0e+003 * 2.5685 0.0005 0.7439 0.0000 fval = 4.8164e-006

exitflag = ％表明得到的数据是发散的％ -1 output =

iterations: 3 cgiterations: 0 algorithm: 'lipsol' lambda =

ineqlin: [0x1 double] eqlin: [3x1 double] upper: [4x1 double] lower: [4x1 double]

附录8（用MATLAB求模型六目标最优化问题）

clear all;clc for i=50:10:300 f=[]; A=[]; b=[]; Aeq=[]; beq=[]; lb=zeros(1,1);

[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) end

输出结果如下： x ＝

370 361 351 342 332 323 313 304 294 285 275 266 256 247 228 218 209 199 190 180 171 161 152 142 133

附录9（用C语言求不同M值时各年的增量） main() {int x[70],n,d;

int sum,sum1,sum2; clrscr(); sum=0; sum1=0; sum2=0;

printf(\ scanf(\

printf(\ x[1]=0.75*x[0]; for(n=2;n<=60;n++) {x[n]=0.95*x[n-1]; }

for(n=0;n<=60;n++) {sum=sum+x[n];

printf(\ \ /*打印出各年龄段的数量*/ }

for(n=13;n<=60;n++) {sum1=sum1+x[n]; }

printf(\ /*用总和来控制，使之等于M在第t年数量*/ printf(\/*每年的增量约等于13—16岁雌性个体的数量

*/

}

附录10（求不同M值增量d的关系） clear all;clc

t=[1 2 3 4 5 6 7 8 9 10]; t1=[1 2]; t2=[1 2 3]; t5=[1 2 3 4]

t3=[1 2 3 4 5 6]; t4=[1 2 3 4 5 6 7 8]; M1=[665 707 757]; M2=[594 638 699 729]; M3=[526 562 602 640 689 741];

M4=[456 489 521 557 596 638 682 729]; M5=[388 414 444 472 507 538 581 623 665 707]; M6=[322 342 365 390 418 445 472 507 543 581]; M7=[250 267 286 305 325 348 371 401 426 456]; M8=[179 191 205 219 235 251 268 286 307 327]; M9=[113 124 134 141 152 161 173 185 199 212]; f=inline('k(1)*t+k(2)','k','t'); k1=[0.5 30];

m1=lsqcurvefit(f,k1,t2,M1) z1=f(m1,t2) k2=[0.5 30];

m2=lsqcurvefit(f,k2,t5,M2) z2=f(m2,t5) k3=[0.5 30];

m3=lsqcurvefit(f,k3,t3,M3) z3=f(m3,t3) k4=[0.5 30];

m4=lsqcurvefit(f,k4,t4,M4) z4=f(m4,t4) k5=[0.5 30];

m5=lsqcurvefit(f,k5,t,M5) z5=f(m5,t) k6=[0.5 30];

m6=lsqcurvefit(f,k6,t,M6) z6=f(m6,t)

k7=[0.5 30];

m7=lsqcurvefit(f,k7,t,M7) z7=f(m7,t) k8=[0.5 30];

m8=lsqcurvefit(f,k8,t,M8) z8=f(m8,t) k9=[0.5 30];

m9=lsqcurvefit(f,k9,t,M9) z9=f(m9,t)

plot(t2,M1,'b.'),hold

on,plot(t2,z1,'r-'),plot(t5,M2,'b.'),hold

on,plot(t5,z2,'r--'),plot(t3,M3,'b.'),hold on,plot(t3,z3,'y-'),plot(t4,M4,'b.'),hold on,plot(t4,z4,'y--'),plot(t,M5,'b.'),hold on,plot(t,M6,'b.'),hold

on,plot(t,z5,'g-'),hold

on,plot(t,z6,'g--'),plot(t,M7,'b.'),hold

on,plot(t,z8,'m--'),hold

on,plot(t,z7,'m-'),plot(t,M8,'b.'),hold on,plot(t,M9,'b.'),hold on,plot(t,z9,'k-') title('M取不同的值时每年增量的分布情况')

其运行结果如下：

Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m1 =

46.0000 617.6667 z1 =

663.6667 709.6667 755.6667 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m2 =

46.6000 548.5000 z2 =

595.1000 641.7000 688.3000 734.9000 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m3 =

42.6857 477.2667 z3 =

519.9524 562.6381 605.3238 648.0095 690.6952 733.3810 Optimization terminated successfully:

Norm of the current step is less than OPTIONS.TolX m4 =

38.8810 408.5357 z4 =

Columns 1 through 7

447.4167 486.2976 525.1786 564.0595 602.9405 641.8214 680.7024 Column 8 719.5833

Optimization terminated successfully:

Norm of the current step is less than OPTIONS.TolX m5 =

35.6424 337.8667 z5 =

Columns 1 through 7

373.5091 409.1515 444.7939 480.4364 516.0788 551.7212 587.3636 Columns 8 through 10 623.0061 658.6485 694.2909 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m6 =

28.6121 281.1333 z6 =

Columns 1 through 7

309.7455 338.3576 366.9697 395.5818 424.1939 452.8061 481.4182 Columns 8 through 10 510.0303 538.6424 567.2545 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m7 =

22.8061 218.0667 z7 =

Columns 1 through 7

240.8727 263.6788 286.4848 309.2909 332.0970 354.9030 377.7091 Columns 8 through 10 400.5152 423.3212 446.1273 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m8 =

16.4364 156.4000 z8 =

Columns 1 through 7

172.8364 189.2727 205.7091 222.1455 238.5818 255.0182 271.4545 Columns 8 through 10 287.8909 304.3273 320.7636 Optimization terminated successfully:

First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected m9 =

10.7636 100.2000 z9 =

Columns 1 through 7

110.9636 121.7273 132.4909 143.2545 154.0182 164.7818 175.5455 Columns 8 through 10 186.3091 197.0727 207.8364

附录11（M、D关系的稳定性检验）

x=[2000 3000 4000 5000 6000 7000 8000 9000 10000];

y=[10.7636 16.4364 22.8061 28.6121 35.6424 37.5357 40.4000 46.6000 46.0000]; f=inline('k(1)*x+k(2)','k','x'); k=[0.5 30];

m=lsqcurvefit(f,k,x,y) z=f(m,x)

plot(x,y,'.'),hold on,plot(x,z,'-') 其输出结果如下：

>> Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m =

0.0046 4.0893 z =

Columns 1 through 7

13.2742 17.8666 22.4591 27.0516 31.6440 36.2365 40.8290 Columns 8 through 9 45.4214 50.0139

附录12（问题四稳定性分析）

>> x=[2000 3000 4000 5000 6000 7000 8000 9000 10000]'; >> X=[ones(9,1)x];

Warning: Future versions of MATLAB will require that whitespace, a comma, or a semicolon separate elements of a matrix. Please type

\>> Y=[10.7636 16.4364 22.8061 28.6121 35.6424 38.8810 42.6857 46.6000 46.0000]';

>> [b,bint,r,rint,stats]=regress(Y,X); >> b,bint,stats b = 3.9010 0.0047 bint =

-1.5746 9.3766 0.0039 0.0055 stats =

0.9616 175.1021 0.0000 >> rcoplot(r,rint) >> z=b(1)+b(2)*x z = 13.2832 17.9743 22.6653 27.3564 32.0475 36.7386 41.4296 46.1207 50.8118

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

本文来源：https://www.bwwdw.com/article/dfw.html