# ************************************************************************ # PORTRAIT de PHASE du MODELE de COMPETITION de LOTKA VOLTERRA * # ************************************************************************ # 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) # Modèle de compétition de Lotka-Volterra # --------------------------------------- # dN1/dt=c1*N1*(1-(N1+alpha1*N2)/K1) # dN2/dt=c2*N2*(1-(N2+alpha2*N1)/K2) # # N1: population 1, N2: population 2 # # on fait le changement de variables: # x1=N1/K1, x2=N2/K2, r1=c1*K1, r2=c2*K2 # a1=alpha1*K2/K1, a2=alpha2*K1/K2 # et on tombe sur le système # # dx1/dt= r1*x1*(1-x1-a1*x2) # dx2/dt= r2*x2*(1-x2-a2*x1) library(deSolve); # fonction second membre du système différentiel correspondant au modèle de compétition de Lotka Volterra Competition <- function(t, x, parms) { with(parms, { dx1=r1*x[1]*(1-x[1]-a1*x[2]) dx2=r2*x[2]*(1-x[2]-a2*x[1]) return(list(c(dx1,dx2))) }) } # Calcul de la jacobienne de la fonction second membre Jacob <- function(x,parms) { with(parms, { J11 <- r1*(1-x[1]-a1*x[2])+r1*x[1]*(-1) J12 <- r1*x[1]*(-a1) J21 <- r2*x[2]*(-a2) J22 <- r2*(1-x[2]-a2*x[1])+r2*x[2]*(-1) J <- rbind(c(J11,J12), c(J21,J22)) return(J) }) } # paramètres du modèle parms <- list( r1=0.1, r2=0.2, # Tester les trois combinaisons de paramètres (a1,a2) ## Jeu de paramètres 1 #a1=0.7, #a2=0.8 # Jeu de paramètres 2 a1=1.2, a2=1.5 ## Jeu de paramètres 3 #a1=0.7, #a2=1.5 ) # points d'équilibre du modèle: calcul du point d'équilibre, de la jacobienne et de ses valeurs propres pour chacun des 4 points d'équilibre # E0 # valeur du point d'equilibre Xstar0=c(0,0); # appel de la fonction Jacob pour calculer la matrice jacobienne JXstar0 <- Jacob(Xstar0,parms) # calcul des valeurs propres eigen(JXstar0)$values # E1 Xstar1=c(1,0); JXstar1 <- Jacob(Xstar1,parms) eigen(JXstar1)$values # E2 Xstar2=c(0,1); JXstar2 <- Jacob(Xstar2,parms) eigen(JXstar2)$values # E3 # Xstar3 = point d'équilibre permettant la coexistence des 2 espèces M=rbind(c(1, parms[['a1']]), c(parms[['a2']],1)) Xstar3=solve(M)%*%c(1,1); JXstar3 <- Jacob(Xstar3,parms) eigen(JXstar3)$values # tracé du plan de phase espece 1 vs espece 2 plot(Xstar0[1],Xstar0[2],xlab="espece 1", ylab = "espece 2", type = "l", lwd = 2, xlim=c(0,1.2), ylim=c(0,1.2)) # affichage des points d'équilibre dans le plan de phase decal <- 0.05 points(Xstar0[1],Xstar0[2]) text(Xstar0[1]+decal,Xstar0[2]+decal,"E0") points(Xstar1[1],Xstar1[2]) text(Xstar1[1]+decal,Xstar1[2]+decal,"E1") points(Xstar2[1],Xstar2[2]) text(Xstar2[1]+decal,Xstar2[2]+decal,"E2") points(Xstar3[1],Xstar3[2]) text(Xstar3[1]+decal,Xstar3[2]+decal,"E3") # Affichage des valeurs propres de la jacobienne en les 4 points d'équilibre eigen(JXstar0)$values eigen(JXstar1)$values eigen(JXstar2)$values eigen(JXstar3)$values # Tracé du portrait de phase intéractif # -------------------------------------- # vecteur des instants auxquels la solution va etre calculee times <- seq(from = 0, to = 500, by = 1) while (1) { # boucle infinie print('Choisissez un nouvelle condition initiale en cliquant sur le plan de phase') x0 <- locator(1) # calcul de la solution partant de cette condition initiale out <- ode(as.numeric(x0), times, Competition, parms) # trace de la solution dans le plan de phase lines(out[,"1"], out[,"2"]) }