% Simulation du modele Lotka Volterra (Murray p63)
% Integration par la methode de Runge Kutta
% On peut discuter l'influence de alpha et gamma
%

clf;
hold on;

U='';
V='';

% Définition du pas de temps
% Constantes de dynamique intrinsèque et populations initiales

dt = 0.01;
Tmax=100;
alpha=1;
gamma=1;
u=2;
v=2;
I=0;

% Calcul de l'évolution des populations

for t=0:dt:Tmax,

        I=I+1;
        u1=u;
        v1=v;
        u2=u1+dt*(alpha*u1-u1*v1)/2;
        v2=v1+dt*(-gamma*v1+u1*v1)/2;
        u3=u1+dt*(alpha*u2-u2*v2)/2;
        v3=v1+dt*(-gamma*v2+u2*v2)/2;
        u4=u1+dt*(alpha*u3-u3*v3);
        v4=v1+dt*(-gamma*v3+u3*v3);

        u=u+dt*((alpha*u1-u1*v1)/6 + (alpha*u2-u2*v2)/3 + (alpha*u3-u3*v3)/3 + (alpha*u4-u4*v4)/6);
        v=v+dt*((-gamma*v1+u1*v1)/6 + (-gamma*v2+u2*v2)/3 + (-gamma*v3+u3*v3)/3 + (-gamma*v4+u4*v4)/6);

if rem(I,10)==0
        U(I/10) = t + i*u;
        V(I/10) = t + i*v;
        W(I/10) = u + i*v;
end

end

% Affichage 

plot(U);
hold on
plot(V,'green');
title('Modele de Lotka Volterra par methode d Euler');
ylabel('u et v');
xlabel('temps');

% Affichage de la trajectoire dans le plan de phase

pause
clf;
plot(W);
title('Espace des phases : u v');
xlabel('u');
ylabel('v');
