Euler implicite et semi-implicite sur R

mharambo · August 2018

Bonjour à tous !
Je voudrais utiliser un schéma d'[large]E[/large]uler semi-implicite ou implicite pour mon modèle, seulement, j'ai beau lire les équations mathématiques, je ne comprends pas du tout comment l'on code ça.

Je n'ai pas envie d'utiliser un package (ex: deSolve) car je voudrais par la suite programmer mon modèle sous un autre logiciel où il n'y aura pas de package pour m'aider. J'ai donc besoin de comprendre comment l'on fait.

J'ai mis un code en exemple avec un shéma d'Euler explicite (beaucoup plus simple à comprendre). Mais vous pouvez me présenter un autre code si vous le souhaitez.
En vous remerciant,
Bonne journée

[En toute occasion, Leonhard Euler (1707-1783) prend une majuscule. AD]

#Donnees
#-------
parameters<-c(borne=0.4,fA1=0.08,mA2=0.02,KA2=1000,fA2=0.06,mA3=0.01)#parametres
state.vars <- c(A1=1000, A2=0, A3=0)#variable d'etat
dt<-0.01#pas de temps
ndays<-1000#temps

#fonction du modele: 
#-------------------
modele <- function(state.vars, parameters){
  
 invisible(mapply(assign, names(state.vars), state.vars, MoreArgs = list(envir = environment())))
 invisible(mapply(assign, names(parameters), parameters, MoreArgs = list(envir = environment())))
   
  #Equation du modele
   dA1dt = borne*A3  - A1* fA1
   dA2dt = fA1*A1 - A2*(mA2*(1.0+A2/KA2)+ fA2)
   dA3dt = fA2*A2 - A3*mA3

   #Euler explicite:
###########################
   A1 <- A1 + dA1dt * dt ##
   A2 <- A2 + dA2dt * dt ##
   A3 <- A3 + dA3dt * dt ##
###########################
   
 c("A1"=A1, "A2"=A2, "A3"=A3)
}


#Simulation: 
#-----------
data.Euler.exp <- NULL
collection <- data.frame(t(state.vars))

for (i in 1:ndays) { 
  state.vars <- modele(state.vars, parameters)
  collection <- rbind(collection, state.vars)
  
} #fin de la boucle temps

data.Euler.exp <- rbind(data.Euler.exp, collection)

#graph:
#-----
plot(data.Euler.exp$A1)  
plot(data.Euler.exp$A2)  
plot(data.Euler.exp$A3)

Saturne · August 2018

Bonjour,

Quels points te posent problème ?

mharambo · August 2018

Bonjour,
et bien je pense que mon problème vient du fait que je ne comprends pas le schéma.
Comment avoir y_n+1=y_n+hf_{(tn+1, yn+1)} alors que je ne connais pas f_{(tn+1, yn+1)} ??

Rescassol · August 2018

Bonjour,

$y_{n+1}=y_n+hf(t_{n+1}, y_{n+1})$ est une équation en $y_{n+1}$ à résoudre par le moyen que tu voudras (Newton, fsolve ...)

Cordialement,

Rescassol

Saturne · August 2018

C'est expliqué sur Wikiipedia.

Une méthode simple pour résoudre cette équation, c'est de partir de $y_{n+1}^{[0]} = y_n$, puis de calculer $y_{n+1}^{[1]} = y_n + h \times f(t_{n+1}, y_{n+1}^{[0]})$, puis $y_{n+1}^{[2]} = y_n + h \times f(t_{n+1}, y_{n+1}^{[1]})$, etc. jusqu'à ce que ce procédé converge.

Avec R, on peut faire comme ça par exemple :

f <- function(t, y){
  t * y
}

t <- 1 # t_{n+1}
h <- 0.1
y_n <- 2
# calcul de y_{n+1}
y0 <- y_n
y1 <- y_n + h * f(t, y0)
while(abs(y1-y0) > 1e-6){ # tant que le procédé n'a pas convergé
  y0 <- y1
  y1 <- y_n + h * f(t, y0)
}
y_nPlus1 <- y1 # voilà y_{n+1}

Ou alors, prendre un critère de convergence sur la différence relative entre deux termes successifs:

while(abs(y1-y0) > 1e-6 * abs(y0)){

J'ai pris $1\textrm{e-6}$ ici mais on peut prendre plus petit pour être plus précis.

Comprends-tu ? Est-ce que le reste pose un autre problème ?

Rescassol · August 2018

Bonjour,

Par exemple, en Python:

def Euler_implicite(f,xmin,xmax,xpas,y0):
    
    y, taby = y0, [y0]
    
    for x in np.arange(xmin,xmax-xpas,xpas):
        F = lambda z: y+xpas*f(z,x+xpas)-z
        y = fsolve(F,y)
        taby.append(y[0])
        
    return taby

Cordialement,

Rescassol

mharambo · August 2018

Bonjour,

Déjà merci beaucoup pour vos réponses, ça m'a aidé à mieux comprendre.
J'étais totalement passé à côté des exposants [0], [1], [2] de y.

Du coup j'ai modifié le script de départ avec le schéma d'Euler implicite (en espérant avoir bien compris). Il y a probablement de meilleures façons de le faire mais peut être que ça pourra en aider d'autres.

Merci encore !

Bonne journée ;-)

#Donnees
#-------
parameters<-c(borne=0.4,fA1=0.08,mA2=0.02,KA2=1000,fA2=0.06,mA3=0.01)#parametres
state.vars <- c(A1=1000, A2=0, A3=0)#variable d'etat
dt<-0.01#pas de temps
ndays<-1000#temps

#fonction du modele: 
#-------------------
modele <- function(state.vars, parameters){
  
  invisible(mapply(assign, names(state.vars), state.vars, MoreArgs = list(envir = environment())))
  invisible(mapply(assign, names(parameters), parameters, MoreArgs = list(envir = environment())))
  
  M=c(A1=0, A2=0, A3=0)
  
  #Equation du modele
  M["A1"] = borne*A3  - A1* fA1
  M["A2"] = fA1*A1 - A2*(mA2*(1.0+A2/KA2)+ fA2)
  M["A3"] = fA2*A2 - A3*mA3
  
 return(M)
 
  
}

#Euler implicite:
#-------------------

Euler_imp<-function(state.vars, parameters){
  
  invisible(mapply(assign, names(parameters), parameters, MoreArgs = list(envir = environment())))
  
  y0 <- state.vars #y0_n+1=y_n
  y1 <- y0 + dt * modele(y0, parameters) #y^1_n+1=y_n+h*f(t_n+1,y0_n+1)
  while(sum(abs(y1-y0) > 1e-6, na.rm=TRUE)>=1){
    y0 <- y1
    y1 <- state.vars + dt * modele(y0, parameters)#y^2(n+1)=y(n)+h*f(t_n+1,y1_n+1)
  }
return(y1)
}


#Simulation: 
#-----------
data.Euler.imp <- NULL
collection <- data.frame(t(state.vars))

for (i in 1:ndays) { 
  
  state.vars <-  Euler_imp(state.vars, parameters)
  collection <- rbind(collection, state.vars)
  

} #fin de la boucle temps

data.Euler.imp <- rbind(data.Euler.imp, collection)

#graph:
#-----
plot(data.Euler.imp$A1)  
plot(data.Euler.imp$A2)  
plot(data.Euler.imp$A3)

Euler implicite et semi-implicite sur R

Réponses

Bonjour!

Catégories

In this Discussion

Qui est en ligne 28