Résolution d'un système d'edo avec python

ccapucine · October 2021

Bonjour Paul
Tu peux m’aider à améliorer mon code? Stp
Ou bien Bissam ou quelqu’un d’autre?
À partir de t=5 il donne des valeurs à 0.2 au dessus. Je ne comprends pas pourquoi.
Merci beaucoup d’avance.

paul18 · October 2021

Je n'ai fait que reprendre le code en essayant de l'adapter pour supprimer les "append", donc je ne maîtrise pas trop le résultat et pas sûr que cela soit juste.
J'espère toutefois qu'il y a des "choses" utiles à reprendre

import matplotlib.pyplot  as plt
import numpy as np
import os

path = PATH = str(os.path.abspath(''))

# Constantes
tau=2.
a=0.01
b=80000.
beta=0.1
sigma=1.

# 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

# pre-dimensionnement du tableau & initialisation
U = np.zeros(M+1)
I = np.copy(U)
V = np.copy(U)
t = np.copy(U)

U[0]=1.
I[0]=0.
V[0]=4.5

for n in range(1,M+1): # on commense au premier
    if n>ntau:
        V[n] = (V[n-1] + h*b*I[n-ntau])/(1+h*sigma)
    else:
        V[n] = V[n-1]/(1+h*sigma) 
    U[n] = U[n-1]/(1 + a*h*V[n-1])
    I[n] = (I[n-1] + a*h*U[n]*V[n])/(1+beta*h)
    t[n] = t[n-1] + h

# Figure
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(t, U, color = 'b', linewidth=1.5, linestyle ='-', label = r"$U$")
# ax.plot(t, V, color = 'g', linewidth=1.5, linestyle ='-', label = r"$V$")
ax.plot(t, I, color = 'r', linewidth=1.5, linestyle ='-', label = r"$I$")

ax.minorticks_on()
ax.grid(which='major', linestyle='-', linewidth='1.', color='black')
ax.grid(which='minor', linestyle=':', linewidth='0.5', color='blue')
ax.set_xlabel(r"t", fontsize=20)
ax.set_ylabel(r"??",  fontsize=20)

# plt.legend(loc="lower right", title=r"Legend", frameon=True, fontsize=20)
plt.legend(loc="upper right", title=r"Legend", frameon=True, fontsize=14, facecolor = "white", edgecolor = "black")

plt.savefig(PATH + '/plot.png',dpi=300)

fichier = 'toto.txt'
tableau = np.transpose(np.vstack((t,U,V,I)))
np.savetxt(path + '/' + fichier, tableau, delimiter=';')

ccapucine · October 2021

@paul18. Stp pour calculer log10(V) qu'est ce que vous suggérez? et afficher une colonne pour les valeurs de t et une autre pour les valeur de log10(V)

paul18 · October 2021

V étant déjà un vecteur, il suffit de faire

logV = np.log10(V)

ccapucine · October 2021

Pardon. Ca n'affiche pas toute les valeurs de log10(V), ca met des pointillés. Est-ce que c'est possible d'avoir une colonne pour t et une autre pour log10V? Stp

paul18 · October 2021

ccapucine a écrit:

ca met des pointillés + colonne pour t et une autre pour log10V

désolé mais je ne comprends pas : où les pointillés ? et ajouter où les colonnes ?

paul18 · October 2021

c'est en dehors de la boucle qu'il faut calculer $log10(V)$; le reste du code reste inchangé

for n in range(1,M+1): # on commence à t>0
    if n>ntau:
        V[n] = (V[n-1] + h*b*I[n-ntau])/(1+h*sigma)
    else:
        V[n] = V[n-1]/(1+h*sigma) 
    U[n] = U[n-1]/(1 + a*h*V[n-1])
    I[n] = (I[n-1] + a*h*U[n]*V[n])/(1+beta*h)
    t[n] = t[n-1] + h

    
# passage de V à log10(V)
V = np.log10(V)

ensuite tu mets ce que tu veux dans ton "tableau)

tableau = np.transpose(np.vstack((t,U,V,I)))
ou juste
tableau = np.transpose(np.vstack((t,V)))

ccapucine · October 2021

Voici une capture de ce qui s'affiche, ce ne sont pas toutes les valeurs qui sont affichées 128030

paul18 · October 2021

ah ok

Ton éditeur (IDE) ne peut pas afficher les 250_000 lignes (ou colonnes), donc il ajoute des pointillés et n'affiche que les premières et dernières x valeurs. Pour t'en convaincre, il suffit d'éditer le fichier txt qui a été créé.

Personnellement j'utilise Spyder comme IDE (sous Anaconda), et Vi comme éditeur, mais Excel fera très bien son office ici puisque c'est un format genre "csv"

ccapucine · October 2021

On obtient exactement les mêmes résultats.
A partir de t=5 c'est plus grand de 0.2. Je suis complétement perdue et je ne sais plus quoi faire pour améliorer le code et enlever ces 0.2 en plus