1.a- The system of equation and theory

1.b- The discrete model

1.c- The results

Matlab script

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)=0Hence :

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

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 :

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

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;