大象群落的稳定发展数模之欧阳文创编

封一
答卷编号（竞赛组委会填写）：
答卷编号（竞赛组委会填写）：
论文题目大象群落的稳定发展
参赛队员：
1. 黄立敬电话：
2. 陈光电话：
3. 陈灵电话：
大象群落的稳定发展
摘要
本文根据非洲某国的国家公园近两年内从公园运出的大象的大致年龄和性别的统计情况，探讨大象的合理的存活率并推测当前的年龄结构，针对不同情况给出如何进行避孕注射以达到控制大象数量的目的。

首先，充分利用给出的近两年来运出的大象的数量与性别统计表，分析近两年来的大象群落的情况，建立一个线性方程组的数学模型，通过求解方程组得到年龄在2岁到60岁之间的大象的总数，并且求出了存活率为：98.9718%；因为假设公园内2岁到60岁之间的大象占总大象的比例等于运出的2岁到60岁之间的大象占总移出大象的比例，所以通过一些比例之间的关系得到这个大象群落的当前的年龄结构（见表1）。

然后，建立一个按年龄分组的种群增长的差分方程模型，运用第一问求出的各年龄段大象的存活率以及繁殖率，
求解当前大象群落对应的Leslie矩阵的特征根，发现该特征根大于1，根据Leslie矩阵的稳定性理论知道：如果不进行避孕注射该大象种群将无限增长（如果环境允许）；据此，利用Leslie矩阵稳定的充要条件求出应该保持多大的繁殖率才能使种群保持稳定，求解的主要思路是：特征根取为1、把繁殖率当成未知数，将此时的各年龄段的存活率代入方程VI，求解这个以繁殖率为未知数的方程可以得到要使种群保持稳定繁殖率的取值；根据需要避孕掉母象所生的幼象的数目等于注射避孕药使得母象没有繁殖幼象的数目这一条件建立一个方程，最后求得每年注射避孕药的母象头数为：1393（头）。

最后，假设被转移的大象只考虑处于1—60岁之间，这样可以认为转移后的大象看成每年多死了这么多头大象，即意味着死亡率将增加，存活率将减少；仍然按照解决第二问的模型，只需将此时不同的各年龄段大象的存活率代入那个以繁殖率为未知数的方程（方程VI），求出应该保持多大的繁殖率才能使种群保持稳定。

考虑到求解的数据比较多，采取计算机模拟的方法来确定移出大象后所需要进行避孕的母象头数（见表2），为了检验计算机模拟的正确性，用理论去验证。

模拟的思路方法见计算机模拟流程图—图2。

关键字：线性方程组、差分方程模型、Leslie矩阵、计算机模拟
问题重述
位于非洲某国的国家公园中栖息着近11000头大象。

管理者要求有一个健康稳定的环境以便维持这个11000头大象的稳定群落。

管理者逐年统计了大象的数量，发现在过去的20年中，整个大象群经过一些偷猎枪杀以及转移到外地还能保持在11000头的数量，而其中每年大约有近600头到800头是被转移的。

由于近年来，偷猎被禁止，而且每年要转移这些大象也比较困难，现决定采取避孕注射法以维持大象数量的平衡。

我们已知此公园近两年内从这个地区运出的大象的大致年龄和性别的统计。

根据这些信息我们需要解决以下问题：
1.探讨年龄在2岁到60岁之间的象的合理的存活率的模
型，推测这个大象群落的当前的年龄结构。

2.估计每年有多少母象要注射避孕药，可以使象群固定在
11000头左右。

这里不免有些不确定性，是否能估计这种不确定性的影响。

3.假如每年转移50至300头象到别处，那么上面的避孕措
施将可以有怎样的改变？
问题假设
1、假设大象的性别比近似认为1：1，并且采用措施维持这个性别比；
2、假设母象可以怀孕的年龄为11岁—60岁、最高年龄为70岁，70岁的死亡率为100%，并且61—70岁的大熊的头数
呈线性递减；
3、假设大象在各年龄段中的分布率不变，即年龄结构不变，并采用各种措施维持这一结构；
4、假设被转移的大象只考虑处于1—60岁之间，转移后的大象看成每年多死了这么多头大象；
5、假设0岁大象能够活到1岁的比例为75%；
符号说明
X：表示一年中大象的头数（i=0表示0岁大象的头数，i
i=1表示1--60岁大象头数，i=2表示61—70岁大象的头数）；
p：表示存活率（0p表示0岁大象的存活率，1p表示1—i
60岁大象的存活率，
p表示61岁—70岁大象的存活
2
率）；
x k：表示时段k第i年龄组的大象数量；
()
i
b：第i年龄组每个（母象）个体在1个时段内平均繁殖i
的数量；
s：第i年龄组的存活率；
i
L： Leslie矩阵；
：L矩阵的那个唯一正特征根；
1
n：表示移出大象的头数；
问题分析
对于问题一，利用给出的近两年来运出的大象的数量与性别统计表，可以分析近两年来的大象群落的情况，比如移
出的各个年龄段的大象占移出的总的大象的头数的比例是多少，还可以根据两年移出大象后大象总数都是11000来建立方程，用于求解存活率。

对于问题二，因为考虑的是公园在未来很长一段时间的大象种群控制问题，所以可以建立一个按年龄分组的种群增长的差分方程模型，根据差分方程的 Leslie矩阵的特征根，结合 Leslie矩阵的稳定性理论对当前大象种群的情况进行分析。

为了保持大象种群的稳定，必须使得Leslie矩阵的最大特征根为1，而这样，特征根取为1、把繁殖率当成未知数，将此时的各年龄段的存活率代入方程特征方程，求解这个以繁殖率为未知数的方程可以得到要使种群保持稳定繁殖率的取值；根据需要避孕掉母象所生的幼象的数目等于注射避孕药使得母象没有繁殖幼象的数目这一条件建立方程来求解应该对多少头母象进行避孕。

对于问题三，由于假设被转移的大象只考虑处于1—60岁之间，故可以认为转移后的大象看成每年多死了这么多头大象，即意味着死亡率将增加，存活率将减少。

按照解决第二问的模型，只需将此时不同的各年龄段大象的存活率代替原来的存活率，就可以求出此时应该保持多大的繁殖率才能使种群保持稳定。

为了方便，可用采用计算机模拟的方法来确定移出的大象在哪个年龄段，考虑到计算机模拟的不确定性，可以对模拟结果进行检验。

探讨大象的存活率和当前大象的年龄结构
下面将根据给出的近两年来运出的大象的数量与性别统
计表，分析近两年来的大象群落的情况，建立一个线性方程
组数学模型，通过求解方程组得到年龄在2岁到60岁之间
的大象的存活率，并给出大象各年龄所占的比例，进而得到
这个大象群落的当前的年龄结构。

1、线性方程组模型的建立
（1）首先，计算一年中大象的头数。

大象群是由0岁，1—60岁，61岁—70岁组成 ,且稳定
在11000头。

设0岁的头数为X0，1—60岁大象头数为X1，
61岁—70岁大象头数为X2。

所以得到第一个方程：
X0+X1+X2=11000 （I）
（2）其次，考虑到前一年大象的总数等于前两年存活
下来的大象加上新生的幼儿再减去运出的大象数。

设0岁大象的存活率为
p，1—60岁大象的存活率为
p，61岁—70岁大象的存活率为2p。

则经过一年后，新生的
1
大象存活下来的头数为X0⨯
p；1到60岁的大象存活下来的
头数为X1
p⨯；61岁——70岁的大象能存活下来的头数为
1
X2
p⨯，因此得到第二个方程：
2
（X0⨯
p+ X11p⨯+ X22p⨯）+ X0-622=11000
（II）
联立（I）、（II）得到方程组：
0120011220X +X +X =11000X + X + X + X -622=11000 p p p ⎧⎨⨯⨯⨯⎩
（*）
2、模型的求解
根据近两年来运出的大象的数量与性别统计表，得到如下分析结果：
（1）计算0岁的大象头数
由表中统计，1岁—10岁的大象占1岁—60岁的大象比
例为：（67/620+169/876）/2=15.05%
所以得到：11岁—60岁能生小象的母象占1岁—60岁的大象比例为：
（1-15.032%）⨯0.5 =42.48%
因为能生小象的母象每3.5年生一头小象，且双胞胎的机会为1.35%，相当于每年生0.2896头 ，所以0岁的大象占1岁—60岁的大象比例为：
0.4248⨯0.2896=0.12303
这样0岁的大象共有：
0X =0.12303⨯1 X
（III ）
（2）计算60岁—70岁的大象头数
从表中计算运出的59岁的大象占运出的总大象比率为：
（14/622+22/876）/2=0.0238
由于运出的大象都是1岁—60岁的，所以0.0238也可
看为59岁的大象占1—60岁的大象的头数比例，得到60岁的大象占的比例为0.02381p ⨯，由假设可以知道：
61岁—70岁的大象头数为：
2X =1/2⨯10⨯0.02381p ⨯⨯X 1
（ IV ）
60岁——70岁的大象经过一年能存活下来的头数为：
2 211X =(1/2)90.0238X p p ⨯⨯⨯⨯⨯
（V ）
（3）、将（III ）、（V ）和（IV ）两个式子代入上面方程组（*）得：
又由假设知道，0岁大象的存活率为0p =75%代入上述方程
组，然后用Mathematica 解之得：
再依次将1X 、1s 代入（III ）、（V ）和（IV ）求得：
所以， 0岁大象的总头数为1091（头）；1—60岁的大象的存活率为98.9719%，总头数为8865（头）；61岁—70岁的大象头数为1091（头）。

把0—70岁的大象分为八个年龄段，由假设知道，各个年龄段占总数可以用各个年龄段移出的头数除以移出的总头数来衡量。

下面以1—10年龄段的大象头数计算为例：
前一年总共移出622头，其中1—10岁移出为67头；前两年总共移出876头，其中1—10岁移出169头。

故1—10年龄段的大象头数可以这样计算：
11X =671698865 [()/2]622 876
⨯+=1332（头） 其他的年龄段用同样的方法计算，得到如下表（附饼形图）：
表1（大象年龄结构）
图1（大象年龄结构饼图）
3、结果分析
（1）由结果可以知道，2—60岁大象的存活率为98.9718%，这与题目给出的大于95%是相一致的，所以可以认为结果是合理的；
（2）从图1可以看出，各个年龄段的大象所占的比例基本上是一样的，21—30岁和41—50岁的大象比例相对比较大，因为这段大象正处于年龄的黄金时期。

由此，可以认为求出的大象年龄也是合理的。

估计每年注射避孕药的母象头数
为了估计每年注射避孕药的母象头数，首先建立一个按年龄分组的种群增长的差分方程模型；然后用Leslie 矩阵稳定的充要条件分析如果不进行避孕注射种群的增长情况；最后仍然利用Leslie 矩阵稳定的充要条件求出应该保持多大的繁殖率才能使种群保持稳定，进而利用一个方程求出每年注射避孕药的母象头数。

1、按年龄分组的种群增长的差分方程模型的建立
记时段k 第i 年龄组的大象数量为()i x k ，k=0，1，2，…，i=1，2，…n,第i 年龄组的繁殖率为i b ，即第i 年龄
组每个（母象）个体在1个时段内平均繁殖的数量，第i 年龄的存活率为i s ，我们这里假设i b 和i s 不随时段k 变化，在稳
定的环境下这个假设是合理的。

i b 和i s 可由统计资料获得。

()i x k 的变化规律由以下的基本事实得到：时段k+1第1年龄组种群数量是时段k 各年龄组繁殖数量之和，即
时段k+1第i+1年龄组的种群数量是时段k 第i 年龄组存活下来的数量，即
记时段k 种群按年龄组的分布向量为
由繁殖率i b 和存活率i s 构成的矩阵为
根据Leslie 矩阵的性质可以得到如下的定理：
定理1：L 矩阵有唯一的正特征根1λ，且它是单重根的，
1λ对应正特
征向量 *12111
22111
1...1,,,...,T n n s s s s s s x λλλ--⎡⎤=⎢⎥⎣⎦ L 矩阵的其他n-1个特征根k λ都是满足 1,2,3,...,k k n λλ≤= 该定理表明L 矩阵的特征方程
只有一个正根，并且易知，**1Lx x λ=
2、如果不进行避孕注射种群的增长情况
（1）建立Leslie 矩阵
首先，由第一问的求解知道，0岁的大象的存活率为0.75；1—60岁大象的存活率为0.989718；根据假设61—70
岁大象头数是线性递减的，而且到70岁所有的大象都死完
了，所以很容易求出存活率为（1-0.1）
1s ⨯=0.9⨯0.989718=0.8907；11—60岁大象的繁殖率为
0.1448。

然后根据上面的矩阵L 建立起如果不进行避孕注射种群的
增长的Leslie 矩阵如下所示：
7171L ⨯=
00...00.14480.1448......0.14480...000.750...000......00...0000.9897...000......00 (000)
0...000......00...0000...0.989700......00 (000)
0...00.98970......00...00.........000.9897......00...00...
......000......00...0000...000......00...0000...000......0.98970 (000)
0...000......00.9897...0000...000......00...0000...000......00...0.98970⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦这
是一个71⨯71的矩阵。

（2）讨论7171L ⨯的特征根，分析种群增长规律
用Matlab 软件求得特征根为R=1.0481，根据定理1知
道，如果不进行避孕注射，该大象种群将无限增长下去（如
果环境允许），所以要进行避孕注射。

3、求出每年注射避孕药的母象头数
根据Leslie 矩阵的性质知道，要保持种群稳定，必须
使得特征根r=1，即使得下面式子成立：
010201011220121.........1n n n n b b s b s s s b s s s b s s s s ---+++++=
（VI ）
具体到本题来就是使得如下成立：
解这个方程求出要保持大象种群的稳定，繁殖率应该为
0b =0.0377
保持大象种群数量不变的繁殖率b 0与没采取避孕时的繁
殖率b 有一定的差距，所以需要避孕掉具有（b- b 0）繁殖率
母象所生的幼象。

假设每年要避孕0n 头大象，由于一次注射
可以使得一头成熟的母象在两年内不会受孕，所以每年实际上共有2
n 头大象处于避孕期。

这样根据需要避孕掉具有（b- b0）繁殖率母象所生的幼象的数目等于注射避孕药使得母象没有繁殖幼象的数目这个条件得到一个方程：
解之得
n=1393
所以每年注射避孕药的母象头数为：1393（头）
4、分析不确定因素的影响
（1）最初一两年避孕母象发情期增多，与未避孕母象产生竞争求偶的公象，使部分能怀孕的母象不能怀孕。

而避孕的母象每月发情一次，会扰乱了正常求偶的母象，这样会造成未避孕母象的繁殖率出现下降，避孕的母象数量应该减少。

（2）随着时间的增长，如果持续使用避孕药，会使象的年龄结构发生变化，象的结构呈老龄化，所以随着时间的增长，要保证象群的稳定，避孕药的使用量必定会逐年减少直至禁用。

考虑转移大象时母象的避孕策略
被转移的大象只考虑处于1—60岁之间，转移后的大象看成每年多死了这么多头大象，即意味着死亡率将增加，存活率将减少。

下面首先通过计算机模拟来确定移出大象后所需要进行避孕的母象头数；然后用理论去验证计算机模拟的正确性。

1、通过计算机模拟确定需要进行避孕的母象头数
产生了n（可自己指定）个0~1的随机数。

具体算法如
下页图所示。

图2（计算机模拟流程图）
下面议n=100为例进行计算机模拟。

令n=100,进行10次计算机模拟，得到当运出大象头数为100头时，要使大象头数稳定在11000头需避孕的大象头
数如下：
第一次模拟：1164（头）
第二次模拟：1195（头）
第三次模拟：1192（头）
第四次模拟：1167（头）
第五次模拟：1206（头）
第六次模拟：1163（头）
第七次模拟：1207（头）
第八次模拟：1190（头）
第九次模拟：1173（头）
第十次模拟：1178（头）
10次模拟得到需避孕的大象头数的平均数为1184头。

因此可以认为当运出大象头数为100头时，要使大象头数稳定在11000头需避孕的大象头数为1184头。

同理，可以得到题目中要求的当运出大象50—300头要使大象头数稳定在11000头需避孕的大象头数分别为：
表2（移出大象与对应得避孕母象头数）
图3（运出大象数目与应该避孕母象数目的关系）
从该图的变化趋势可以看出为了使大象总数稳定在
11000头而应该避孕的母象数目与运出大象的头数近似成线
性关系。

通过最小二乘进行拟合二次函数得到移出大象头数
x 与避孕母象头数y 关系：20.00344 1.75581396.2y x x =--+
画出图像为：
图4（原始数据关系与拟合关系对比
图）
2、用理论去验证计算机模拟的正确性
假设所转移的大象是有目的的挑选的，即挑选大象时是
按照第i （i=1，2，3，4，5，6）年龄组所占比例进行的，
这也是符合情理。

设每年被转移的大象共有M 头，1——60
岁的大象共有W 头。

设第i 年龄组占1——60岁的比例为i w ，第i 年龄组的
存活率1i s ，则有
则要保持大象群落稳定，如第二问的做法有：
当移出大象头数M=50时，解得b=0.0455。

又由第三问，把这个方程
(1/2)85%X (0.1448-0.0377)=2n 0.144810⨯⨯⨯⨯的0n 换为n 、0.0377换
为0.0455，求得n=1291。

所以，现在只需避孕1291头，由计算机模拟得到的是1297头，非常地接近。

同样的道理可以验证当M=60，70，80…的时候也是与计算机模拟很接近的。

由此，可以说明我们用计算机模拟的方法是有效的。

模型的评价和改进方法
1、模型的优点
（1）本文解决问题的模型都是比较简单的，但是这并不影响得到的结果的准确性，因为这些简单的模型都有很强的理论依据；
（2）在求解第二问的时候，充分利用Leslie矩阵稳定性理论来求解应该让多少母象进行避孕注射，这些理论在差分方程中都是经典的理论，经得起许多事实的考验；（3）第三问的求解中运用了计算机模拟方法来模拟移出大象属于哪个年龄段，这样不仅求解方便、简洁（只需要把算法程序写好就可以得到结果），得到的结果与实际也更接近；
（4）第三问用计算机模拟得到数据后，又用理论去验证，这样使得结果更具有说服力；
2、模型需要改进的地方
（1）因为假设了大象性别是严格地1：1关系，而实际中不一定那么地严格是这样，所以如果能够把各个年龄
段大象的性别比例分别计算，那么模型的结果可能更
接近实际；
（2）在进行计算机模拟时，最开始的随机数的产生个数只有几十个，这几十个随机数不能很好的反映各个年龄
段的大象所占的比重，这样势必会对结果造成一定的
误差；
参考文献：
1、姜启源数学模型（第三版） [M].北京：高等教育出版社 2003
2、赵静数学建模与数学实验（第2版）[M].北京：高等教育出版社 2003
3、周晓阳数学实验与Matlab [M].武汉：华中科技大学出版社 2002
4、郑谏当代数学的若干理论与方法 [M].上海：华东理
工大学出版社 2002
5、李尚志数学建模竞赛教程 [M].江苏：江苏教
育出版社 1996
6、北峰数模网/
2006-8-13
附录：
1、计算机模拟Matlab程序代码
function [bys,b0]=moni(n)
b=[12 16 10 12 17 14];
c=b/sum(b);
b=cumsum(c)';
a=rand(n,1);
k1=0;k2=0;k3=0;k4=0;k5=0;k6=0;
i=1;
for i=1:n
if(a(i)<=b(1))
k1=k1+1;
elseif(a(i)>b(1)&a(i)<=b(2))
k2=k2+1;
elseif(a(i)>b(2)&a(i)<=b(3))
k3=k3+1;
elseif(a(i)>b(3)&a(i)<=b(4))
k4=k4+1;
elseif(a(i)>b(4)&a(i)<=b(5))
k5=k5+1;
else(a(i)>b(5)&a(i)<=b(6))
k6=k6+1;
end
end
s10=1-(11000*0.12*(1-0.9897)+k1)/(11000*0.12); s20=1-(11000*0.16*(1-0.9897)+k2)/(11000*0.16);
s30=1-(11000*0.1*(1-0.9897)+k3)/(11000*0.1);
s40=1-(11000*0.12*(1-0.9897)+k4)/(11000*0.12);
s50=1-(11000*0.17*(1-0.9897)+k5)/(11000*0.17);
s60=1-(11000*0.14*(1-0.9897)+k6)/(11000*0.14);
b0=1/(s10^10*(1+s20+s20^2+s20^3+s20^4+s20^5+s20^6+s2
0^7+s20^8+s20^9)+s10^10*s20^10*(1+s30+s30^2+s30^3+s3
0^4+s30^5+s30^6+s30^7+s30^8+s30^9)+...
s10^10*s20^10*s30^10*(1+s40+s40^2+s40^3+s40^4+s40^5+
s40^6+s40^7+s40^8+s40^9)+s10^10*s20^10*s30^10*s40^10
*(1+s50+s50^2+s50^3+s50^4+s50^5+s50^6+s50^7+s50^8+s5
0^9)+...
s10^10*s20^10*s30^10*s40^10*s50^10*(1+s60+s60^2+s60^
3+s60^4+s60^5+s60^6+s60^7+s60^8+s60^9))*4/3;
bys=0.5*0.85*8864*(0.1448-b0)/(2*0.1448);
2、求大象年龄各个年龄头数及饼图代码
a1=[ 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 1
4 10 16 21 13 10 12...
6 3 6 9 13 10 3 6
21 15 4 13 10 32 14 0];
a2=[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];
a11=sum(a1([1:10]))/sum(a1,2);
a12=sum(a1([11:20]))/sum(a1,2);
a13=sum(a1([21:30]))/sum(a1,2);
a14=sum(a1([31:40]))/sum(a1,2);
a15=sum(a1([41:50]))/sum(a1,2);
a16=sum(a1([51:60]))/sum(a1,2);
a21=sum(a2([1:10]))/sum(a2,2);
a22=sum(a2([11:20]))/sum(a2,2);
a23=sum(a2([21:30]))/sum(a2,2);
a24=sum(a2([31:40]))/sum(a2,2);
a25=sum(a2([41:50]))/sum(a2,2);
a26=sum(a2([51:60]))/sum(a2,2);
c1=(a11+a21)/2*8865
欧阳文创编
c2=(a12+a22)/2*8865
c3=(a13+a23)/2*8865
c4=(a14+a24)/2*8865
c5=(a15+a25)/2*8865
c6=(a16+a26)/2*8865
d=[1091 c1 c2 c3 c4 c5 c6 1044]; pie(d);
欧阳文创编。