Chapter 2 : Model involving logistic growth


2.a- The system of equation and theory
2.b- Numerical exploration
Matlab script

2.a- The system of equations.

We notice that the population growth of both animal is based on a logistic growth which is a bit more realistic than Malthusian growth because involving a carrying capacity Kv and Kp. That permits to the model to include an overpopulation . The natural  demograghy is negatively infuenced for preys by the interaction with the  predator through the  coefficient of predation impact K1 and positively for predators through the predation efficiency coefficient K2.

Stability analysis.

If we are looking for equilibrium states we find :
 
 

As the expression of eigen values for Xe1 is quite complicated we ommited it.

The expression of the jacobian matrix DF(Xe2) and related eigen values  gives :

    
The eigen values are real and the first is always positive if we assume that parameter are positive reals. Thus Xe2 is an unstable equilibrium.
 

As we assume that the parameters are all positive,  equilibrium Xe3 is stable if a < K1Kp .

The eigen values are real positive , this mean Xe4 as an unstable equilibrium .
 
 

BACK
TOP
NEXT

 

2.b-Numerical exploration.
 

Bifurcation from Xe1 to Xe3.
 

-1<a-K1Kf<3  : We can see on the plot bellow the spiral trajectories corresponding to the equilibrium Xe1, the * corresponding to Xe1(a). Xe2=(0,Kf=20)
Note that the point (0,20) is an equilibrium point for Xe1 and Xe3.
 

a=[1, 2, 3, 4, 5]; b=0.5; K1=0.1; K2=2; Kr=2000; Kf=20; R(0)=100; F(0)=10

 

Check if (Kv,0) is really unstable :

F(1)=1.e-8;R(1)=2000;0.1<a<2;b=0.5;K1=0.005;K2=0.005;Kr=2000;Kf=20; -> -2<eigen values (Xe2) <0.5

Here we can notice that all the trajectories converge to the equilibrium 1 (dotted line). Since we do not observe any stabilty for the (Kv,0) equilibrium .

Vary the a parameter : 0.5<a<10

F(1)=10;R(1)=50;b=0.01;K1=0.005;K2=0.005;Kr=2000;Kf=20;dt=0.2s


 
 

Vary the b parameter : 0.001<b<0.5

F(1)=10;R(1)=50;a=0.5;K1=0.005;K2=0.0005;Kr=2000;Kf=20;dt=2s

 

In these numerical exploration we observed only the atraction to the first equilibrium Xe1.
 
 

BACK
TOP
NEXT

 
 

MATLAB  SCRIPT   TO BE COPIED IN A M-FILE
 

%******************************************************
%*PREDATOR-PREY POPULATION DYNAMICS
%*
%*LOTKA VOLTERRA EQUATIONS INVOLVING LOGISTIC GROWTH
%*
%*TIME SCHEME : RUNGE-KUTA OF 4th ORDER
%******************************************************
% PREDATORS=FOXES ->ARRAY F
%PREYS=RABBITS ->ARRAY R
clear all;
help lotg
 

T=150;dt=0.2;
N=T/dt
K=0:dt:T-dt;
a0=0.5;b=0.01;K1=0.005;K2=0.005;Kr=2000;Kf=20;
range=20;Na0=5;
F=zeros(1,N);R=zeros(1,N);
F(1)=10;R(1)=50;
a=(a0):a0*(range-1)/Na0:(a0*range);
[arow,acol]=size(a);

sprintf('SIMULATION TIME :%d',T)
sprintf('TIME STEP FOR RK4 :%f',dt)
sprintf('REQUEST AN AMOUNT OF ITERATIONS OF :%d',N)
sprintf('LOGISTIC GROWTH COEFFICIENT FOR PREY(RABBITS) :%f<%f\n',a(1),a(acol))
sprintf('LOGISTIC GROWTH COEFFICIENT FOR PREDATOR(FOXES) :%f<%f\n',b)
sprintf('CARRYING CAPACITY FOR PREY(RABBITS) :%d\n',Kr)
sprintf('CARRYING CAPACITY COEFFICIENT FOR PREDATOR(FOXES) :%d\n',Kf)
sprintf('PREDATION IMPACT COEFFICIENT FOR PREY(RABBITS) :%f\n',K1)
sprintf('PREDATION IMPACT COEFFICIENT FOR PREDATOR(FOXES) :%f\n',K2)

sprintf('EQUILIBRIUM POPULATION FOR PREY(RABBITS) :%f<R<%f\n',-b*Kr*(K1*Kf-a(1))/(a(1)*b+K1*K2*Kf*Kr),-b*Kr*(K1*Kf-a(acol))/(a(acol)*b+K1*K2*Kf*Kr))
sprintf('EQUILIBRIUM POPULATION FOR PREDATOR(FOXES) :%f<%f\n',a(1)*(1+b*(K1*Kf-a(1))/(a(1)*b+K1*K2*Kf*Kr))/K1,a(acol)*(1+b*(K1*Kf-a(acol))/(a(acol)*b+K1*K2*Kf*Kr))/K1)
 

% LOOP ON a -values
for i=1:acol

  % LOOP ON TIME
   for k=1:N-1

      R1=R(k)+dt*(a(i)*R(k)*(1-R(k)/Kr)-K1*R(k)*F(k));
      F1=F(k)+dt*(b*F(k)*(1-F(k)/Kf)+K2*R(k)*F(k));
      R2=R(k)+dt*(a(i)*R1*(1-R1/Kr)-K1*R1*F1);
      F2=F(k)+dt*(b*F1*(1-F1/Kf)+K2*R1*F1);
      R3=R(k)+dt*(a(i)*R2*(1-R2/Kr)-K1*R2*F2)/2;
      F3=F(k)+dt*(b*F2*(1-F2/Kf)+K2*R2*F2)/2;
      R4=R(k)+dt*(a(i)*R3*(1-R3/Kr)-K1*R3*F3);
      F4=F(k)+dt*(b*F3*(1-F3/Kf)+K2*R3*F3);

      R(k+1)=R(k)+(2*(R2-R(k))+4*(R3-R(k))+2*(R4-R(k))+dt*(a(i)*R4*(1-R4/Kr)-K1*R4*F4))/6;
      F(k+1)=F(k)+(2*(F2-F(k))+4*(F3-F(k))+2*(F4-F(k))+dt*(b*F4*(1-F4/Kf)+K2*R4*F4))/6;

  end
% GRAPHIC OUTPUT
   subplot(1,3,1);
   h=plot(F,R,'b-');
   xlabel('Foxes')
   ylabel('Rabbits')
   hold on
      zoom
   subplot(1,3,2);
   h=plot(F,K,'r-');
   xlabel('Foxes')
   ylabel('Time')
   hold on
zoom
   subplot(1,3,3);
   h=plot(R,K,'b-');
   xlabel('Rabbits')
   ylabel('Time')
zoom
   set(h,'linewidth',1)
   %set(h,'markersize',2)
   hold on
   zoom
end
% line of population convergence
subplot(1,3,1);
x=-b*Kr*(K1*Kf-a)./(a*b+K1*K2*Kf*Kr)
y=a.*(1+b*(K1*Kf-a)./(a*b+K1*K2*Kf*Kr))/K1
plot(y,x,'k-')

zoom