cic.m

一个关于数学建模的资料,是纽约大学的程序和文档

% Cicadas and predators. 
% Why do the cicadas emerge all at once?
% Why at prime number of years' intervals?

% nc     Number of species of cicadas. 
% x(j,k) Population of cicadas of species k and age j.
% nx(k)  Life-span of the cicadas of species k.

% np     Number of species of predators. 
% y(j,k) Population of predators of species k and age j.
% ny(k)  Life-span of the predators of species k.

% Initialization:

clear nx ny x y x2 y2 xp yp
clear ax rx dxy ay0 ay1 ry0 ry1

nc=3;
nx=[12 13 15];

x=zeros(max(nx),nc);

np=3;
ny=[2 3 5];

y=zeros(max(ny),np);

% Parameters:

% K      Capacity of the soil.
% ax     Yearly survival rate of the cicadas under earth.   
% rx     Number of offspring per cicada.
% dxy    Number of cicadas eaten by each predator.

% ay=ay0+ay1*x     Yearly survival rate of the predators.
% ry=ry0+ry1*x     Number of offspring per predator.

K=sum(nx);
Nxmin=min(nx);
Nymin=min(ny);
Nymax=max(ny);

% cr How much more do the predators eat before reproducing.

cr=4;
cr=cr-1;

% Approximated expected values of xn and yt:

xnexp=(K*0.9^Nxmin)/Nxmin; 
ytexp=(Nymin+cr)*np/Nymin;

% Initialization with random populations of cicadas and birds,
% with parameters chosen so as to give approximately equal
% chances of survival to all species.

for k=1:nc,
  xmin(k)=0.1;
  for j=1:nx(k),
    x(j,k)=Nxmin*rand/nx(k);
  end
  x(1,k)=max(x(1,k),xmin(k));
  ax(k)=1-0.1*Nxmin/nx(k);
  dxy(k)=0.01;
  rx(k)=ax(k)^(-nx(k))+dxy(k)*nx(k)*ytexp/(Nxmin*xnexp);;
end

for k=1:np
  ymin(k)=0.1;
  for j=1:ny(k),
    y(j,k)=Nymin*rand/ny(k);
  end
  y(1,k)=max(y(1,k),ymin(k));
  ay0(k)=0.75;
  ay1(k)=0.025;
  ry1(k)=0.3*ny(k)/(Nymax*xnexp);
  ry0(k)=(ay0(k)+ay1(k)*xnexp)^(-ny(k))-ry1(k)*xnexp;
end

eps=0.00000001; % A small value, to avoid dividing by zero.

% Final year of the run:

tf=10000;

xp=zeros(1:tf,1:nc); % Variables for plotting.
yp=zeros(1:tf,1:np);

% Main loop.

for t=1:tf,

%   yt: Total number of birds, weighted by how much they eat.

    yt=sum(sum(y));
    for k=1:np,
      yt=yt+cr*y(ny(k),k);
    end

%   xn: total number of adult cicadas.

    xn=0;
    for k=1:nc,
      xn=xn+x(nx(k),k);
    end

% Update of the cicadas' population:

  for k=1:nc,

%   xy * x(nx(k),k) measures the amount of cicadas of
%                   species k which are eaten by the birds.
%                   This number is proportional to yt, the total
%                   number of birds, and to x(nx(k),k)/xn, the
%                   proportion of cicadas of species k among the
%                   total population of adult cicadas.

    xy=dxy(k)*yt/(xn+eps);

    x2(2:nx(k),k)=ax(k)*x(1:nx(k)-1,k);

    x2(1,k)=max((rx(k)*ax(k)-xy)*x(nx(k),k),xmin(k));

  end

% Number of new nymphs that fit under earth:
% Kt is the total number of cicadas; of these,
% only K can stay. 

  Kt=sum(sum(x2));

  if(Kt > K)
    K1=sum(x2(1,:));
    Kavail=K-(Kt-K1);
    prop=Kavail/K1;
    x2(1,:)=prop*x2(1,:);
  end

% Update of the predators' population:
% Their survival rate ay depends on the
% total number of adult cicadas xn.

  for k=1:np,

    ay=min(ay0(k)+ay1(k)*xn,1);
    ry=ry0(k)+ry1(k)*xn;
    
    y2(2:ny(k),k)=ay*y(1:ny(k)-1,k);
    y2(1,k)=max(ry*ay*y(ny(k),k),ymin(k));

  end

  x=x2;
  y=y2;

% Storage:

  xp(t,1:nc)=x(1,1:nc);
  yp(t,1:np)=y(1,1:np);

end

subplot(321)
plot(xp(:,1),'.'); xlabel('t'); ylabel('12 year cicadas');
subplot(323)
plot(xp(:,2),'.'); xlabel('t'); ylabel('13 year cicadas');
subplot(325)
plot(xp(:,3),'.'); xlabel('t'); ylabel('13 year cicadas');

subplot(322)
plot(yp(:,1),'.'); xlabel('t'); ylabel('2 year birds');
subplot(324)
plot(yp(:,2),'.'); xlabel('t'); ylabel('3 year birds');
subplot(326)
plot(yp(:,3),'.'); xlabel('t'); ylabel('5 year birds');