Programation d'un system d´équa diff pour lsoda (package ode

[Johan Ransayer]

Bonjour

Je travail sur un model de dynamique des populations. J'ai besoin d'ecrire une fonction qui m'ecrive et resolve un systeme d'equations differentielles, a partir d'une matrice d'interaction decrivant un reseau trophique.

J'ai reussi a ecrire cette fonction pour un nombre connus d'especes. Voici mon script pour un reseau de 5 especes:

Code : Tout sélectionner

##----------------------------------------------------------------##
##    Write and solve the ODE model for a 5 species community     ##
##----------------------------------------------------------------##

Solve.fdwb.5sp <- function(y, pop.init = c(n1=10,n2=10,n3=10,n4=10,n5=10), step = 1000)
{
## The argument y must be a product of the function Draw.int.strength()

library(odesolve)

  # Creation of the ODE system:

qim <- y[[4]]                         # Interaction parameters

a11=qim[1,1]
a12=qim[1,2]
a13=qim[1,3]
a14=qim[1,4]
a15=qim[1,5]
a21=qim[2,1]
a22=qim[2,2]
a23=qim[2,3]
a24=qim[2,4]
a25=qim[2,5]
a31=qim[3,1]
a32=qim[3,2]
a33=qim[3,3]
a34=qim[3,4]
a35=qim[3,5]
a41=qim[4,1]
a42=qim[4,2]
a43=qim[4,3]
a44=qim[4,4]
a45=qim[4,5]
a51=qim[5,1]
a52=qim[5,2]
a53=qim[5,3]
a54=qim[5,4]
a55=qim[5,5]


k <- y[[3]]                          # Carrying capacity

K1=k[1]
K2=k[2]
K3=k[3]
K4=k[4]
K5=k[5]

  LV.model <- function(t, x, parms)
  {
  n1 <- x[1]
  n2 <- x[2]
  n3 <- x[3]
  n4 <- x[4]
  n5 <- x[5]

  # Populations dynamic ODEs

  dn1 <- (1-(n1/K1)) * (a11*n1) + (a12*n2*n1) + (a13*n3*n1) + (a14*n4*n1) + (a15*n5*n1)           
  dn2 <- (1-(n2/K2)) * ((a21*n1*n2) + (a22*n2)) + (a23*n3*n2) + (a24*n4*n2) + (a25*n5*n2)
  dn3 <- (1-(n3/K3)) * ((a31*n1*n3) + (a32*n2*n3) + (a33*n3)) + (a34*n4*n3) + (a35*n5*n3)
  dn4 <- (1-(n4/K4)) * ((a41*n1*n4) + (a42*n2*n4) + (a43*n3*n4) + (a44*n4)) + (a45*n5*n4)
  dn5 <- (1-(n5/K5)) * ((a51*n1*n5) + (a52*n2*n5) + (a53*n3*n5) + (a54*n4*n5)) + (a55*n5)

  return(list(c(dn1, dn2, dn3, dn4, dn5)))
  }

t <- seq(from=1, to=step, length=step)                               # Time sequence

parms <- c(a11,a12,a13,a14,a15,
           a21,a22,a23,a24,a25,
           a31,a32,a33,a34,a35,
           a41,a42,a43,a44,a45,
           a51,a52,a53,a54,a55,
           K1,K2,K3,K4,K5)

pop.dyn <- lsoda(y=pop.init, times=t, func=LV.model, parms)           #Solver

return(pop.dyn)
}


Voici un exemple des objets y[[4]] et y[[3]] utilise' dans cette fonction:

Code : Tout sélectionner


y[[4]]  # Matrice d'interaction entre 5 sp

          [,1]       [,2]       [,3]       [,4]       [,5]
[1,] 0.3085765 0.00000000 -0.3634194 -0.5988247 -0.4239148
[2,] 0.0000000 0.20881228  0.0000000 -0.2723391  0.0000000
[3,] 0.1205536 0.00000000  0.0000000 -0.3395979  0.0000000
[4,] 0.8772238 0.01030920  0.7149582  0.0000000 -0.2277903
[5,] 0.3504805 0.00000000  0.0000000  0.2576328 -0.1433585

y[[3]]  # Vecteur contenant la capacite' biotique de chaque sp

[1] 100.000000 100.000000  10.000000   4.641589   2.154435


Mon probleme est que je voudrais generaliser cette fonction pour pouvoir l'appliquer a des communautes d'especes de differentes tailles.

Si j'ai bien compris l'enjeu du probleme, il me faut ecrire une fonction qui sache elle meme ecrire un script du model dans un fichier texte (avec le bon nombre de parametres, d'equations, ...) puis qui rappele ce script pour en resoudre le systeme d'equation.

Je voudrais donc savoir (1) si ce raisonement est bon ?
(2) Si oui, est ce que quelqu'un saurait me donner au moins une piste sur la facon de programmer une telle chose.
(2 bis) Si non, comment faire ?

J'espere que quelqu'un saura m'aider
Merci d'avance
Johan

[Eric Pagot]

J'ai refais une petite fonction qui permet de s'affranchir du nombre de populations.
n est le vecteur de longueur variable (il est possible de reprendre 5 en mettant c(10,10,10,10,10))

Solve.fdwb <- function(y, pop.init = n, step = 1000)
{
## The argument y must be a product of the function Draw.int.strength()

library(odesolve)

# Creation of the ODE system:

qim <- y[[4]] # Interaction parameters
k <- y[[3]] # Carrying capacity

LV.model <- function(t, x, parms)
{
# Populations dynamic ODEs
n<-x
dyn<-NULL
for (i in 1:length(k)){
dn<-paste("dn",i,sep="")
assign(dn,(1-n[i]/k[i])*sum(qim[i,]*n[i]))
dyn<-c(dyn,eval(parse(text=dn)))
}
return(list(dyn)))
}

t <- seq(from=1, to=step, length=step) # Time sequence

parms <- c(qim,k)

pop.dyn <- lsoda(y=pop.init, times=t, func=LV.model, parms) #Solver

return(pop.dyn)
}

[Johan Ransayer]

Merci beaucoup Eric, cette fonction debloque pas mal de mes soucis ( je vais quand meme devoir la retoucher un peu, notament parce que ma mise en equation est un peu plus complique' que ca)

J'avais essaye' quelque chose dans ce gout, mais qui ne marchait pas :

Code : Tout sélectionner


...

 LV.model<-function(pop.init,t,parms)
  {
    for(i in 1:n.sp)
    {
    assign(paste("n",i,sep=""),pop.init[i])
    }

    for(i in 1:n.sp){
      for(j in 1:n.sp)
      {
      assign(paste("dn",i,sep=""),sum(qim[i,j]*(paste("n",j,sep=""))))
      }
    }
  ...
  return(...)
  }

...


Donc je pensais que cette methode etait a abandonner.
Merci de m'avoir remis sur les rails aussi rapidement et efficacement ! Je m'apercoit maintenant que mon raisonement n'etait peut-etre pas faux, mais beaucoup trop complique' !

Johan Ramsayer

[Johan Ransayer]

Bonjour,

Je remonte ce sujet car je n'arrive toujours pas à obtenir le même résultat avec la fonction généralisée (Solve.fdwb(), voir ci dessous) qu'avec la fonction pour 5 éspèces (Solve.fdwb.5sp()).

J'ai beau retourner le problème dans tous les sens, je ne comprends pas d'où viens la difference.

Il peut s'agir d'une erreur dans la version "calcul matriciel" des équations, mais j'ai verifié 10 fois, à la main, elles semblent bien donner le même résultat que les equations explicites.
Donc il doit s'agir d'une erreur de programation, mais je n'arrive pas à la trouver.


Si vous voulez essayer, voici de quoi utiliser les deux fonctions:


1) Je vous remet ma fonction d'origine, pour 5 éspèces seulement (y'a une petite erreur dans celle du premier message)

Code : Tout sélectionner

##---------------------------------------------------------------##
##    Write and solve the ODE model for a 5 species community    ##
##---------------------------------------------------------------##

Solve.fdwb.5sp<-function(x,pop.init=c(n1=10,n2=10,n3=10,n4=10,n5=10),step=1000)
{
## The argument x must be a product of the function Draw.int.strength()

library(odesolve)

  # Creation of the Ordinary Differential Equations model :

qim<-x[[4]]                         # Interaction parameters

a11=qim[1,1]
a12=qim[1,2]
a13=qim[1,3]
a14=qim[1,4]
a15=qim[1,5]
a21=qim[2,1]
a22=qim[2,2]
a23=qim[2,3]
a24=qim[2,4]
a25=qim[2,5]
a31=qim[3,1]
a32=qim[3,2]
a33=qim[3,3]
a34=qim[3,4]
a35=qim[3,5]
a41=qim[4,1]
a42=qim[4,2]
a43=qim[4,3]
a44=qim[4,4]
a45=qim[4,5]
a51=qim[5,1]
a52=qim[5,2]
a53=qim[5,3]
a54=qim[5,4]
a55=qim[5,5]


k<-x[[3]]                          # Carrying capacity

K1=k[1]
K2=k[2]
K3=k[3]
K4=k[4]
K5=k[5]

  LV.model<-function(t,x,parms)
  {
  n1<-x[1]
  n2<-x[2]
  n3<-x[3]
  n4<-x[4]
  n5<-x[5]

  dn1<-(1-(n1/K1)) * (a11*n1) + (a12*n2*n1) + (a13*n3*n1) + (a14*n4*n1) + (a15*n5*n1)            # Populations dynamic ODEs
  dn2<-(1-(n2/K2)) * ((a21*n1*n2) + (a22*n2)) + (a23*n3*n2) + (a24*n4*n2) + (a25*n5*n2)
  dn3<-(1-(n3/K3)) * ((a31*n1*n3) + (a32*n2*n3)) + (a33*n3) + (a34*n4*n3) + (a35*n5*n3)
  dn4<-(1-(n4/K4)) * ((a41*n1*n4) + (a42*n2*n4) + (a43*n3*n4)) + (a44*n4) + (a45*n5*n4)
  dn5<-(1-(n5/K5)) * ((a51*n1*n5) + (a52*n2*n5) + (a53*n3*n5) + (a54*n4*n5)) + (a55*n5)

  return(list(c(dn1,dn2,dn3,dn4,dn5)))
  }

t<-seq(from=1,to=step,length=step)                               # Time sequence

parms<-c(a11,a12,a13,a14,a15,
         a21,a22,a23,a24,a25,
         a31,a32,a33,a34,a35,
         a41,a42,a43,a44,a45,
         a51,a52,a53,a54,a55,
         K1,K2,K3,K4,K5)

pop.dyn<-lsoda(y=pop.init,times=t,func=LV.model,parms)           #Solver

return(pop.dyn)
}


2) Voici la version généralisée réécrite grâce à l'aide d'Eric

Code : Tout sélectionner

##---------------------------------------------------------------##
##    Write and solve the ODE model for a N species community    ##
##---------------------------------------------------------------##

Solve.fdwb<-function(y,pop.init=c(),step=1000)
{
## The argument y must be a product of the function Draw.int.strength()

library(odesolve)

qim<-y[[4]]                                                # Interaction parameters (Quantitative Interaction Matrix)
k<-y[[3]]                                                  # Carrying capacities for each species
n.sp<-length(k)

  LV.model<-function(t,x,parms)
  {
  n<-x
  dyn<-NULL
    for (i in 1:n.sp)
    {
    dn<-paste("dn",i,sep="")

      for (j in 1:n.sp)
      {
        if(qim[i,i]>0)                                     # Differential equation for the Primary-producer species (positive on the diagonal)
        {
        assign(dn, (((1-(n[i]/k[i]))* (diag(qim))[i]*n[i] )+ sum( (upper.tri(qim)*qim)[i,j] *n[i] *n[j] )) )
        }

        else                                               # Differential equation for the Non-primary-producer species (negative on the diagonal)
        {
        assign(dn,  ((1-(n[i]/k[i])) * sum( (lower.tri(qim)*qim)[i,j] *n[i] *n[j] ) + (sum( (upper.tri(qim)*qim)[i,j] *n[i] *n[j] ) +(diag(qim))[i]*n[i])) )
        }

      }
    dyn<-c(dyn,eval(parse(text=dn)))
    }
  return(list(dyn))
  }

t<-seq(from=1,to=step,length=step)                          # Time sequence
parms<-c(qim,k)                                             # Parameters

pop.dyn<-lsoda(y=pop.init,times=t,func=LV.model,parms)      # Solver

return(pop.dyn)
}


3) Voici la matrice d'une communauté pour faire tourner les deux fonctions

Code : Tout sélectionner

# Description de la communauté d'éspèces (y[[3]]:K capacités biotique théorique, y[[4]]: matrice d'interaction, y[[1]] et y[[2]] inutile ici.)

y<-list()

y[[3]]<-c(50.000000,50.000000,5.000000,5.000000,2.320794)

y[[4]]<-matrix(c(0.83712919,0.00000000,0.37591823,0.26398157,0.45241972,0.00000000,0.99060780,0.00000000,0.63598754,0.51360748,-0.49823316,0.00000000,-0.47475615,0.00000000,0.00000000,-0.21853796,-0.37361636,0.00000000,-1.58329479,0.72355146,-0.01944298,-0.51651342,0.00000000,-0.58447542,-1.79052702),5,5)

y


4) Resultat avec Solve.fdwb.5sp() :

Code : Tout sélectionner

pd<-Solve.fdwb.5sp(y)

  plot(x=pd[,"time"],y=pd[,"n1"],type="l",col="green",xlab="Time",ylab="Population size",ylim=c(0,max(pd[1000,-1])))
  lines(pd[,"time"],pd[,"n1"],col="green")
  lines(pd[,"time"],pd[,"n2"],col="blue")
  lines(pd[,"time"],pd[,"n3"],col="purple")
  lines(pd[,"time"],pd[,"n4"],col="orange")
  lines(pd[,"time"],pd[,"n5"],col="red")


5) Resultat avec Solve.fdwb() :

Code : Tout sélectionner

pd2<-Solve.fdwb(y,pop.init=c(n1=10,n2=10,n3=10,n4=10,n5=10))

plot(x=pd2[,"time"],y=pd2[,2],type="l",col="green",xlab="Time",ylab="Population size",ylim=c(0,50),xlim=c(0,1000))
  lines(pd2[,"time"],pd2[,2],col="green")
  lines(pd2[,"time"],pd2[,3],col="blue")
  lines(pd2[,"time"],pd2[,4],col="purple")
  lines(pd2[,"time"],pd2[,5],col="orange")
  lines(pd2[,"time"],pd2[,6],col="red")


Si vous essayez, vous verrez que l'on obtient:

*Avec la premiere fonction, un graphe où les 5 sp coexistent et sont en équilibre apres une période d'oscillations.

*Avec la fonction généralisée, un graphe où les deux éspèces Végetales (ici producteurs primaires = taux de croissance intrinsèque positif sur la diagonale et ne mangent aucune autre sp)
atteignent leur K théorique très rapidement, et toutes les autres éspèces disparaissent.

Dans le second cas, avec la même communauté, tout ce passe comme si la consomation de végétaux par les prédateurs ne permetait pas aux predateurs d'augmenter leur pop.
Car il y a bien consomation des végétaux : Si l'on zoom sur le début des courbes,

Code : Tout sélectionner

plot(x=pd2[,"time"],y=pd2[,2],type="l",col="green",xlab="Time",ylab="Population size",ylim=c(0,20),xlim=c(0,100))
     lines(pd2[,"time"],pd2[,2],col="green")
     lines(pd2[,"time"],pd2[,3],col="blue")
     lines(pd2[,"time"],pd2[,4],col="purple")
     lines(pd2[,"time"],pd2[,5],col="orange")
     lines(pd2[,"time"],pd2[,6],col="red")


on voit que l'éspèce bleue a commencé a décroitre avant de remonter, ce qui signifie qu'elle a été consommée.


Voila, je ne sais pas si tout ces détails sont utiles ou pas, mais si quelqu'un a une idée sur la raison pour laquelle ces deux fonctions ne donnent pas les mêmes résultats, ca m'aiderais vraiment beaucoup.
J'ai vraiment besoin que cette fonction marche et je n'ai plus aucune idée.

Merci d'avance et désolé pour la longueur du message.
Johan Ramsayer

[Timothee Poisot]

Pourquoi ne pas utiliser le package simecol?

Il est typiquement conçu pour ce genre de problèmes, et il gère très bien les populations de taille différentes. Je l'utilise pour ça tous les jours (simul. avec des spéciations/extinctions).

Je peux t'envoyer un script pour que tu voies à quoi ça ressemble.

[Johan Ransayer]

Ah oui, ca peut m'interesser aussi. Ca serait vraiment gentil si tu pouvais m'envoyer un script, j'aimerais bien voir à quoi ca ressemble.
Merci
Johan