Resolution d'equation differentielle sous R

[fantomex]

bonjour, me voila heurter depuis des jours et des jours...
je dois resoudre le systeme d'equations suivantes à l'aide du logiciel R,
mais je n'y arrive pas.. je ne vois pas où est l'erreur!!!!!

les hypothèses de departs sont les suivantes:
- 6 espèces de phytoplancton (Ni ), 3 ressources ( Rj).
- dNi/dt = Ni (µi(R1,...,Rk)-mi), i=1,...,n (1)
- dRj/dt = D (Sj-Rj) - ∑ Cji µi (R1,...,Rk) Ni, j=1,...,k (2)
- µi(R1,...Rk) = min ( (riR1 / K1i + R1),...,(riRk / Kki+Rk))
-Valeurs des paramètres :
S1=6, S2=10, S3=14
r=1, D=0.25, m=D

(1 0,9 0,3 1,04 0,34 0,77 )
K=(0,3 1 0,9 0,71 1,02 0,76 )
(0,9 0,3 1 0,46 0,34 1,07 )

(0,04 0,07 0,04 0,1 0,03 0,02 )
C=(0,08 0,08 0,1 0,1 0,05 0,17 )
(0,14 0,1 0,1 0,16 0,06 0,14 )

Les matrices K et C comprennent en ligne les ressources (3 lignes = 3 ressources) et en colonnes les espèces (6 colonnes = 6 espèces).

- Valeurs initiales pour les variables :
R1=S1, R2=S2, R3=S3
Ni=0.1+i/100 pour les trois première espèces présentes à t=0.
N4=0.1
N5=0.1
N6=0.1

- Introduction des espèces dans le système :
Les 3 premières espèces sont présentes dès t0.
La 4ième espèce est introduite à t=1000, la 5ième à t=2000, la 6ième à t=5000.


Pour l'instant, je voudrai juste a reussir avec les trois especes initialements presente (sans se preoccuper du reste).

voici le script sur lequel je bloque :


S<-c(6,10,14)

r<-1 #taux de croissance spécifique maximum
D<-0.25 #taux de renouvellement du système
m<-D #taux de mortalité spécifique

#Création matrice Kij, constante de demie-saturation en ressource j pour l'espèce i
K<-matrix(NA,3,6)
K[1,]<-c(1.00,0.90,0.30,1.04,0.34,0.7 7)
K[2,]<-c(0.30,1.00,0.90,0.71,1.02,0.7 6)
K[3,]<-c(0.90,0.30,1.00,0.46,0.34,1.0 7)

#Création matrice C, teneur en ressource j dans l'espèce i
C<-matrix(NA,3,6)
C[1,]<-c(0.04,0.07,0.04,0.10,0.03,0.0 2)
C[2,]<-c(0.08,0.08,0.10,0.10,0.05,0.1 7)
C[3,]<-c(0.14,0.10,0.10,0.16,0.06,0.1 4)

#valeurs initiales des ressources
R<-S

#valeurs initiales des espèces
N<-c(0.11,0.12,0.13)

#Temps des 3 1ères espèces
Temps<-seq(0,15000,1000)


#equation différentielle

equ<-function(t,y,parms){
for(i in 1:3){
dN1_dt<-y[1]*(min(r*y[4:6]/(K[,1]+y[4:6]))-m)
dN2_dt<-y[2]*(min(r*y[4:6]/(K[,2]+y[4:6]))-m)
dN3_dt<-y[3]*(min(r*y[4:6]/(K[,3]+y[4:6]))-m)
dR1_dt<-D*(S[1]-y[1])-sum(C[1,]*min(r*y[4:6]/(K[,i]+y[4:6]))*y[1:3])
dR2_dt<-D*(S[2]-y[2])-sum(C[2,]*min(r*y[4:6]/(K[,i]+y[4:6]))*y[1:3])
dR3_dt<-D*(S[3]-y[3])-sum(C[3,]*min(r*y[4:6]/(K[,i]+y[4:6]))*y[1:3])
return<-list(c(dN1_dt,dN2_dt,dN3_dt,dR 1_dt,dR2_dt,dR3_dt))}}

#resolution equation differentielle
Res.equ<-lsoda(equ,y=c(N[1],N[2],N[3],R[1],R[2],R[3]),times=Temps,parms=NULL)
Res.equ
-----

[invitec5eb4b89]

Bon, j'ai essayé de débugger ton code, mais c'est vraiment fastidieux : essaye avec un problème plus simple !!
De plus, est-ce que ce n'est pas dangereux d'avoir "parms = NULL" ? Il va peut être falloir nommer chaque paramètre ("k11", "k12" etc.).

[invite986312212]

bonjour,

j'ai l'impression que ta fonction "equ" n'est pas correcte: tu y crées une variable "return" dans une boucle (donc ses valeurs successives sont écrasées), mais elle ne renvoie rien.

[fantomex]

Heu... merci de votre aide en tout cas, en fait, le graphique que j'aimerai obtenir devrait etre sensiblement ressemblant a celui là (je ne considere que 3 especes) :

[A voir en vidéo sur Futura]

[invitec5eb4b89]

Je seconde la remarque d'Ambrosio.
En simplifiant le problème (j'ai enlevé les min, gardé deux espèces et deux ressources) :

rm(list=ls())
require(odesolve)

r <- 1    #taux de croissance spécifique maximum
D <- 0.25 #taux de renouvellement du système
m <- D    #taux de mortalité spécifique

# Constantes kij de demie-saturation en ressource j pour l'espèce i
#            cij de teneur en ressource j pour l'espèce i
p <- c(k11 = 1.00, k12 = 0.90, k21 = 0.30, k22 = 1.00,
       c11 = 0.04, c12 = 0.07, c21 = 0.08, c22 = 0.08)  

# Valeurs initiales des ressources
R <- S <- c(6,10)

# Valeurs initiales des espèces
N <- c(0.11,0.12)

# Temps des 2 1ères espèces
Temps <- seq(0,100)

# Equation différentielle

equ <- function(t,y,parms){
  dN1 <- y[1]*( r*y[3]/( p["k11"]+y[3] )-m)
  dN2 <- y[2]*( r*y[4]/( p["k22"]+y[4] )-m)
  dR1 <- D*(S[1]-y[1])-p["c11"]*( r*y[3]/( p["k11"]+y[3] ))*y[1]
  dR2 <- D*(S[2]-y[2])-p["c22"]*( r*y[4]/( p["k22"]+y[4] ))*y[2]
  list(c(dN1,dN2,dR1,dR2))
}

# Resolution equation differentielle
r <-lsoda(y=c(N[1],N[2],R[1],R[2]),times=Temps,equ,parms=p)
r
matplot(r[,2:3],type='l')

Ça ne résout pas ton problème, mais au moins tu peux repartir d'un code qui apparemment fonctionne.

[fantomex]

Merci beaucoup, je devrai m'en sortir a present!
merci les mathematiciens!

[fantomex]

Bonsoir, je sollicite de nouveau votre aide avec mon dernier script...
Tout marche, sauf la fin, en effet, R me dit que j'ai 21 derivés a ma fonction! je ne comprend pas du tout! si j'enleve "dR3_dt" (par exemple) a la ligne "list", il me dira que j'ai 15 derivés. Il me compte 6 derivés par dR _dt....
Je ne comprend plus

Voici le script :



rm(list=ls())
library(odesolve)


r<-1 #taux de croissance spécifique maximum
D<-0.25 #taux de renouvellement du système
m<-D #taux de mortalité spécifique


#Création matrice Kij, constante de demie-saturation en ressource j pour l'espèce i
K<-matrix(NA,3,6)
K[1,]<-c(1.00,0.90,0.30,1.04,0.34,0.7 7)
K[2,]<-c(0.30,1.00,0.90,0.71,1.02,0.7 6)
K[3,]<-c(0.90,0.30,1.00,0.46,0.34,1.0 7)


#Création matrice C, teneur en ressource j dans l'espèce i
C<-matrix(NA,3,6)
C[1,]<-c(0.04,0.07,0.04,0.10,0.03,0.0 2)
C[2,]<-c(0.08,0.08,0.10,0.10,0.05,0.1 7)
C[3,]<-c(0.14,0.10,0.10,0.16,0.06,0.1 4)


#valeurs initiales des ressources
R<-S<-c(6,10,14)


#valeurs initiales des espèces
N1<-0.11
N2<-0.12
N3<-0.13

#Temps des 3 1ères espèces
Temps<-seq(0,15000,1000)


#equation différentielle

equ<-function(t,y,parms){
dN1_dt<-y[1]*((min((r*y[4])/(K[1,1]+y[4]),(r*y[5])/(K[2,1]+y[5]),(r*y[6])/(K[3,1]+y[6])))-m)
dN2_dt<-y[2]*((min((r*y[4])/(K[1,2]+y[4]),(r*y[5])/(K[2,2]+y[5]),(r*y[6])/(K[3,2]+y[6])))-m)
dN3_dt<-y[3]*((min((r*y[4])/(K[1,3]+y[4]),(r*y[5])/(K[2,3]+y[5]),(r*y[6])/(K[3,3]+y[6])))-m)
dR1_dt<-D*(S[1]-y[4])-(C[1,]*(min((r*y[4])/(K[1,1]+y[4])+(r*y[5])/(K[2,1]+y[5])+(r*y[6])/(K[3,1]+y[6])))*y[1])
dR2_dt<-D*(S[2]-y[5])-(C[2,]*(min((r*y[4])/(K[1,2]+y[4])+(r*y[5])/(K[2,2]+y[5])+(r*y[6])/(K[3,2]+y[6])))*y[2])
dR3_dt<-D*(S[3]-y[6])-(C[3,]*(min((r*y[4])/(K[1,3]+y[4])+(r*y[5])/(K[2,3]+y[5])+(r*y[6])/(K[3,3]+y[6])))*y[3])
list(c(dN1_dt,dN2_dt,dN3_dt,dR 1_dt,dR2_dt,dR3_dt))
}


#resolution equation differentielle
Res.equ<-lsoda(y=c(N1,N2,N3,R[1],R[2],R[3]),times=Temps,equ,parms=NULL)
Res.equ

[invitec5eb4b89]

Pour les dRi_dt, le minimum n'est pas au bon endroit, ce qui fait que ce sont des vecteurs.

[fantomex]

Merci de votre réponse, mais je ne comprend pas, pour les "dNi_dt", où il ya également cette fonction "min", il n'y a pas de probleme. Alors pourquoi m'en pose t-il pour les dRi_dt ???
Je n'ai toujours pas trouver l'erreur... comment faire ????
J'ai vraiment du mal avec "R"...

[invitec5eb4b89]

Le problème ne vient pas de la fonction min : comme C[i,] est un vecteur, le multiplier par un scalaire (même si c'est le minimum de quelque chose) rend un vecteur.

[invitec5eb4b89]

...
equ<-function(t,y,parms){
dN1_dt<-y[1]*((min((r*y[4])/(K[1,1]+y[4]),(r*y[5])/(K[2,1]+y[5]),(r*y[6])/(K[3,1]+y[6])))-m)
dN2_dt<-y[2]*((min((r*y[4])/(K[1,2]+y[4]),(r*y[5])/(K[2,2]+y[5]),(r*y[6])/(K[3,2]+y[6])))-m)
dN3_dt<-y[3]*((min((r*y[4])/(K[1,3]+y[4]),(r*y[5])/(K[2,3]+y[5]),(r*y[6])/(K[3,3]+y[6])))-m)
dR1_dt<-D*(S[1]-y[4])-(C[1,]*(min((r*y[4])/(K[1,1]+y[4])+(r*y[5])/(K[2,1]+y[5])+(r*y[6])/(K[3,1]+y[6])))*y[1])
dR2_dt<-D*(S[2]-y[5])-(C[2,]*(min((r*y[4])/(K[1,2]+y[4])+(r*y[5])/(K[2,2]+y[5])+(r*y[6])/(K[3,2]+y[6])))*y[2])
dR3_dt<-D*(S[3]-y[6])-(C[3,]*(min((r*y[4])/(K[1,3]+y[4])+(r*y[5])/(K[2,3]+y[5])+(r*y[6])/(K[3,3]+y[6])))*y[3])
list(c(dN1_dt,dN2_dt,dN3_dt,dR 1_dt,dR2_dt,dR3_dt))
}

Au temps pour moi, après relecture des équations différentielles (1er message), ce n'est pas le minimum qui est mal placé, mais la somme
qui manque :

equ<-function(t,y,parms){
dN1_dt<-y[1]*((min((r*y[4])/(K[1,1]+y[4]),(r*y[5])/(K[2,1]+y[5]),(r*y[6])/(K[3,1]+y[6])))-m)
dN2_dt<-y[2]*((min((r*y[4])/(K[1,2]+y[4]),(r*y[5])/(K[2,2]+y[5]),(r*y[6])/(K[3,2]+y[6])))-m)
dN3_dt<-y[3]*((min((r*y[4])/(K[1,3]+y[4]),(r*y[5])/(K[2,3]+y[5]),(r*y[6])/(K[3,3]+y[6])))-m)
dR1_dt<-D*(S[1]-y[4])-sum(C[1,]*(min((r*y[4])/(K[1,1]+y[4])+(r*y[5])/(K[2,1]+y[5])+(r*y[6])/(K[3,1]+y[6])))*y[1])
dR2_dt<-D*(S[2]-y[5])-sum(C[2,]*(min((r*y[4])/(K[1,2]+y[4])+(r*y[5])/(K[2,2]+y[5])+(r*y[6])/(K[3,2]+y[6])))*y[2])
dR3_dt<-D*(S[3]-y[6])-sum(C[3,]*(min((r*y[4])/(K[1,3]+y[4])+(r*y[5])/(K[2,3]+y[5])+(r*y[6])/(K[3,3]+y[6])))*y[3])
list(c(dN1_dt,dN2_dt,dN3_dt,dR1_dt,dR2_dt,dR3_dt))
}