Python pour une edo
Bonjour,
j'essaye de programmer le système d'edo suivant avec Python:
\begin{align*}
U'(t)&=-a_1 U V_1 - a_2 U V_2;
\\
I_1'(t)&=a_1 U V_1 - \beta_1 I_1;
\\
I_2'(t)&=a_2 U V_2 - \beta_2 I_2;
\\
V_1'(t)&=b_1 I_{1_{\tau_1}}-\sigma_1 V_1;
\\
V_2'(t)&= \varepsilon I_{1_{\tau_1}}+b_2 I_{2_{\tau_2}}-\sigma_2 V_2,
\end{align*} L'objectif est de calculer $\log(V_1)$ et $\log(V_2)$ où $V_1$ et $V_2$ sont les solutions des deux dérnières équations.
Voici ce que j'ai écrit, mais j'ai l'erreur suivante.
RuntimeWarning: divide by zero encountered in log10 W2= np. log10(V2)
Je ne comprends pas cette erreur.
Voici le code.
Je vous remercie d'avance
j'essaye de programmer le système d'edo suivant avec Python:
\begin{align*}
U'(t)&=-a_1 U V_1 - a_2 U V_2;
\\
I_1'(t)&=a_1 U V_1 - \beta_1 I_1;
\\
I_2'(t)&=a_2 U V_2 - \beta_2 I_2;
\\
V_1'(t)&=b_1 I_{1_{\tau_1}}-\sigma_1 V_1;
\\
V_2'(t)&= \varepsilon I_{1_{\tau_1}}+b_2 I_{2_{\tau_2}}-\sigma_2 V_2,
\end{align*} L'objectif est de calculer $\log(V_1)$ et $\log(V_2)$ où $V_1$ et $V_2$ sont les solutions des deux dérnières équations.
Voici ce que j'ai écrit, mais j'ai l'erreur suivante.
RuntimeWarning: divide by zero encountered in log10 W2= np. log10(V2)
Je ne comprends pas cette erreur.
Voici le code.
Je vous remercie d'avance
import matplotlib.pyplot as plt import numpy as np import os path = PATH = str(os.path.abspath('')) # Constantes tau=2. a1=10. a2=10. b1=5. b2=5. beta1=0. beta2=0. sigma1=0.4 sigma2=0.4 eps=0.1 # 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 # pre-dimensionnement du tableau & initialisation U = np.zeros(M+1) I1 = np.copy(U) I2 = np.copy(U) V1 = np.copy(U) V2 = np.copy(U) t = np.copy(U) U[0]=1. I1[0]=0. I2[0]=0. V1[0]=1. V2[0]=0. for n in range(1,M+1): # on commense au premier if n>ntau: V1[n] = (V1[n-1] + h*b1*I1[n-ntau])/(1+h*sigma1) V2[n]= (V2[n-1]+eps*h*I1[n-ntau]+b2*h*I2[n-ntau])/ (1+sigma2*h) else: V1[n] = V1[n-1]/(1+h*sigma1) V2[n]= V2[n-1]/(1+h*sigma2) U[n]=U[n-1]/(1+a1*h*V1[n]+a2*h*V2[n]) I1[n] = (I1[n-1] + a1*h*U[n]*V1[n])/(1+beta1*h) I2[n] = (I2[n-1] + a2*h*U[n]*V2[n])/(1+beta2*h) t[n] = t[n-1] + h #fichier = 'u.txt' #tableau = np.transpose(np.vstack((t,U))) #np.savetxt(path + '/' + fichier, tableau, delimiter=' ') # passage de V1 à log10(V1) W1 = np.log10(V1) # passage de V2 à log10(V2) W2 = np.log10(V2) #fichier = 'i.txt' #tableau = np.transpose(np.vstack((t,I))) #np.savetxt(path + '/' + fichier, tableau, delimiter=' ') #fichier = 'v.txt' #tableau = np.transpose(np.vstack((t,V))) #np.savetxt(path + '/' + fichier, tableau, delimiter=' ') fichier = 'w1.txt' tableau = np.transpose(np.vstack((t,W1))) np.savetxt(path + '/' + fichier, tableau, delimiter=' ') fichier = 'w2.txt' tableau = np.transpose(np.vstack((t,W2))) np.savetxt(path + '/' + fichier, tableau, delimiter=' ')
Réponses
-
Bonjour.
De manière générale, les logarithmes n'aiment pas quand leur argument est nul, ce qui est le cas pour V2 à son initialisation, si j'ai bien lu.
À bientôt.Cherche livres et objets du domaine mathématique :
Intégraphes, règles log et calculateurs électromécaniques.
-
Oui mais c'est ça la valeur initiale. On ne peut pas la changer.
Que faire dans ce cas? Comment modifier l'algorithme pour que ça marche?
Merci d'avance. -
Si je ne m'abuse, c'est une résolution numérique.
Quid de prendre une valeur initiale de l'ordre de $10^{-6}$ en simple précision ou $10^{-12}$ en double précision ?
Si vraiment on ne veut pas avoir à chipoter avec cela, il doit aussi y avoir moyen de reprendre eps...
Une façon de faire :def machineEpsilon(func=float): machine_epsilon = func(1) while func(1)+func(machine_epsilon) != func(1): machine_epsilon_last = machine_epsilon machine_epsilon = func(machine_epsilon) / func(2) return machine_epsilon_last
>>> import numpy >>> machineEpsilon(numpy.float64)
À bientôt.Cherche livres et objets du domaine mathématique :
Intégraphes, règles log et calculateurs électromécaniques.
-
Je ne comprends pas que fait cette partie du code? Où est ce qu'on la met?Car ça ne fonctionne pas
-
pour $V_2$ de 0 à ntau, le numérateur est toujours égal à zéro, donc $log(V_2)$ sera toujours impossible. On peut "biaiser" en remplaçant les valeurs nulles par une autre ($10^{-12}$ par exemple) avant de passer en log ce qui évitera le problème.
Par exemple:# passage de V1 à log10(V1) ValeursNulles = np.where(V1 == 0.) V1[ValeursNulles] = 1.E-12 W1 = np.log10(V1) # passage de V2 à log10(V2) ValeursNulles = np.where(V2 == 0.) V2[ValeursNulles] = 1.E-12 W2 = np.log10(V2)
-
J'aurais dû préciser que le code donné fournit juste la plus petite valeur positive pour l'implémentation machine. En bref c'est juste pour avoir la valeur la plus petite acceptable par le log.
La recommandation de paul18 est plus simple.
À bientôt.Cherche livres et objets du domaine mathématique :
Intégraphes, règles log et calculateurs électromécaniques.
-
pour ma culture personnelle, j'ai cherché l'information pour numpy :
np.finfo(np.float64).eps Out[1]: 2.220446049250313e-16 np.finfo(np.float64).tiny Out[2]: 2.2250738585072014e-308 np.log10(np.finfo(np.float64).tiny) Out[3]: -307.6526555685888
-
Merci pour l'info.
Comme quoi j'étais encore loin du compte.
À bientôt.
Cherche livres et objets du domaine mathématique :
Intégraphes, règles log et calculateurs électromécaniques.
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
- 65 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