Résolution d'un système d'edo avec python
Bonjour,
je souhaite résoudre le système ci dessous avec Python.
Ma difficulté est dans la présence de $I_{\tau}$.
S'il n y avait pas $I_{\tau}$ alors je sais faire directement avec Python,
mais la présence de $I_{\tau}$ me perturbe beaucoup.
Je vous remercie d'avance pour votre aide.
On considère le système suivant : pour tout $t \in [0,T]$ :
\begin{align*}
U'&=-a U V,
\\
I'&=aUV-\beta I,
\\
V'&= b I_{\tau} - \sigma V,
\end{align*} où
$$
I_{\tau}(t)= \begin{cases}0 \ &\mbox{si} \ t \leq \tau, \\ I(t-\tau)\ &\mbox{si} \ t > \tau,
\end{cases}
$$ avec les conditions initiales
$$
U(0)=u_0, \ I(0)=0, \ V(0)=v_0,
$$ $a,b,\beta,\sigma$ et $\tau$ sont donnés, par exemple
$u_0=1,v_0=10,a=10,b=5,\beta=1,\tau=10$.
je souhaite résoudre le système ci dessous avec Python.
Ma difficulté est dans la présence de $I_{\tau}$.
S'il n y avait pas $I_{\tau}$ alors je sais faire directement avec Python,
mais la présence de $I_{\tau}$ me perturbe beaucoup.
Je vous remercie d'avance pour votre aide.
On considère le système suivant : pour tout $t \in [0,T]$ :
\begin{align*}
U'&=-a U V,
\\
I'&=aUV-\beta I,
\\
V'&= b I_{\tau} - \sigma V,
\end{align*} où
$$
I_{\tau}(t)= \begin{cases}0 \ &\mbox{si} \ t \leq \tau, \\ I(t-\tau)\ &\mbox{si} \ t > \tau,
\end{cases}
$$ avec les conditions initiales
$$
U(0)=u_0, \ I(0)=0, \ V(0)=v_0,
$$ $a,b,\beta,\sigma$ et $\tau$ sont donnés, par exemple
$u_0=1,v_0=10,a=10,b=5,\beta=1,\tau=10$.
Connectez-vous ou Inscrivez-vous pour répondre.
Réponses
Merci
"Résoudre" me parait un peu trop espérer pour un tel système d'équations différentielles non linéaire et comportant en plus un "retard".
Tout au plus pourra-t-on obtenir une solution approchée.
La question est-elle : comment programmer en Python pour obtenir une solution approchée ?
Dans ce cas, veux-tu le faire "à la main" ou bien en utilisant une bibliothèque existante ?
Bonjour,
j'ai travaillé et j'ai fait le code ci-dessous.
Avec Trinket il fonctionne mais les graphiques ne sont pas jolis
Si j’exécute ça avec Python 3.10 sous windows alors ça donne l'erreur : 4NameError: name 'plt' is not defined
Comment faire pour afficher les jolis graphes ? Svp
import numpy as np from math import pi from math import sin from math import cos from math import log from math import sqrt import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t)=nb non inféctées au temps t # i(t)=nb infectes au temps t # v(t)= nb cellules du virus au temps t # a #b #beta #sigma #du/dt=-auv #i'=auv-beta i #v'=b i(t-tau)-sigma v #unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0 V=10. tau=10 a=10 b=5 beta=0.1 sigma=0.4 t=0. temps.append(t) #pas de temps h=0.01 ntau=int(tau/h) #temps final T=500. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1-h*a*V) In=I*(1-beta*h)+a*h*U*V Vn=V-h*sigma*V if n > ntau : Vn=Vn-h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print (U, I , V) plt.plot(temps,UU,'blue') plt.plot(temps,VV,'red') plt.plot(temps,VV,'green') plt.show()Avec Python 3.8.3 dans Pyzo, il n'y a pas d'erreur.
Cordialement,
Rescassol
Si tu exécutes le tout, ça marche sans problème.
Est-ce qu'il y a comment l'executer en ligner comme trinket mais là où on obtient de belles courbes? Svp
Que faire? Svp
Par ailleurs, il y a une faute de frappe à la fin : devrait probablement être :
Enfin, il y a également une erreur de signe dans le calcul de $V_n$ lorsque $n>n_{\tau}$. Ce devrait être :
if n > ntau : Vn=Vn+h*b*II[n-ntau]Avec ces changements, on obtient quelque chose... qui pourrait être correct si la valeur du pas avait été mieux choisie.
Mais là, à la première itération, $U(1)$ devient nul car $haV(0)=1$.
Je ne pense pas que ce soit le comportement attendu... mais diminuer le pas ne fait pas beaucoup mieux.
Le choix des autres paramètres n'est sans doute pas correct non plus.
Ca me donne un graphe pour t de 0 à 10 seulement or que le T final est 500.
Que faire pour avoir un graphique pour t au delà de 10?
Comment l'améliorer pour obtenir un meilleur résultat? Svp
Il y a un truc bizarre: l'initialisation de V est 0.45 mais sur le graphe ça commence par 0!
Et le $U$ aussi est négatif. Ce n'est pas bon
Voici le code corrigé.
import numpy as np from math import pi from math import sin from math import cos from math import log from math import sqrt import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t)=nb non inféctées au temps t # i(t)=nb infectes au temps t # v(t)= nb cellules du virus au temps t # a #b #beta #sigma #du/dt=-auv #i'=auv-beta i #v'=b i(t-tau)-sigma v #unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0 V=0.45 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.01 ntau=int(tau/h) #temps final T=500 #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1-h*a*V) In=I*(1-beta*h)+a*h*U*V Vn=V-h*sigma*V if n > ntau : Vn=Vn+h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print (U, I , V) #plt.plot(temps,UU,'blue') #plt.plot(temps,II,'red') plt.plot(temps,VV,'green') plt.show()si on veut le graphe de log V(t) alors qu’elle quelle commande utiliser ? Svp
Cordialement.
Peut-être que la méthode de résolution pourrait être un peu améliorée, genre utiliser Runge-Kutta RK4 plutôt que ce qui me semble être Euler explicite. Certes il faut compliquer un peu le code...
En Python, il y a la fonction quad et ses avatars.
Cordialement,
Rescassol
Je pense que le problème le plus crucial est la valeur des paramètres $\tau, a, b, \beta$ et $\sigma$ (en faisant en particulier attention aux unités !!) ainsi que les valeurs initiales de $U, I$ et $V$.
Comme rien n'est vraiment clair dans ce que tu nous as fourni, on ne sait pas quelles sont les unités de tes différentes grandeurs (en particulier, sont-elles comparables, pour être mises sur un même graphique ) et on ne sait pas à quoi correspondent tes différents paramètres ni sur quoi tu te bases pour les choisir.
Il n'y a pas de miracle : quand on remplit "au pif", on obtient "n'importe quoi"...
Pouvez-vous m'aider à l'améliorer avec RK4? Svp.
Pour les paramètres ils sont bons j'en suis sûre. Je ne cherche pas à mettre tous le monde sur le même graphique. Chaque solution dans un graphique à part
import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t) = nb non infectées au temps t # i(t) = nb infectées au temps t # v(t) = nb cellules du virus au temps t # a # b # beta # sigma # u' = -auv # i' = auv - beta i # v' = b i(t-tau)-sigma v # unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0. V=0.45 tau=2 a=0.01 b=800 # J'ai changé la valeur car elle donnait des résultats incohérents. beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.01 ntau=int(tau/h) #temps final T=100. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1-h*a*V) In=I*(1-beta*h)+a*h*U*V Vn=V*(1-h*sigma) if n > ntau : Vn=Vn+h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print(U, I, V) f, (ax1, ax2, ax3) = plt.subplots(3,1) # On fait 3 graphes en colonne ax1.plot(temps,UU,'blue') ax2.plot(temps,II,'red') ax3.plot(temps,VV,'green') f.show()Voici le code. Quand je rajoute la ligne pour afficher log V, ça marque erreur.
Comment afficher le graphe de log V? Svp
Comment améliorer ce code avec RK4? En sachant que les données des paramètres sont corrects
Je vous remercie d'avance
import numpy as np from math import pi from math import sin from math import cos from math import log from math import sqrt import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t)=nb non inféctées au temps t # i(t)=nb infectes au temps t # v(t)= nb cellules du virus au temps t # a #b #beta #sigma #du/dt=-auv #i'=auv-beta i #v'=b i(t-tau)-sigma v #unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0 V=0.45 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.001 ntau=int(tau/h) #temps final T=25 #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1-h*a*V) In=I*(1-beta*h)+a*h*U*V Vn=V-h*sigma*V if n > ntau : Vn=Vn+h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print (U, I , V) #plt.plot(temps,UU,'blue') #plt.plot(temps,II,'red') #plt.plot(temps,VV,'green') VL = map(log, VV) plt.plot(temps, VL, 'green') plt.show()L'utilisation d'un schéma d'Euler implicite corrige aisément ce type d'erreur.
En revanche, je ne sais pas appliquer un schéma de Runge-Kutta avec retard.
Quelles modifications faire dans le cas de Euler implicite? Svp. Est-ce qu'il suffit de changer les signes des équations?
On met I-I_n au lieu de I_n-I?
Sinon pour le plot de log(V) ça ne marche pas. Que faire? Svp
Merci beaucoup d'avance
Comment faire le plot de log(V) ? Svp. Ça m'aidera beaucoup à voir si c'est bon.
import numpy as np from math import pi from math import sin from math import cos from math import log from math import sqrt import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t)=nb non inféctées au temps t # i(t)=nb infectes au temps t # v(t)= nb cellules du virus au temps t # a #b #beta #sigma #du/dt=-auv #i'=auv-beta i #v'=b i(t-tau)-sigma v #unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0 V=0.45 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.001 ntau=int(tau/h) #temps final T=25 #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1+h*a*V) In=I*(1+beta*h)-a*h*U*V Vn=V+h*sigma*V if n > ntau : Vn=Vn-h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print (U, I , V) plt.plot(temps,UU,'blue') #plt.plot(temps,II,'red') #plt.plot(temps,VV,'green') plt.show()Voilà un exemple de fonction implémentant Euler implicite sur une équation de la forme $y'=f(y,x)$.
D'autre part, c'est "odeint" que je voulais dire dans mon précédent message et non "quad".
def f(y,x): return np.sin(x**2+y) 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 def ExoI_2(): xmin, xmax, xpas, y0 = 0, 10, 0.01, 0 tabx = np.arange(xmin,xmax,xpas) taby = Euler_implicite(f,xmin,xmax,xpas,y0) plt.figure(linewidth=10,facecolor = 'gold',edgecolor='red') plt.plot(tabx,taby,'r',linewidth=3) plt.title("Méthode d'Euler implicite",fontsize=40,fontweight='bold',color='m') plt.show()Cordialement,Rescassol
le souci est que je n'arrive pas à faire le plot jusqu'à t=25
Merci de m'aider à améliorer le code. Ci joint aussi le plot de log V qu'il faut obtenir (en bleu)
Je vous remercie d'avance pour votre aide.
import numpy as np from math import pi from math import sin from math import cos from math import log from math import sqrt import matplotlib.pyplot as plt #modele sir # tau = le délai # u(t)=nb non inféctées au temps t # i(t)=nb infectes au temps t # v(t)= nb cellules du virus au temps t # a #b #beta #sigma #du/dt=-auv #i'=auv-beta i #v'=b i(t-tau)-sigma v #unite=heures UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0 V=0.45 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.001 ntau=int(tau/h) #temps final T=10000 #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) dim=0. for n in range(M) : Un=U*(1+h*a*V) In=I*(1+beta*h)-a*h*U*V Vn=V+h*sigma*V if n > ntau : Vn=Vn-h*b*II[n-ntau] U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) print (U, I , V) #plt.plot(temps,UU,'blue') #plt.plot(temps,II,'red') plt.plot(temps,VV,'green') plt.gca().set_yscale("log") plt.show()J'ai obtenu les équations ci-dessous, et les graphes correspondants.
import matplotlib.pyplot as plt from numpy import log #L'avantage de la fonction log de numpy, c'est qu'elle est déjà vectorisée... UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0. V=0.45 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.01 ntau=int(tau/h) #temps final T=25. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) for n in range(M) : Vn = (V + (0 if n <= ntau else h*b*II[n-ntau]))/(1+h*sigma) Un = U/(1 + a*h*Vn) In = (I + a*h*Un*Vn)/(1+beta*h) U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) t=t+h temps.append(t) if U < 0 or I < 0 or V < 0: #Si l'une des valeurs devient négative, c'est qu'il y a un bug !! break print(U, I, V) # Affichage de 4 graphes superposés, avec la même échelle de temps. f, (ax1, ax2, ax3, ax4) = plt.subplots(4,1) ax1.plot(temps,UU,'blue') ax2.plot(temps,II,'red') ax3.plot(temps,VV,'green') ax4.plot(temps,log(VV),'green') f.show()Le "retard" devient encore plus flagrant si on prend une valeur plus élevée pour $\tau$.
Svp pour les graphes, je souhaite afficher uniquement la courbe de log(V) et pas les 4 graphes. Du coup il faut enlever la commande f, (ax1, ax2, ax3, ax4) = plt.subplots(4,1), par quoi la remplacer? Svp pour avoir uniquement un seul graphe celui de log(V)
Merci pour toute votre aide
on remarque que les valeurs de ce graphique sont deux fois plus grandes que le graphique qu'on devait obtenir (le deuxième graphe joint). Comment l'expliquer? Svp
Quelles modifications faire ? Svp
Tu remplaces ta ligne
par :
puis tu remplaces log(VV) par log10(VV) dans ton code.
Alain
Merci beaucoup à tous!:-)
est-ce que c'est possible avec Python d'enregistrer les valeurs de log10(V) par exemple, dans un fichier .txt?
Si oui, comment? Svp
Merci d'avance
J'utilise la commande suivante: print(t,log10(V)), pour afficher les valeurs de t et les valeurs de $\log10(V)$ qui y correspondent. Ce que je souhaite, c'est de sauvegarder ces données sur un fichier .txt. Je ne trouve pas ce genre d'exemple dans la documentation. Si tu sais comment faire, merci de m'aider
Tu exagères, il suffit de taper "Python fichier" dans google et tu as ce que tu cherches.
Cordialement,
Rescassol
en haut de ton fichier tu dois importer le module json donc tu ajoutes un
puis là où tu veux écrire dans ton fichier tu fais un
with open('nom_de_fichier.json', 'w') as f: json.dump(log10(V).tolist(), f)Quelques commentaires sur le code ci-dessus pour mieux comprendre :
1) Il faut savoir que vu que tu veux sauvegarder une liste, le format approprié de ton fichier est le json et pas le txt. Donc ton fichier dois se nommer "nom_de_fichier.json" comme indiqué ci-dessus.
2) Du dois appeler la méthode tolist sur ton objet log10 pour le transformer en liste car log10 est une fonction de numpy qui renvoie un objet de type ndarray et la méthode json.dump attend une liste plutôt.
3) Bref si tu n'as jamais fait de Python les points 1) et 2) risquent d'être du chinois :-D
import matplotlib.pyplot as plt from numpy import log10 #L'avantage de la fonction log de numpy, c'est qu'elle est déjà vectorisée... import json UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0. V=4.5 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.0001 ntau=int(tau/h) #temps final T=25. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) for n in range(M) : Vn = (V + (0 if n <= ntau else h*b*II[n-ntau]))/(1+h*sigma) Un = U/(1 + a*h*Vn) In = (I + a*h*Un*Vn)/(1+beta*h) U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) #print(t,log10(V)) with open('log.json', 'w') as f: json.dump(log10(V).tolist(), f) t=t+h temps.append(t) if U < 0 or I < 0 or V < 0: #Si l'une des valeurs devient négative, c'est qu'il y a un bug !! #breakEt l'indentation ?
Cordialement,
Rescassol
Il faut plutôt écrire dans ton fichier toute la liste de valeurs à la fin de la boucle.
Essaie ça :
import matplotlib.pyplot as plt from numpy import log10 #L'avantage de la fonction log de numpy, c'est qu'elle est déjà vectorisée... import json UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0. V=4.5 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.0001 ntau=int(tau/h) #temps final T=25. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) for n in range(M) : Vn = (V + (0 if n <= ntau else h*b*II[n-ntau]))/(1+h*sigma) Un = U/(1 + a*h*Vn) In = (I + a*h*Un*Vn)/(1+beta*h) U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) #print(t,log10(V)) t=t+h temps.append(t) if U < 0 or I < 0 or V < 0: #Si l'une des valeurs devient négative, c'est qu'il y a un bug !! break with open('log.json', 'w') as f: json.dump(log10(VV).tolist(), f)mais il y a un souci. C'est que je souhaite avoir un fichier avec deux colonnes: une colonne avec les valeurs de t et une autre colonne avec les valeurs de log10(V(t)). Mais avec cette méthode o obtient un fichier avec uniquement les valeurs de log10(V). Est-ce qu'il y a une solution à ce problème? Stp
J'ai une dernière question. Comment il est possible d'améliorer la précision du code ci-dessous? (Que Bissam m'a aidé à écrire). A partir de t=5 , les valeurs sont supérieures à ce qu'elles devraient être de 0.2. J'ai essayé en diminuant le $h$ mais sans succès.
import matplotlib.pyplot as plt from numpy import log10 #L'avantage de la fonction log de numpy, c'est qu'elle est déjà vectorisée... import json UU=[] II=[] VV=[] temps=[] #initialisation U=1. I=0. V=4.5 tau=2 a=0.01 b=80000 beta=0.1 sigma=1 t=0. temps.append(t) #pas de temps h=0.001 ntau=int(tau/h) #temps final T=25. #M=nb de pas de temps M=int(T/h)+1 UU.append(U) II.append(I) VV.append(V) for n in range(M) : Vn = (V + (0 if n <= ntau else h*b*II[n-ntau]))/(1+h*sigma) Un = U/(1 + a*h*Vn) In = (I + a*h*Un*Vn)/(1+beta*h) U=Un I=In V=Vn UU.append(U) II.append(I) VV.append(V) print(t,log10(V)) t=t+h temps.append(t) if U < 0 or I < 0 or V < 0: #Si l'une des valeurs devient négative, c'est qu'il y a un bug !! breakPourquoi vas-tu chercher Gnuplot au lieu de faire directement ton graphique avec Python ?
Cordialement,
Rescassol
Sinon tu as une idée de comment améliorer la précision de mon code? Stp
Tu peux faire toutes tes courbes, sur le même graphique ou non, et les comparer, avec Python.
Pour la précision, essaie d'autres valeurs du pas que $h=0.001$.
De toutes façons, tu emploies une méthode approchée, il y aura toujours des erreurs.
Cordialement,
Rescassol
import numpy as np import os path = PATH = str(os.path.abspath('')) n = 100 print(f"log = {np.log10(n)}") fichier = 'toto.txt' tableau = np.random.random( (n,5) ) np.savetxt(path + '/' + fichier, tableau, delimiter=';')