Euler implicite et semi-implicite sur R
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]
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)
Réponses
-
Bonjour,
Quels points te posent problème ? -
Bonjour,
et bien je pense que mon problème vient du fait que je ne comprends pas le schéma.
Comment avoir yn+1=yn+hf(tn+1, yn+1) alors que je ne connais pas f(tn+1, yn+1) ?? -
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 -
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 ? -
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 -
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)
Connectez-vous ou Inscrivez-vous pour répondre.
Bonjour!
Catégories
- 163.2K Toutes les catégories
- 9 Collège/Lycée
- 21.9K Algèbre
- 37.1K Analyse
- 6.2K Arithmétique
- 53 Catégories et structures
- 1K Combinatoire et Graphes
- 11 Sciences des données
- 5K Concours et Examens
- 11 CultureMath
- 47 Enseignement à distance
- 2.9K Fondements et Logique
- 10.3K Géométrie
- 64 Géométrie différentielle
- 1.1K Histoire des Mathématiques
- 68 Informatique théorique
- 3.8K LaTeX
- 39K Les-mathématiques
- 3.5K Livres, articles, revues, (...)
- 2.7K Logiciels pour les mathématiques
- 24 Mathématiques et finance
- 314 Mathématiques et Physique
- 4.9K Mathématiques et Société
- 3.3K Pédagogie, enseignement, orientation
- 10K Probabilités, théorie de la mesure
- 773 Shtam
- 4.2K Statistiques
- 3.7K Topologie
- 1.4K Vie du Forum et de ses membres