x composé de M+2 blocs identiques chacun comportant xi de i = 0 à i = N+2
y composé de M+2 blocs chacun comportant N+2 fois la valeur tj = j * Delta t
Z composé M+2 blocs chacun comportant les valeurs de u(xi,tj) pour tj fixé et i de 0 à N+2
Ces vecteurs seront utilisés pour le traçage du plot 3D de u(x,t).
Voici mon code :
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
# Saisie des inputs par l'utilisateur
c = -1
e = 1
L = 1
T = 1
U0t = 2
ULt = 1
Ux0 = 1.5
N = input ("Saisissez le nombre des noeuds de l'espace")
N = int (N)
M = input ("Saisissez le nombre des noeuds du temps")
M = int (M)
# Calcul du pas de l'espace
h = L/(N+1)
# Discrétisation de l'intervalle [0,L] selon le pas h
X = np.linspace (0.0,L,N+2)
# Calcul du pas du temps
t = T/(M+1)
# Remplissage de A matrice tridiagonale symétrique de taille N*N
def tridiag(P, Q, R, k1=-1, k2=0, k3=1):
N = len(Q)
return (sps.spdiags(P,-1,N,N) + sps.spdiags(Q,0,N,N) + sps.spdiags(R,1,N,N))
P = (1/h**2)*np.ones(N)
Q = (c - 2/h**2)*np.ones(N)
A = sps.csc_matrix(tridiag(P,Q,P))
# Remplissage de B vecteur colonne de taille N
B = np.zeros(N)
B[0] = (1/h**2)*U0t
B[N-1] = (1/h**2)*ULt
# Initialisation de Uj à l'instant tj=0
U = Ux0*np.ones(N)
# Nouvelle A
A = sps.csc_matrix((e/t)*np.eye(N) - A)
# Intialisation de x, y = t et Z = u (x,t) à l'instant y = t = 0
x = X
y = np.zeros(N+2)
z = np.insert (U,0,[U0t])
z = np.append (z,[ULt])
Z = z
# Calcul de Uj+1 à partir de Uj
for j in range (1,M+2):
# Insertion du jème bloc des xi
x = np.append(x,X)
# Insertion du jème bloc de tj
y = np.append(y,[j*t*np.ones(N+2)])
# Calcul de Uj à partir de Uj-1
V = (e/t)*U + B
U = spsolve(A,V)
U = np.ravel(U)
# Insertion de u(0,tj) et u(L,tj)
z = np.append (U,[ULt])
z = np.insert (z,0,[U0t])
# Insertion du jème bloc de u(xi,tj)
Z = np.append (Z,z)
# Traçage du plot 3D de u(x,t)
Z = np.reshape(Z,((N+2)*(M+2),1))
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection="3d")
ax.plot_wireframe(x, y, Z, color='green')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')
ax.plot_surface (x, y, Z, rstride=1, cstride=1, cmap='winter', edgecolor='none')
plt.show()
Pour N = 10 et M =10, le code marche parfaitement. Mais, pour N = 100 et M =10, j'obtiens un Memory error.
Ton code marche sans problème, même pour des valeurs bien plus grandes de N et M (essayé avec N=10000, M=100) jusqu'à la ligne qui précède :
ax.plot_wireframe(x, y, Z, color='green')
que j'ai désactivée dans mon code ainsi que toute la suite.
Après c'est le rendu 3D qui est trop gourmand en ressources. Visiblement, x, y et Z ont chacun une taille de MxN donc avec N=100 et M=10 le code 3D de matplotlib va boucler sur 10**9 ce qui n'est pas hors d'atteinte mais matplotlib n'utilise pas OpenGl. Mais clairement c'est le rendu qui pose problème, il suffit déjà de voir le temps que ça prend pour N et M petits et de toute façon je ne suis même pas sûr que ce rendu aurait vraiment un sens à moins d'avoir une résolution très très élevée.
J'ai lu que pour faire de la visualisation 3D lourde, valait mieux utiliser d'autres outils comme Mayavi
Ton code marche sans problème, même pour des valeurs bien plus grandes de N et M (essayé avec N=10000, M=100) jusqu'à la ligne qui précède :
ax.plot_wireframe(x, y, Z, color='green')
que j'ai désactivée dans mon code ainsi que toute la suite.
Après c'est le rendu 3D qui est trop gourmand en ressources. Visiblement, x, y et Z ont chacun une taille de MxN donc avec N=100 et M=10 le code 3D de matplotlib va boucler sur 10**9 ce qui n'est pas hors d'atteinte mais matplotlib n'utilise pas OpenGl. Mais clairement c'est le rendu qui pose problème, il suffit déjà de voir le temps que ça prend pour N et M petits et de toute façon je ne suis même pas sûr que ce rendu aurait vraiment un sens à moins d'avoir une résolution très très élevée.
J'ai lu que pour faire de la visualisation 3D lourde, valait mieux utiliser d'autres outils comme Mayavi
- Edité par PascalOrtiz il y a environ 1 heure
Bonjour,
Merci pour ta réponse.
J'ai désactivé les lignes dont t'as parlé et j'ai fait un essai avec N = 10000 et M = 100, mais j'obtiens toujours un Memory error.
Se pourrait-il que ça soit un problème de RAM (la mienne est de 4GB) ?
Pour ce qui est du plot 3D, je verrai du côté de Mayavi, merci pour l'info
Difficile de dire, a priori, si ton système n'est pas trop occupé, je dirais que ça devrait passer. Ci-dessous, le code que j'ai executé ainsi que l'occupation mémoire que me donne mprof
#!/usr/bin/python3
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
# Saisie des inputs par l'utilisateur
c = -1
e = 1
L = 1
T = 1
U0t = 2
ULt = 1
Ux0 = 1.5
N=10000
M=100
# Calcul du pas de l'espace
h = L/(N+1)
# Discrétisation de l'intervalle [0,L] selon le pas h
X = np.linspace (0.0,L,N+2)
# Calcul du pas du temps
t = T/(M+1)
# Remplissage de A matrice tridiagonale symétrique de taille N*N
def tridiag(P, Q, R, k1=-1, k2=0, k3=1):
N = len(Q)
return (sps.spdiags(P,-1,N,N) + sps.spdiags(Q,0,N,N) + sps.spdiags(R,1,N,N))
P = (1/h**2)*np.ones(N)
Q = (c - 2/h**2)*np.ones(N)
A = sps.csc_matrix(tridiag(P,Q,P))
# Remplissage de B vecteur colonne de taille N
B = np.zeros(N)
B[0] = (1/h**2)*U0t
B[N-1] = (1/h**2)*ULt
# Initialisation de Uj à l'instant tj=0
U = Ux0*np.ones(N)
# Nouvelle A
A = sps.csc_matrix((e/t)*np.eye(N) - A)
# Intialisation de x, y = t et Z = u (x,t) à l'instant y = t = 0
x = X
y = np.zeros(N+2)
z = np.insert (U,0,[U0t])
z = np.append (z,[ULt])
Z = z
# Calcul de Uj+1 à partir de Uj
for j in range (1,M+2):
# Insertion du jème bloc des xi
x = np.append(x,X)
# Insertion du jème bloc de tj
y = np.append(y,[j*t*np.ones(N+2)])
# Calcul de Uj à partir de Uj-1
V = (e/t)*U + B
U = spsolve(A,V)
U = np.ravel(U)
# Insertion de u(0,tj) et u(L,tj)
z = np.append (U,[ULt])
z = np.insert (z,0,[U0t])
# Insertion du jème bloc de u(xi,tj)
Z = np.append (Z,z)
# Traçage du plot 3D de u(x,t)
Z = np.reshape(Z,((N+2)*(M+2),1))
Je viens de réaliser que j'avais installé la version 32 bits de Python (alors que j'ai un pc sous Windows 64 bits), ça marche bien maintenant pour N = 10000 et M = 100 (sans plot 3D) avec la version 64 bits. Merci pour ton temps
× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
Découverte Python Doc Tkinter Les chaînes de caractères
Découverte Python Doc Tkinter Les chaînes de caractères
Découverte Python Doc Tkinter Les chaînes de caractères
Découverte Python Doc Tkinter Les chaînes de caractères