Chapter 1 : The basic model of Lotka-Volterra : 2 species version.


1.a- The system of equation and theory
1.b- The discrete model
1.c- The results
Matlab script

1.a- The system of equations.


The Lotka-Volterra  (LV)  model describes interactions between two species in an ecosystem, a predator and a prey.
Here we note P and V are respectively the number of living predators and victims.
 

Note all the coefficients are positive.

We note that the growth of the prey population is a Malthusian law  (aV <=> exponential growth) . So a is the natural growth rate of the prey population This exponential growth is negatively infuenced by the interaction with the  predator through the  coefficient of predation impact K1.
For the predator we assume that D is the death rate of predator in the absence  of food (preys). That is the predator population would have an exponential decay( exp(-Dt) ) if the preddation would not have a positive impact on the population growth through the predation efficiency coefficient K2.

The system can be written as : 
If we look for balance points we have to solve F(X)=0

Hence :

If look for stability, we calculate the Jacobian matrix of the system :








Then for the first equilibrium point we find the two eigen values :
 
 










Since the real part of eigen values is null we are in a marginal case of the study of linear stabilty. Trajectoty would cycling around equilibrium point.

For the second equilibrium point we find the two eigen values :

We have here two real eigen values of which one is positive, hence this point is unstable.
 
 
 

BACK
TOP
NEXT

1.b- The discrete model

We use a 4th order Runge-Kutta for the time integration scheme. The scheme is very precise against first order Euler scheme for instance.
On the plots below , we can see that with a 1s time step the Euler scheme makes the trajectory to "explode" although Runge-Kuta sheme makes it converge to (20,30) (green point)

D=0.1;a=0.03;K1=0.001;K2=0.005;dt=1s

 

But evenwith dt=0.1s we have the same behavior:

D=0.1;a=0.3;K1=0.001;K2=0.005;dt=0.1s

Zooming :

We can see that the  use of a 4th ordre Runge-Kuta scheme for time discretization permits to save time (vs Euler 1st order) especially because of its precision.

4th orderRunge Kuta scheme :
Let us assume we have a fisrt order ODE system :

We compute the successive internal steps of the scheme:

F(Xn)=(X1-Xn)/dt       ->X1=Xn+dtF(Xn)
F(X1)=2(X2-Xn)/dt    ->X2=Xn+dtF(X1)/2
F(X2)=2(X3-Xn)/dt    ->X3=Xn+dtF(X2)/2
F(X3)=(X4-Xn)/dt      ->X4=Xn+dtF(X3)

Xn+1=Xn+dt(F(X1)+2F(X2)+2F(X3)+F(X4))/6
A quick and uncomplete simplification yields to:
Xn+1=Xn+(2(X2-Xn)+4(X3-Xn)+2(X4-Xn)+dtF(X4))/6

 
 
BACK
TOP
NEXT

 

1.c- The results

For concreteness let us assume that fox is the predator and rabbit the prey. Intial point in phase space is a red * and the equilibrium point corresponding to the parameter set is a cyan o.
 
On this first example we assume aD<1:

F(1)=10;R(1)=50;D=0.1;a=0.3;K1=0.01;K2=0.001;dt=1s
 

We can see that the trajectory in the phase space is cycling around the equilibrium state and converge to this point. This mean that this equilibruim point is stable for such seytof parameters. We can observe on the time history plot that population of both foxes and rabbits follows a quasi-periodic evolution. This behavior is  unusual  in ecosystems but it has been observed that for instance the snowshoe hare and the Canadian lynx population follows this behaviour.
 

Second example aD>1:

F(1)=10;R(1)=50;D=1.3;a=1.3;K1=0.01;K2=0.001;dt=0.1s

It should be noticed that even with a modulus of eigen values grater than 1 the first balance state is still stable.

Third example : aD<1 and initial conditions are close to the second equilibrium i.e. (0,0).

F(1)=1.e-4;R(1)=1.e-3;D=0.1;a=0.3;K1=0.01;K2=0.001;
zoom->
 

Here we can observe that even with an initial point close to (0,0) the trajectory converge to the first equilibrium. This confirm the unstable property of the (0,0) point. Note that the more the intial condition is close to (0,0) the smaller the time step must be set.
 
 
 

 
 
BACK
TOP
NEXT


MATLAB  SCRIPT   TO BE COPIED IN A M-FILE

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

% INPUT : MANUAL OR INTERACTIVE

% INTERACTIVE INPUT :
%T = input('Time lenght of the simulation = ');
%dt = input('Time step = ');

%a = input('a:natural growth rate of rabbits in the abscence of predation=');
%b = input('K1:predation impact coefficient(for prey)=');
%c = input('D:natural deth rate of foxes in abscence of food=');
%e = input('K2:coeffient of efficiency of predation(for predator)=');

% MANUAL INPUT
T=200;
dt=0.05;
T/dt
K=0:dt:T-dt;
F(1)=10;R(1)=50;
D=0.5;a=0.01;K1=0.001;K2=0.005;

sprintf('Equilibrium population for preys :%d\n',D/K2)
sprintf('Equilibrium population for predators :%d\n',a/K1)

% LOOP ON TIME

for k=1:T/dt-1

% RUNGE-KUTA INTEGRATION

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

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

end
 

% SORTIES GRAPHIQUES
subplot(1,2,1)
plot(K,R,'k-');
%title('nb rabbit vs time')
hold on;
plot(K,F,'r-');
title('number Foxes (red) & Rabbits vs time')
xlabel('Time')
ylabel('Population')
hold off
subplot(1,2,2)
plot(R,F,'b-');
hold on
%EQUILIBRIUM STATE=GREEN O
plot(a/K2,D/K1,'go');
%INITIAL STATE=RED *
plot(R(1),F(1),'r*');
xlabel('Rabbits')
ylabel('Foxes')
title('number of Foxes vs Rabbits')
zoom;