In: Advanced Math
this is my Matlab code for lemur and plant population simulation. But i didn't get the right graph.
R(1,1)=35; % Initial pop of males
R(2,1)=35; % Initial pop of females
P(3,1)=200000; % Iniitial pop of plants
r(1)=1.01; % Growth rate of male pop
r(2)=1.01; % Growth rate of female pop
r(3)=1.10; % Growth rate of plants
K(1)=800; % Carrying capacity for male pop
K(2)=600; % carrying capacity for female pop
K(3)=7500000; % carrying capacity for plants
E=65; % Plants per animal
A=150000; % Animal 'Plant needs'
NW=200; % Weeks for simulation
for m=2:NW
P(3,m)=(r(3)*P(3,m-1)*((K(3)-P(3,m-1))/K(3))) -(E*R(m));
R(m-1)=R(1,m-1)+ R(2,m-1);
R(1,m)=(r(1)*R(1,m-1)*(K(1)-R(1,m-1))*P(3,m-1))/(K(1)*A);
R(2,m)=(r(2)*R(2,m-1)*(K(2)-R(2,m-1))*P(3,m-1))/(K(2)*A);
end
subplot(1,2,1)
plot(1:NW,R(1,1:NW),1:NW,R(2,1:NW))
xlabel('Month')
ylabel('population')
legend('male','female')
subplot(1,2,2)
plot(1:NW,P(3,1:NW))
xlabel('Month')
ylabel('Plant population')
i guessed there are two problems with this code.
(i) There is direct assignment to the third row of P matrix. The first two rows of P matrix are zeros. So there is no need of the first two rows. Instead of P(3,m) you can use R(3,m).
(ii) How you are updating and using R(m) and R(m-1) ?? ... Doing R(m) in MATLAB considers R(m,1) by default. If you mean this, then it is right... otherwise you've to look into this.
However, I've corrected the first issue and the code is listed below :
---------------
R(1,1)=35; % Initial pop of males
R(2,1)=35; % Initial pop of females
R(3,1)=200000; % Initial pop of plants
r(1)=1.01; % Growth rate of male pop
r(2)=1.01; % Growth rate of female pop
r(3)=1.10; % Growth rate of plants
K(1)=800; % Carrying capacity for male pop
K(2)=600; % carrying capacity for female pop
K(3)=7500000; % carrying capacity for plants
E=65; % Plants per animal
A=150000; % Animal 'Plant needs'
NW=200; % Weeks for simulation
for m=2:NW
R(3,m)=(r(3)*R(3,m-1)*((K(3)-R(3,m-1))/K(3))) -(E*R(m));
R(m-1)=R(1,m-1)+ R(2,m-1);
R(1,m)=(r(1)*R(1,m-1)*(K(1)-R(1,m-1))*R(3,m-1))/(K(1)*A);
R(2,m)=(r(2)*R(2,m-1)*(K(2)-R(2,m-1))*R(3,m-1))/(K(2)*A);
end
subplot(1,2,1)
plot(1:NW,R(1,1:NW),1:NW,R(2,1:NW))
xlabel('Month')
ylabel('population')
legend('male','female')
subplot(1,2,2)
plot(1:NW,R(3,1:NW))
xlabel('Month')
ylabel('Plant population')
------------------------------
I hope you will find it useful.
Thank You!