Forum de mathématiques - Bibm@th.net
Vous n'êtes pas identifié(e).
- Contributions : Récentes | Sans réponse
#1 18-10-2021 18:57:06
- ccapucine
- Membre
- Inscription : 19-05-2018
- Messages : 195
Résoudre un système d'edo avec Python
Bonjour,
je souhaite résoudre le système d'edo ci dessous en utilisant Python
et avoir le graphe de $log(V)$ où $V$ est la solution de la troisième équation
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,
$$
où $u_0=1, v_0=0.45, a=0.01, b=80000, \beta=0.1, \sigma=1,\tau=2$
J'ai essayé le code suivant mais sans pouvoir faire le plot de log(V)
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')
plt.show()
Hors ligne
#3 19-10-2021 08:34:39
- Zebulor
- Membre expert
- Inscription : 21-10-2018
- Messages : 2 220
Re : Résoudre un système d'edo avec Python
Bonjour,
peut être que Fred n'a pas tout vu .. si on déroule l'escalier je vois un print à la fin ..
...
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()
Hors ligne
#4 19-10-2021 08:46:14
- Fred
- Administrateur
- Inscription : 26-09-2005
- Messages : 7 348
Re : Résoudre un système d'edo avec Python
Salut Zebulor,
Si, si, j'avais vu (c'est d'ailleurs pour cela que je demandais si les valeurs obtenues pour U,I et V étaient cohérentes).
Je que je voulais dire, c'est que le plt.show() ne sert à rien...
F.
Hors ligne
#5 22-10-2021 10:46:56
- ccapucine
- Membre
- Inscription : 19-05-2018
- Messages : 195
Re : Résoudre un système d'edo avec Python
Bonjour,
voici le code. L'objectif étant d'afficher les valeurs de log10(V)
Le souci est que, à partir de t=5, ce code donne des valeurs à 0.2 au dessus.
Je n'arrive pas à régler ce problème. Je vous remercie d'avance de m'aider à améliorer ce code
import matplotlib.pyplot as plt
import numpy as np
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.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 0.1*h*b*II[n-ntau]))/(1+0.1*h*sigma)
Un = U/(1 + a*0.1*h*Vn)
In = (I + a*0.1*h*Un*Vn)/(1+beta*0.1*h)
U=Un
I=In
V=Vn
UU.append(U)
II.append(I)
VV.append(V)
print(t,np.log10(V))
#print(t,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
Hors ligne