2024 MCM数学建模美赛2024年A题复盘，思路与经验分享：资源可用性与性别比例 | 七鳃鳗性别比例变化对生态系统的影响（二）

本文链接：https://blog.csdn.net/2302_80902795/article/details/136290979

书接上回

在上一篇，我们得出了思路，利用logistics增长模型，来建立七鳃鳗性别比例对生态系统的影响模型。那我们就要先知道logistics模型长什么样的，是干什么用的。我这里简单介绍一下，如果想深入了解，你们可以自己去搜索了解。

2024 MCM数学建模美赛2024年A题复盘，思路与经验分享：资源可用性与性别比例 | 审题与选题（一）-CSDN博客

logistics增长模型

logistics模型是一个微分方程，这是logistics模型的最基础形态，N(t)是种群数量，r是种群的自然增长率，t是时间。公式左边dN/dt是种群数量随时间的变化趋势，右边是自然增长率乘以当前种群数量。这个公式合起来就是说种群数量随时间的变化趋势等于自然增长率乘以当前种群数量。然后，我们通过这个微分方程，把N(t)解出来。这个N(t)对应的曲线就是，种群数量的增长曲线。大家不要看到微分方程就害怕，实际上，这个微分方程不需要我们自己解，我们都是通过代码和软件让计算机求解的。

如果你们还有高中生物的知识，那么这个N(t)对应的曲线就是，高中生物说的种群增长的J型曲线。长这样。

有J型曲线，必然有S型曲线对吧。也就是说这个种群不能无限的增长，它还是受到环境的限制。有一个环境允许的最大种群数量。所以我们的公式就变成了这样。

这里多了一个参数Q，Q就是最大种群数量。这样一改，这个微分方程解出来就是S型曲线了。长这样。

这就是logistics增长模型的基础形态了，我们看看养龙那篇论文是如何对这个模型进行修改的。

龙羊生态系统模型

这个模型，上来先假设，在一个只有羊的生态系统中，羊的种群数量在一开始已经达到了最大值。也是就是N0=K。

然后就是，龙作为入侵者进入这个生态系统，对羊的种群数量进行了消耗。然后把公式写成这样。

这个公式就是基础的logistics增长模型减了一个a，a就是龙每天要吃掉羊的数量。这个a他们是怎么算出来的呢。就是他们估计一头龙一天要消耗900公斤肉，而一头羊的体重在100公斤，所以a就等于9。

然后，他们联立这两条公式，就求出来羊的种群数量的变化。

这里，他们假设羊的一开始的数量K=5000，在龙进入生态系统后，羊的数量稳定在了3500。

七鳃鳗-湖鳟生态系统模型

仿照上面那个龙羊生态系统模型，我们同样建立一个七鳃鳗-湖鳟模型，湖鳟是七鳃鳗的食物。同样是在一个只有湖鳟的生态系统中，湖鳟已经达到最大值，即N0=K。然后，我们还考虑了湖鳟的自然死亡d*N(t)。d是自然死亡率。A是七鳃鳗一天要消耗的湖鳟鱼。

这里的重点是七鳃鳗的性别比例是如何影响湖鳟鱼的数量的。也就是，七鳃鳗的雄性和雌性有什么区别。

通过查资料我们发现，一般来说雄性七鳃鳗比雌性七鳃鳗的体重要小一些。这就导致了雄性七鳃鳗消耗的资源更小，一天要吃掉的湖鳟就更少。而雌性要消耗的资源更多，一天要吃的湖鳟就更多。这里我们定义alpha为雄性七鳃鳗的比例。打个比方，如果alpha=0.5，就是雄性七鳃鳗的数量占总体的50%，也就是说雄性比雌性的比例为1：1。如果alpha=0.7，就是雄性七鳃鳗的数量占总体的70%，也就是说雄性比雌性的比例为7：3。

再说回这个A，A就是一个七鳃鳗种群一天要吃掉的湖鳟数量。因为七鳃鳗存在性别差异，所以如果种群中雄性七鳃鳗的数量更多，那么这个七鳃鳗种群一天要消耗的湖鳟就更少。反之，如果雌性七鳃鳗比较多，那么一天要消耗的湖鳟就更多。这里我们不考虑，雌性比雄性多的情况。也就是alpha>=0.5。

那么我们一个雄性七鳃鳗和一个雌性七鳃鳗一天要摄入具体的量，就能计算出这个A了。

$A = 0.023\alpha\cdot 100 + 0.029(1-\alpha) \cdot 100$

这里，我们定义雄性七鳃鳗一天的消耗量是0.023，雌性为0.029。七鳃鳗的种群数量为100。通过这个公式我们就能计算出一个数量为100的七鳃鳗种群一天要消耗的湖鳟数量了。

接着，定义湖鳟鱼的初始数量为1000，湖鳟鱼种群增长率为0.03，死亡率为0.01。联立两个方程，就能求出湖鳟鱼的变化曲线了。而只要改变alpha就能得到再性别比例不同的情况下，湖鳟鱼的变化曲线。这样也就可以性别比例对生态系统的影响了。

这个就是湖鳟鱼种群数量的变化曲线，这种趋于稳定，而alpha的值越大，最后稳定的值就高，也符合我们之前的预期。

这张图可以看的更清楚些。

可看到好像最后的结果比较接近，那我们可以得出，性别比例对生态系统的影响不大吗。是不能的。当我们把初始种群数量设置在700的时候，我们就得到了下面这幅图。

这时，由于alpha的不同，最后的结果就很不相同。这说明当K=1000时，生态相同比较稳定，性别比例对生态系统有影响，但影响比较小。而当K=700时，这时生态系统相比之前就更脆弱了。那么这时，性别比例对生态系统的影响就更大。

代码

这里我们使用matlab的ode45来求解这个微分方程。不了解怎么使用的同学可以自己搜索学习，我就不解释了。

% 参数定义
r = 0.03; % 增长率
K = 700; % 环境承载能力
N0 = K; % 初始种群数量
d = 0.01;%死亡率
alpha = 0.5;
alpha_lis = [alpha];
N_l = 100;%七鳃鳗数量
a = 0.023*alpha*N_l + 0.029 * (1-alpha) * N_l;%捕食量
w = 4;
W = a/w;
% 时间跨度
tspan = [0, 700]; % 从0到100个时间单位

% 定义微分方程
%odefun = @(t, N) r*N(1) * (1 - N(1)/K) - a - d*N;

for i=1:6
    % 定义微分方程
    odefun = @(t, N) r*N(1) * (1 - N(1)/K) - a - d*N;
    % 使用ode45求解微分方程
    [t, N(:,i)] = ode45(odefun, tspan, N0);
    alpha = alpha + 0.07;
    alpha_lis = [alpha_lis,alpha];
    a = 0.023*alpha*N_l + 0.029 * (1-alpha) * N_l;
    W = a/w;
end
% 绘制结果
figure
plot(t, N);
xlabel('Time');
ylabel('Population N');
title('Population Dynamics');
%text('K=1000','FontSize','right');

% 在曲线末尾的x值处添加文本
text(t(end), N(end,1), ['N= ', num2str(N(end,1)),'   ','K=',num2str(K)], 'HorizontalAlignment', 'right');

yline(0,'r')
% 设置图例位置
legend(['alpha = ',num2str(alpha_lis(1))], ...
   ['alpha = ',num2str(alpha_lis(2))], ...
   ['alpha = ',num2str(alpha_lis(3))], ...
   ['alpha = ',num2str(alpha_lis(4))], ...
   ['alpha = ',num2str(alpha_lis(5))], ...
    ['alpha = ',num2str(alpha_lis(6))]);
%grid on;

figure
plot(alpha_lis(:,1:6),(700-N(end,:))/700,'r-*');
title('Stability of salmon population after lamprey invasion')
ylabel('Salmon population')
xlabel('sex ratio')