# 注：能力不是很多，网上东西确实看不懂，尽自己最大的理解来写了，真是论文翻烂了都，唯一用上的就草-羊-狼这个简单模型，不多说了，这个点还在敲，只能先把自己理解的敲完了，搞一个模型出来，后面啥生态系统分析就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)
注：至于后面，没想好怎么去解读自己得的结果，只能从平稳解和变化趋势上瞎扯一点，网上的那些参考答案，唉（尤其是建模老哥）