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