24美赛A题食物链捕食者＜-七鳃鳗＜-食物，有性别调节和无性别调节比较（详细解读，有参考文献（真有用）和代码）

# 注：能力不是很多，网上东西确实看不懂，尽自己最大的理解来写了，真是论文翻烂了都，唯一用上的就草-羊-狼这个简单模型，不多说了，这个点还在敲，只能先把自己理解的敲完了，搞一个模型出来，后面啥生态系统分析就GG

1

2

3

4

5结果

代码

%食物链是 捕食者<-七鳃鳗<-食物
%食物资源是 食物数量/其最大环境容量
%R是雄性的概率
%雄性占比和食物资源的关系
% R（x） x 食物资源 R雄性占比
R = @(x) (55.09 + 24.72 / (1 + exp(x*10 - 6.638)))/100 % x = 食物数量/其环境容纳量（最大量）


%雄性占比和七鳃鳗种群增长率的关系
%性别比例是雄个数/雌个数

r0 = 0.5;
% 出生人口性别比例 放在下面欧拉差分方程求解里在
% p = R(x)/(1-R(x));
% r = r0/p;
% r = 0.5*(1-R(x))/(R(x));
%最后得到r与x的关系
%食物资源和种群增长速率的模型
%注：食物资源影响雄性比例，雄性比例影响性别比例，性别比例影响种群增长速率


clc
clear
close all
% 差分动态模型
dt = 0.1;    % 时间步长
T = 10;      % 总模拟时间
n = T / dt;  % 总步数,迭代100代
% sigma1 = 0.6;sigma2 = 5;sigma3 = 0.5;sigma4 = 2;
% r1 = 1;r2 = 0.5;r3 = 0.6;
% N1 = 1000;N2 = 300;N3 = 40;
% f = [r1*x(1)*(1-x(1)/N1-sigma1*x(2)/N2);
%      r2*x(2)*(-1-x(2)/N2+sigma2*x(1)/N1-sigma3*x(3)/N3);
%      r3*x(3)*(-1-x(3)/N3+sigma4*x(2)/N2)];
% 食物
k1 = 1000;
r1 = 1;
a1 = 0.6;
% 七鳃鳗
k2 = 300;

R = @(x) (55.09 + 24.72 / (1 + exp(x*10 - 6.638)))/100;
r2 = 0.5;
% r2 = R

a2 = 5;
a3 = 0.5;
% 捕食者
k3 = 40;
r3 = 0.6;
a4 = 2;

for j = 1:2 % 无调节和有调节比较

    % 初始化种群个数
    N1 = zeros(1,n);N1(1) = 1000;
    N2 = zeros(1,n);N2(1) = 200;
    N3 = zeros(1,n);N3(1) = 20;

    for i = 1:n-1
        N1(i+1) = N1(i) + r1*N1(i)*( 1-N1(i)/k1-a1*N2(i)/k2 )*dt;

        if j==1
        N2(i+1) = N2(i) + r2*N2(i)*( -1-N2(i)/k2+a2*N1(i)/k1-a3*N3(i)/k3 )*dt;
        end
        if j==2
        N2(i+1) = N2(i) + 0.5*(1-R(N1(i)/k1))/(R(N1(i)/k1))*N2(i)*( -1-N2(i)/k2+a2*N1(i)/k1-a3*N3(i)/k3 )*dt;
        end
       
        N3(i+1) = N3(i) + r3*N3(i)*( -1-N3(i)/k3+a4*N2(i)/k2 )*dt;
    end
    figure;
    hold on
    box on
    if j==1
        title('无调节')
    else
        title('有调节')
    end
    plot(1:n,N1(1,:),'g',1:n,N2(1,:),'b',1:n,N3(1,:),'r')
    legend("食物","七鳃鳗","捕食者")
    xlabel("时间")
    ylabel("种群数量")
    hold off
end


%% R = @(x) (55.09 + 24.72 / (1 + exp(x*10 - 6.638)))/100 从这来的
% 使用fit和fittype函数
% doc fittype
clc
clear
close all
% 创建原始数据
% x0 = [4 5 8 10]
% y0 = [78 76 60 56]
x0 = [4 5 8 10]
y0 = [78 76 60 56]
hold on
plot(x0,y0,'r*','LineWidth',3)

% fo = fitoptions('Method','NonlinearLeastSquares',...
%                'Lower',[1,0.56],...
%                'Upper',[0.4,0.78],...
%                'StartPoint',[1,0.56]);
ft = fittype('a + b/(1+exp(x + c))','independent','x');

% 原始数据绘图
% Ftype = fittype('a+b*log(x)+c*y','independent',{'x','y'},'dependent',{'z'})
[f,fitStruct]= fit(x0',y0',ft);
%plot(f,[x0',y0'],z0');
 

xlabel('x');
ylim([56,78])
disp('残差平方和');
fitStruct.sse
% 直接可以使用f(x,y)
f(0.4)
f(0.5)
x1 = 4:0.1:10;
plot(x1,f(x1))
hold off

 参考资料分享（真对自己有用的）

1.Lotka-Volterra 猎食者-猎物模型 - 知乎 (zhihu.com)

 2.https://max.book118.com/html/2017/0716/122502653.shtm

3.基于出生性别比例的逻辑斯谛人口模型的混沌性质 (ctbu.edu.cn) 

注：至于后面，没想好怎么去解读自己得的结果，只能从平稳解和变化趋势上瞎扯一点，网上的那些参考答案，唉（尤其是建模老哥）