# ************************************************ # 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")