Discussion :

[HanaChan]

Bonjour,

J'ai besoin de résoudre une équation aux dérivées partielle par la méthode des différences finies : https://drive.google.com/file/d/1QQt...ew?usp=sharing

Je choisis comme a(x) la fonction carré.

Voici mon code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
import os
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
from scipy import misc
*
# Saisie des inputs par l'utilisateur
c = input ("Saisissez la valeur de c")
c = float (c)
e = input ("Saisissez la valeur de epsilon")
e = float (e)
L = input ("Saisissez la longueur de l'intervalle [0,L]")
L = float (L)
T = input ("Saisissez la longueur de l'intervalle [0,T]")
T = float (T)
U0t = input ("Saisissez la valeur de u(0,t)")
U0t = float (U0t)
ULt = input ("Saisissez la valeur de u(L,t)")
ULt = float (ULt)
Ux0 = input ("Saisissez la valeur de u(x,0)")
Ux0 = float (Ux0)
N = input ("Saisissez le nombre des noeuds de l'espace")
N = int (N)
M = input ("Saisissez le nombre des noeuds de 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)
*
# Discrétisation de l'intervalle [0,T] selon le pas t
Y = np.linspace (0.0,T,M+2)
*
*
# Fonction a(x)
def a(x):
****return(x**2)
*
# Remplissage de la matrice tridiagonale A 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 = [- misc.derivative(a,(i)*h)/(2*h) + a((i)*h)/h**2 for i* in range(2,N+2)]
Q = [c - 2*a((i+1)*h)/h**2 for i in range(N)]
R = [misc.derivative(a,(i)*h)/(2*h) + a((i)*h)/h**2 for i in range(N)]
*
A = sps.csc_matrix(tridiag(P,Q,R))
print (A)
*
# Remplissage du vecteur colonne B de taille N
B = np.zeros(N)
B[0] = (- misc.derivative(a,h)/(2*h) + a(h)/h**2)*U0t
B[N-1] = (misc.derivative(a,N*h)/(2*h) + a(N*h)/h**2)*ULt
print (B)
*
# Initialisation de Uj à l'instant tj=0
U = Ux0*np.ones(N)
print(U)
*
# Déclaration de la matrice identité de taille N*N
I = np.eye(N)
*
# Nouvelle A
A = (t/e)*A + I
*
# Nouveau B
B = (t/e)*B
*
# Calcul de Uj+1 à partir de Uj
*
for j in range (1,M+2):
****U = np.dot(A,U)+ B
****print(U)
*
os.system("pause")

Le problème réside dans la boucle For (ou plus précisément dans le calcul de Uj+1 à partir de Uj).

Lorsque j'ai fait un essai pour c = -1, epsilon = 1, longueur de l'intervalle [0,L] = 1, longueur de l'intervalle [0,T] = 1, u(0,t)= 2, u(L,t) = 1, u(x,0) = 1.5, Nombre des nœuds de l'espace N = 2, Nombre des nœuds du temps M = 1, le programme demande tous les inputs puis se ferme rapidement.

Lorsque j'ai fait un autre essai (toujours les mêmes inputs) mais sans boucle (je voulais juste calculer U1 à l'instant tj = t1), le programme marche bien et affiche la valeur correcte de U1. Mais lorsque je retape le même code pour calculer U2 à partir de U1, le programme demande juste les inputs puis se ferme !

Je ne comprends pas non plus pourquoi la boucle n'affiche même pas la première itération sachant qu'elle est correcte.

Merci de votre aide

[wiztricks]

Salut,

Je ne comprends pas non plus pourquoi la boucle n'affiche même pas la première itération sachant qu'elle est correcte.

Probable que vous lanciez votre programme avec un double-clic sur le fichier depuis l'explorateur.
S'il arrive jusqu'à l'instruction os.system('pause'), vous allez voir quelque chose.
S'il se plante avant, vous ne voyez même pas le message d'erreur.

Créez un terminal Windows (avec la commande cmd.exe), vérifier comment lancer l'interpréteur Python, et lancez votre programme avec la commande $ python programme.py. En remplaçant python par la bonne commande

- W

[HanaChan]

Salut,



Probable que vous lanciez votre programme avec un double-clic sur le fichier depuis l'explorateur.
S'il arrive jusqu'à l'instruction os.system('pause'), vous allez voir quelque chose.
S'il se plante avant, vous ne voyez même pas le message d'erreur.

Créez un terminal Windows (avec la commande cmd.exe), vérifier comment lancer l'interpréteur Python, et lancez votre programme avec la commande $ python programme.py. En remplaçant python par la bonne commande

- W

Bonjour,

Merci pour votre réponse J'ai fait ce que vous avez dit et c'est ce que j''ai obtenu :


J'ai ajouté U = np.reshape (U,(N,1)) et B = np.reshape (B,(N,1)) après le calcul de U à l'itération j à partir de U à l'itération j-1, ça marche bien maintenant.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
# Calcul de Uj+1 à partir de Uj
 
for j in range (1,M+2):
	U = np.dot(A,U)+ B
	print(U)
	U = np.reshape (U,(N,1))
	B = np.reshape (B,(N,1))