# ************************************************
#        SIMULATION du MODELE LOGISTIQUE         *
# ************************************************

# pour supprimer tous les objets du workspace courant (en memoire dans R)
rm(list=ls(all = TRUE)) 
# pour s'assurer de la reproductibilité des resultats (alea fixe)
set.seed(1)

# Packages a installer/charger
library(deSolve);

# definition de la fonction logistique qui permet de passer
# r et K en paramètres
logistique = function(x,r,K) {
  return(r*(1-x/K)*x)}

# paramètres du modèle logistique
r=0.5;K=10;x0=1;

# SOLUTION EXACTE
# ***************
# On calcule la solution exacte
tempsExact = seq(0,20,by=0.2); # temps auxquels on veut calculer la solution
solExact = rep(0,length(tempsExact)); # initialisation du vecteur solution
# Boucle sur les temps
for (i in 1:length(tempsExact)) {
  # calcul de la solution au ieme temps
  solExact[i]= x0*K/(x0-exp(-r*tempsExact[i])*(x0-K));}

# tracé de la solution
plot(tempsExact,solExact,col="red",main="Croissance Logistique - schéma d'Euler",
     xlab = "temps", ylab = "Effectif",type='l')
legend(x="topleft",c("Solution exacte"),col=c("red"),lty=c(1),bty="n")

# SOLUTION APPROCHEE
# ******************
# SCHEMA d'EULER EXPLICITE
# ************************
# pas de temps pour la simulation
deltaT = 2.0;

# temps de simulation
tempsDeSimulation=20;

# vecteur des instants auxquels la solution va etre calculée
vecteurTemps = seq(0,tempsDeSimulation,by=deltaT);

# initialisation du vecteur de stockage de la solution
effectif = rep(0,length(vecteurTemps));

# condition initiale
effectif[1]=x0;

# intégration numérique de l'équation via le schéma d'euler
for (i in 2:length(vecteurTemps)) {
  effectif[i] = effectif[i-1]+deltaT*logistique(effectif[i-1],r,K);
}

# tracé de la solution
lines(vecteurTemps,effectif,lty='dashed',col='cyan');
legend(x="topleft",c("Solution exacte","Schema Euler, dt=2"),col=c("red","cyan"),lty=c(1,2),bty="n")

# On refait la meme chose avec un pas de temps plus petit
deltaT = 0.5;
vecteurTemps = seq(0,20,by=deltaT);
effectif = rep(0,length(vecteurTemps));
effectif[1]=x0;
for (i in 2:length(vecteurTemps)) {
  effectif[i] = effectif[i-1]+deltaT*r*(1-effectif[i-1]/K)*effectif[i-1];
}
lines(vecteurTemps,effectif,lty='dashed',col='green');
legend(x="topleft",c("Solution exacte","Schema Euler, dt=2","Schema Euler, dt=0.5"),col=c("red","cyan","green"),lty=c(1,2,2),bty="n")

# et encore la même chose avec un pas encore plus petit
deltaT = 0.1;
vecteurTemps = seq(0,20,by=deltaT);
effectif = rep(0,length(vecteurTemps));
effectif[1]=x0;
for (i in 2:length(vecteurTemps)) {
  effectif[i] = effectif[i-1]+deltaT*r*(1-effectif[i-1]/K)*effectif[i-1];
}
lines(vecteurTemps,effectif,lty='dashed',col="blue");
legend(x="topleft",c("Solution exacte","Schema Euler, dt=2","Schema Euler, dt=0.5","Schema Euler, dt=0.1"),col=c("red","cyan","green","blue"),lty=c(1,2,2,2),bty="n")

# SCHEMA de RUNGE KUTTA d'ORDRE 4
# *******************************
# Comparaison avec la solution donnée par rk4, qui est un autre schéma numérique 
# Cela nécéssite que la fonction 'logistique' soit définie de la manière suivante
logistiqueRK = function(t,x,parms) {
  r = parms[1]; K = parms[2];
  dxsurdt = r*(1-x/K)*x;
  return(list(dxsurdt))
}

# vecteur des instants auxquels la solution va être calculée
tempsrk4=seq(0,20,by=2.0);

# simulation de l'équation différentielle à l'aide du schéma Runge Kutta d'ordre 4,
# codé dans la function rk4
sol=rk4(x0,tempsrk4,logistiqueRK,c(r,K));

dim(sol) # c'est une matrice: la première colonne est le vecteur temps, la seconde le vecteur des solutions

plot(sol[,1], sol[,2],main="Croissance logistique - schéma de Runge Kutta",col="black",lty='dashed',
  xlab="temps",ylab="effectif",type='l');

# Tracé pour comparaison avec la solution exacte
lines(tempsExact,solExact,col="red")
legend(x="topleft",c("Solution exacte","Schema Runge Kutta d'ordre 4"),col=c("red","black"),lty=c(1,1),bty="n")