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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
| 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() |
