蜜蜂数量减少的若干研究数学建模

一、问题重述

1.1 问题背景

蜜蜂的数量对于蜜蜂生产以及生态环境的优化有着极大作用，对于后者主要在于植物的授粉，从而促进植物的繁殖；对于蜜蜂种群的一些基本背景，有以下介绍：蜜蜂可以飞行20公里，但通常停留在离蜂巢6公里以内；一个典型的蜂房包含2万到8万只蜜蜂；一只蜜蜂一天可以访问大约2000朵花甚至更多；由于夏季的繁重工作，大多数蜜蜂都工作到死，导致寿命缩短；蜜蜂的活动水平、花粉消耗和蛋白质含量都会影响它的寿命.此外，对于蜜蜂种群的刻画，有三个基本指标——数量特征、空间特征以及遗传特征，本文将以中华蜜蜂为研究主体，以中华蜜蜂的种群数量特征为研究对象，对三个问题进行解答和分析.

1.2 问题提出

问题一：构建一个数学模型，刻画时间与蜜蜂种群数量之间的关系？换句话说，蜜蜂种群数量是随时间如何变化的？

问题二：在第一题的基础上，进一步的找出对蜜蜂种群有影响的因素，并对影响蜜蜂种群数量的若干因素进行敏感性分析，探究哪一个因素或者哪几个因素对于蜜蜂种群数量的影响最大？

问题三：结合题目给出的蜜蜂活动范围、授粉范围等背景，探究在81000平方米的土地，至少需要多少蜂箱，才能够保证所有蜜蜂均有栖息地。

---

二、问题分析

2.1 对问题一的分析

问题一总共是一个小问，目的是探究时间与蜜蜂种群数量之间的关系。对于这个模型，在查找相关的资料后，在不同时间段中华蜜蜂的出生率、死亡率、迁入率、迁出率，利用简单的生物学模型，建立种群数量与时间之间的常微分方程模型进行求解，在求解过程中分成三种情况——实验蜜蜂，包括在一定时间内未达到环境容纳量以及达到环境容纳量，此时的模型采用单种群增长模型，分别采用Malthus（马尔萨斯）模型以及Logistic(逻辑斯蒂)模型，建立常微分方程进行求解；以及第三张情况，就是对于一般的种群，生活环境为自然环境，此时考虑多种群增长模型，建立Leslie(李斯列)人口模型，建立差分方程组进行求解.

2.2 对问题二的分析

问题二总共是一个小问，目的是探究对蜜蜂种群数量造成影响的主要因素。对于这个问题，我所采用的解决方法是借助第一题得到的时间与蜜蜂种群数量之间的关系，特别是对于建立的差分方程模型，分别寻找对蜜蜂出生率、死亡率、迁入率、迁出率造成影响的因素，然后利用敏感性分析模型，对这些影响因素进行相关性分析并得到结论。

2.3 对问题三的分析

问题三主要是研究至少需要的蜂箱的数量，使得能支持一个20英亩的授粉 (81000平方米)的土地种植受益于授粉的作物。对于此问题，根据题目所给的若干条件——比如蜜蜂可以飞行20公里，但通常停留在离蜂巢6公里以内等等，并通过对蜂箱的规模和每平方米的话的数量进行假设，利用几何分析的方法对问题进行求解。

---

三、模型假设

3.1 假设所研究的蜜蜂种群为典型的中华蜜蜂；

3.2 假设不同时间段的蜜蜂的年龄组配的每一项的取值区间（a1,b1），则取这一项的比例为；

3.3 假设第一个模型的第一种情况中蜜蜂群体是在封闭实验室中进行，除了天敌因素之外，其他条件与自然环境完全一致（即生活环境与自然环境几乎完全一致），因此可以忽略迁入率与迁出率因素；

3.4 假设实验室的空间和其他因素等条件能够使得蜜蜂在一年的时间段内不会达到环境容纳量K;此外，由生物学知识可知在一定条件下环境容纳量几乎为定值，因此假设在实验室中的环境容纳量K和在自然环境中是相等的。

3.5 假设一个月有30天；

3.6 假设在理想条件下蜜蜂的孵化率为1，即出生率＝生育率；

3.7 假设每一个人工蜂箱能够容纳10000只蜜蜂；

---

四、符号说明

---

五、模型分析与求解

5.1问题一分析与求解

针对本题，首先找出不同月段的蜜蜂卵、幼虫、蛹、成虫的年龄组配并进行数据统计，得到的数据如表1

表1 不同月段的蜜蜂年龄组配

5.1.1 模型1分析与求解

首先，对于第一种情况，如果研究的蜜蜂群体是属于实验蜂，也就是生活在封闭实验室中的蜜蜂，且模型建立在实验时间为1年的前提下，此时建立简单的常微分方程模型即可求解问题。由于在这种情况下，根据假设3，蜜蜂的种群数量仅受出生率和死亡率的影响，因此只需要根据实验过程中得到的中华蜜蜂在不同时间段内的出生率和死亡率建立微分方程模型，再根据假设4，可建立Malthus（马尔萨斯）模型：

对任意的t0>0，有

（1）

同时已知N0,解得：

  （2）

由此可以得到在一年的时间内，蜜蜂种群数量N和时间之间的关系式为.根据所查阅的资料：，由此可以解得

（t为天数）。

这是基于最理想的状态，也就是蜜蜂的出生率和死亡率十分稳定且种群数量在一定时间段内没有达到环境容纳量K，但根据生物学知识，这种情况几乎不会发生。因此，我考虑了更加细致的出生率和死亡率，也就是将1年内的不同的四个时间段的出生率和死亡率，则:

（3）

5.1.2 模型2分析与求解

进一步的，若在一定时间内蜜蜂种群达到了环境容纳量K，则此时不仅考虑出生率与死亡率，还要进一步考虑固有增长率r，并根据Logistic模型建立常微分方程进行求解：

（4）

利用Matlab求解（4），解得：

（5）

进一步的，利用Matlab对（5）进行绘图，如图一所示（分别对K,N0,r赋特殊值）：

图一

可见对于实验室环境的进一步考虑下，通过建立Logistic模型进行求解，与

Malthus模型相比，蜜蜂种群数量数量在达到环境容纳量K之后，并不会继续无限制的增长，因此更加符合自然条件下的蜜蜂种群增长情况。而当初始值大于K时，则N在t大于零到某一个T之间恒为减函数，直到N（t）稳定在N附近，在一定的波动范围内。

进一步地，若考虑不同年龄段的蜜蜂的年龄组配与相对应的出生率与死亡率，可以对（5）进行进一步的改进，即考虑考虑不同年龄段的蜜蜂的数量,则可以的到以下的模型：

（6）

其中，

(7)

5.1.3 模型3分析与求解

以上两种情况均基于假设3，但如果对于生活环境为一般自然环境的蜜蜂而言，除了考虑出生率和死亡率，还需要考虑迁入率与迁出率，此时，以上两个简单模型并不能够用于解决这一复杂问题，影响种群数量的因素还有迁入率与迁出率。利用Leslie模型建立差分方程组，对着一问题进行求解：

如上，将蜜蜂分为4个不同年龄阶段——卵、幼虫、蛹、成虫，该模型主要讨论在不同时间内不同年龄段的蜜蜂的分布，对时间进行离散化为t=0，1，2...，定义向量

（i=1..4）（8）

则满足以下的差分方程组：

（9）

定义矩阵L：

（10）

则（9）可以表示为

N(t 1)=LN(t)（11）

从而当L时良定义的时候，则只需要根据初始值N0，便可以求解差分方程（10），从而求解N（t）。

此外，若进一步考虑环境因素，我们可以假设在本地区新一天t 1中的蜜蜂，其主要来源有三个：

（1）第t天留下来的蜜蜂

（2）第t 1天新出生的蜜蜂

（3）第t天净迁入的蜜蜂

利用以上三个关系，可以建立以下的差分方程模型

设第t天的蜜蜂数量为N（t），则由假设条件与相关符号可得：第t 1天的蜜蜂数量变为，则蜜蜂种群数量变化方程的离散形式如下：

（12）

其中是第t天中雌性可生育蜜蜂的比例，区间[a，b]是指蜜蜂的可生育年龄段；则（12）揭示了不同年龄段、雌性蜜蜂比例、时间等对于蜜蜂数量变化的一般规律.

为了简化方程组，我们引入下面的若干记号：

5.1.4 对第一小题的解题心得

第一题只需要找出蜜蜂种群数量与时间之间的关系，因此在解决该题的过程中，我分成了三种情况，按照有简单到复杂，有单一因素（Malthus模型）分析到多因素分析（Lislie模型）.对于此题，我做了很多假设，目的是建模过程中可以简化解决过程.

不过对于此题，我认为题目可以对蜜蜂的种群类型和生活地区进行进一步详尽的说明或者提示，因为在解决该题的过程中，我发现了一个问题，首先是对于同种蜜蜂（例如我选取的中华蜜蜂），在不同地区的不同个体的年龄组配、出生率等等都有一定的差别，因此会造成建模过程中情况十分复杂，需要分若干种情况进行说明，这样子会使得模型十分复杂，如果要解决真正的实际问题.因此在这个问题中，我并没有查找到真实的材料，而只是利用各种符号代表了各个指标，在指标的基础上进行建模.

第二个我觉得需要注意的地方是，“一段时间内”这个表述其实是比较模糊的，因为我们可以选取一天，或者一礼拜，一个月，一个季度等等，而在不同的时间段内，蜜蜂的数量变化规律是有所不同的，举个最简单的生物学知识，尽管自然环境中存在着一个常数K，定义为环境容纳量，这个数也许很大，对于一个初始数量远小于K的蜜蜂种群，在短时间内达到K是几乎不可能的，但是当时间区间T足够大的时候，利用简单的Logistic阻滞增长模型可以知道，一定存在某一个t在时间区间T中，使得蜜蜂种群数量会达到K，那么在接下来的时间段内，蜜蜂种群数量将会在K值附近进行合理的振荡而不会有较大的误差，所以在解决该问题的时候，我是分成了若干种情况进行分析，构建的模型也有所差异.

5.2问题二分析与求解

5.2.1 第二小题分析与求解

针对本题，由于题目所给的数据并不足以进行敏感性分析，因此首先是查阅了相关的资料，例如蜜蜂数量越多，蜂巢内温度越稳定，并能够保持在适温34摄氏度-35摄氏度的范围内，温度太高或太低均不利于其发育和繁殖；我是基于第一题所得到的差分方程模型，对N（t）和不同时期的蜜蜂的年龄组配以及迁入迁出率、不同时间段不同年龄组配的个体的出生率与死亡率进行敏感性分析.

我们还是分为两种情况进行分析——分别对Logistic阻滞增长模型和Leslie模型进行敏感性分析.这里不讨论Malthus是因为这是Logistic模型的一种特殊形式（K无限大），因此为了简化分析，只讨论第二种和第三种模型对的敏感性分析.

首先是对于Logistic阻滞增长模型的分析，由（6）

我们分别对自然增长率、固有增长率、环境容纳量K以及每个年龄组配的比例进行敏感性分析.由K是一个常数，故有：

（14）

由此可以直观看出K与N是成正比的，且N对K的偏导数是一个小于等于1的数，对任意的K，所以可以得到的第一个敏感性分析结论：当实验时环境或者自然环境的条件——包含温度、湿度、栖息地天敌数量等等均有利于蜜蜂的生存时，环境容纳量K会增大，则在短时间内，数量与环境容纳量之间会成正比，也就是环境容纳量对于N（t）是敏感的.

对于年龄组配的敏感性分析，我们采用直接求导函数的方法，也就是对求解导函数，而这通过Matlab计算得到以下的结果：

（15）

因此，N对的偏导数是一个大于1的正数，且随着时间增大，这个偏导数会变大，也就是随着时间增大，年龄组配对于蜜蜂种群数量的影响程度是成正相关关系的，也就是年龄组配对蜜蜂种群数量是敏感的.而且还能够得出年龄组配对蜜蜂种群数量的影响程度要强于环境因素.

另外，对于固有增长率也进行敏感性分析，有：

ans(a, b, c, d) =(K*a*t*exp(log(K/N0) - r*t))/(exp(log(K/N0) - r*t) 1)^2 (K*b*t*exp(log(K/N0) - r*t))/(exp(log(K/N0) - r*t) 1)^2(K*c*t*exp(log(K/N0) - r*t))/(exp(log(K/N0) - r*t) 1)^2(K*d*t*exp(log(K/N0) - r*t))/(exp(log(K/N0) - r*t) 1)^2

由于0＜r＜1，因此可以认为在时间T范围内，当r接近log(K/N0)/T时，种群数量与r之间的相关性是最强的，且当r小于log(K/N0)/T时，N与r之间是负相关关系，当r大于log(K/N0)/T时，N与r之间是正相关关系.

此外，有第一个问题还可以知道，N除了与以上因素相关，还与迁入迁出率相关，现对Leslie模型进行敏感性分析.由第一问可以得到种群数量N（t）与时间之间存在着这样的差分方程关系（13）

我们只需要对进行敏感性分析。显然的，此时N对求偏导数，得到的结果恒为1。

综上所述，年龄组配、迁入迁出率、固有增长率、环境容纳量（环境），对于N(t)影响最大的是固有增长率和年龄组配，其次是迁入迁出率，再者是环境容纳量，也就是环境因素.

5.2.2 对第二小题的解题心得

完成了第二小题的解答，我觉得较大的困难在于，题目缺乏各个指标的数据，所以没办法对每一个影响因素进行数据拟合，从而进一步优化第一小题得到的模型。因为第一小题并没能够细化到每一个因素，而是只抓住了影响种群数量的四个因素——迁入率、迁出率、出生死亡率进行建模和求解，因此在第二小题，利用第一小题得到的模型进行敏感性分析其实是有一定难度的，因为完全是利用代数式进行分析，而缺少确切的数据进行验证.

此外，对于我的做法，我觉得还是有很多不足，例如在分析K的敏感性时，其实可以进一步细化的，例如可以继续分为天敌、温度、湿度等等因素进行分析，这样子能够更加具体的解答题目；因此，对于此题的思路我有以下的一些想法：首先利用已知的数据，对每一个因素与时间之间进行拟合，得到时间与每一个因素之间的函数关系，再利用第一题得到的模型，结合已知的数据进行敏感性分析.

5.3 问题三分析与求解

5.3.1 第三小题分析与求解

针对本题，首先是假设了土地上种植的花卉为虞美人，因为资料显示蜜蜂对于虞美人比较敏感，因此选取虞美人进行假设是合理的。找到资料可以发现，一般地，一平方米的土地可以种植16朵花，并且两两之间均是0.3m*0.3m的间隔；我多采用的是对81000平方米的土地的规格、土地上的花卉的种类以及种植密度进行了一系列假设与模型求解，以下是具体的解答过程.

---

首先证明一个简单的几何命题（最小圆覆盖定理），如图二：

图二

Theorem:对于圆，可以用4个圆进行完全覆盖

在81000平方米的土地中，16朵花的占地面积为0.81平方米，则这片土地中最多可以种植100000朵花。根据假设，由每只蜜蜂最多只能够传授1000朵花，而在每一个蜂巢的影响区域

（定义为蜂巢的影响区域）中种植的花的数量是（18852000/0.81）/16=14546296441，因此在这片区域中至少需要有14546296441/1000=14546297只蜜蜂。由假设3.7，则每一个蜂巢的影响区域内，至多需要的蜂箱的数量为1455个。

利用这个原理对各个假设进行模型求解。首先假设该区域的花卉分布十分特殊，即花卉全部均分布在蜜蜂的有效授粉范围，也就是上，则利用最小圆覆盖定理可以很快就得：

这样子的圆面闭包一共有1个，则最多需要1455个；

此外，可以假设这片土地是宽小于30cm的长方形土地，不妨设长方形的且花朵是在一条直线上，则在内部一个蜂巢x的影响范围是[x-6,x 6]，由已知条件可知这片土地的长为270000m，则花朵的数量最多为270000/0.3 1=900001，因此每一个蜂巢的影响区间上的花的数量最多为12/0.3 1=41，则最多需要900001/41=21952个区间才能够完全覆盖这片土地；此外，每个区间上最多需要21952/10000=3只蜜蜂，就可以实现花朵被完全授粉，因此至多需要21952个蜂箱即可解决问题。

事实上，为了更好的解决以上问题，可以利用贪心算法解决上述问题，具体的代码见附录三.

5.3.2 对第三小题的解题心得

第三小题完全就是一个最优解问题，我认为最好的解决方法是利用图论知识进行建模，由于自己图论知识还只是停留在极少数理论上，没能够深入理解其应用性，所以我还是利用了比较简单的几何分析的方法进行假设求解，从求解的过程看，其实土地的类型和花卉的密度均会影响蜂箱数量的求解结果，因此做好各种假设和预测是非常有必要的。

此外，对于利用图论知识进行建模，我觉得尽管非常高效且准确性可能更高，但是对于普通高中生（指的是广大普通高中，甚至是极少数重点高中的平行班），在高中三年以及之前的学习中并无接触过图论和组合数学，因此采用几何分析的方法解决该问题也是最自然的做法。

---

附录

1.附录一

第一问Matlab代码

1.  
   dsolve('DN=r*N*(1-N/K)','N(0)=N0','t') ;
2.  
    
3.  
   x = -10:0.01:10;y = 1./(1 exp(-x));plot(x,y,'r','linewidth',1.5);xlabel('x');
4.  
    
5.  
   figure(1);x = -10:0.01:10;y = logistic(x);
6.  
    
7.  
   plot(x,y,'r','linewidth',1.5);xlabel('x');hold on;
8.  
    
9.  
   line([0 0],[0 1],'linestyle','--', 'Color','g', 'LineWidth', 1.5);hold on;
10.  
     
11.  
    line([-10 10],[0 1], 'Color','k', 'LineWidth', 1.5);legend('sigmoid','x = 0','y = x');

2.附录二

第二问代码

1.  
   syms a b c d K N0 r t;
2.  
    
3.  
   N(a,b,c,d)=a*K./(exp(log(K./N0)-r*t) 1) b*K./(exp(log(K./N0)-r*t) 1)...
4.  
    
5.  
   c*K./(exp(log(K./N0)-r*t) 1) d*K./(exp(log(K./N0)-r*t) 1);
6.  
    
7.  
   diff(N,a)
8.  
    
9.  
   diff(N,b)
10.  
     
11.  
    diff(N,c)
12.  
     
13.  
    diff(N,d)
14.  
     
15.  
    diff(N,r)

3.附录三

1.  
   function min_ball
2.  
    
3.  
   syms a b c;
4.  
    
5.  
   ooo=[a b c];
6.  
    
7.  
   o=[0 0 0]; %球心
8.  
    
9.  
   R=6; %球的半径
10.  
     
11.  
    n=input('请输入你想要生成的离散点个数:');
12.  
     
13.  
    t=input('请输入你想要生成的点的范围(上限）:');
14.  
     
15.  
    A=rand(n,3);
16.  
     
17.  
    B=(2*A-1)*t;
18.  
     
19.  
    %随机产生n个离散点，储存在pot里面
20.  
     
21.  
    pot=fix(B) 1;
22.  
     
23.  
     
24.  
     
25.  
    %生成最小球
26.  
     
27.  
    for i=1:n
28.  
     
29.  
    if norm(pot(i,:)-o)>R
30.  
     
31.  
    o = pot(i,:);
32.  
     
33.  
    R=0;
34.  
     
35.  
    for j=1:i-1
36.  
     
37.  
    if norm(pot(j,:)-o)>R
38.  
     
39.  
    o = (pot(i,:) pot(j,:))/2;
40.  
     
41.  
    R=norm(pot(j,:)-o);
42.  
     
43.  
    for k=1:j-1
44.  
     
45.  
    if norm(pot(k,:)-o)>R
46.  
     
47.  
    %最小外接球(补充ball函数）
48.  
     
49.  
    �ll函数输入三个点坐标，返回这三个点的最小球的球心坐标
50.  
     
51.  
    x=pot(i,:);
52.  
     
53.  
    y=pot(j,:);
54.  
     
55.  
    z=pot(k,:);
56.  
     
57.  
    if norm(z-y)^2 4*R^2<=(norm(z-x))^2
58.  
     
59.  
    ooo=(z x)/2;
60.  
     
61.  
    R=norm(z-ooo);
62.  
     
63.  
    elseif norm(z-x)^2 4*R^2<=(norm(z-y))^2
64.  
     
65.  
    ooo=(z y)/2;
66.  
     
67.  
    R=norm(z-ooo);
68.  
     
69.  
    else %三个点构成锐角三角形的情况
70.  
     
71.  
    ooo=ballcenter(x,y,z);
72.  
     
73.  
    R=norm(z-ooo);
74.  
     
75.  
    o = ooo;
76.  
     
77.  
    R=norm(pot(k,:)-o);
78.  
     
79.  
    end
80.  
     
81.  
    end
82.  
     
83.  
    end
84.  
     
85.  
    end
86.  
     
87.  
    end
88.  
     
89.  
    end
90.  
     
91.  
    end
92.  
     
93.  
     
94.  
     
95.  
     
96.  
     
97.  
    o=double(o); %保证小数防止mesh出bug
98.  
     
99.  
    R=double(R);
100.  
      
101.  
     [q, w, e]=sphere(30);
102.  
      
103.  
     Q=R*q o(1);
104.  
      
105.  
     W=R*w o(2);
106.  
      
107.  
     E=R*e o(3);
108.  
      
109.  
     subplot(1,3,1)
110.  
      
111.  
     mesh(Q,W,E)
112.  
      
113.  
      
114.  
      
115.  
     hold on
116.  
      
117.  
     %画n个点
118.  
      
119.  
     x1=pot(:,1);
120.  
      
121.  
     y1=pot(:,2);
122.  
      
123.  
     z1=pot(:,3);
124.  
      
125.  
     subplot(1,3,2)
126.  
      
127.  
     scatter3(x1,y1,z1,'r');
128.  
      
129.  
      
130.  
      
131.  
     hold on
132.  
      
133.  
     subplot(1,3,3)
134.  
      
135.  
     mesh(Q,W,E);
136.  
      
137.  
     alpha(0.8) %设置透明度
138.  
      
139.  
     shading flat %去掉那道些线
140.  
      
141.  
     hold on
142.  
      
143.  
     scatter3(x1,y1,z1,'r');

这篇好文章是转载于：学新通技术网