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