真菌群落模型

clc;clear
y1=@(W,W0)-0.0171.*W.^4-0.2146.*W.^3-0.9226.*W.^2-1.4116.*W-0.1605+1.7135.*(W0-0.78);%-0.5到-3
y2=@(T,T0)(T-T0+25.47).^4.*1e-6-0.0001.*(T-T0+25.47).^3+0.0021.*(T-T0+25.47).^2+0.0054.*(T-T0+25.47)+0.1075;%10-30
x=@(t,r,K,x0)K*x0*exp(r*t)/(K+x0*(exp(r*t)-1))
d=0.02;%时间步长
long=365%时间长度
t=0:d:long

W=-1.5%模拟环境潮湿度

n=10;
big=200
W0=-0.5-rand(1,200)%最佳潮湿度
T0=20+10.*rand(1,200)%最佳温度
x0=5.*rand(1,200)%初始个体数
K=1000.*rand(1,200)%模拟环境容纳量



z1=0;increase=[];p_=x0';
for i=t
    T=28+8.*sin(2*pi*i);
    r1=y1(W,W0);r2=y2(T,T0);
    r=r1.^2+r2.^2;r=r.^0.5+rand(1,1);%计算增长值
    alpha=zeros(1,big);
    a=1;
    for ii=r
        b=1;
        for jj=r
            alpha(a,b)=r(b)/r(a);
            b=b+1;
        end
        a=a+1;
    end
    z1=z1+1;%时间列
    z2=0;%种群行
   
       
        for k=1:1:big%得到增长率矩阵
            l=0;
            for m=1:1:big
                l=l-alpha(k,m)*p_(m,z1)/K(m);
            end
            
            increase(k,z1)=r(k)*p_(k,z1)*(1+l);
        end

        pt=(increase(:,z1).*d+1).*p_(:,z1);
        
        p_=[p_ pt];
    
end
p_(:,1)=[]


for i=1:1:big
    plot(t,p_(i,:))
    hold on
end
xlabel('/day')
ylabel('number')
title('Population size changes with time')
P=sum(p_);
hold off
plot(t,P)
xlabel('/day')
ylabel('number')
title('Community size changes with time')
su=[];a=0;
for i=P
    su=[su a];
    a=a+i*d;
end

plot(t,su)