-
Notifications
You must be signed in to change notification settings - Fork 0
/
direct.py
105 lines (84 loc) · 2.86 KB
/
direct.py
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
97
98
99
100
101
102
103
104
105
# -*- encoding: utf-8 -*-
# ligne précedente pour pouvoir avoir les accents dans les commentaires
import numpy as np
from scipy import sparse
import scipy
from scipy.sparse.linalg import spsolve
pi = np.pi
#### fonction construisant la matrice et le second membre ####
def systeme(ksolve,ksch,N):
if ksolve <= 4: #matrice ordinaire
A = np.zeros([N * N, N * N])
else: #ksolve=5: matrice creuse, le format "lil" est le plus rapide en remplissage
A = sparse.lil_matrix((N * N, N * N)) # ne marche pas avec []
#
# Generation du maillage selon les axes x et y
#
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
if ksch == 1: # Remplissage dans le la matrice A dans le cas du schéma à 5 points
for i in range(0, N):
A[i, i] = 1
for i in range(1, N-1):
ii = i*N
A[ii, ii] = 1
for j in range(1, N-1):
ii = i*N + j
A[ii, ii] = 1
A[ii, ii - 1] = -1. / 4.
A[ii, ii + 1] = -1. / 4.
A[ii, ii - N] = -1. / 4.
A[ii, ii + N] = -1. / 4.
ii = i*N + N - 1
A[ii, ii] = 1
for i in range(N*(N-1)-1, N * N):
A[i, i] = 1
else:
for i in range(0, N):
A[i, i] = 1
for i in range(1, N-1):
ii = i*N
A[ii, ii] = 1
for j in range(1, N-1):
ii = i*N + j
A[ii, ii] = 1
A[ii, ii - 1] = -1. / 5.
A[ii, ii + 1] = -1. / 5.
A[ii, ii - N] = -1. / 5.
A[ii, ii + N] = -1. / 5.
A[ii, ii + N + 1] = -1. / 20.
A[ii, ii + N - 1] = -1. / 20.
A[ii, ii - N + 1] = -1. / 20.
A[ii, ii - N - 1] = -1. / 20.
ii = i*N + N - 1
A[ii, ii] = 1
for i in range(N*(N-1)-1, N * N):
A[i, i] = 1
# Remplissage du vecteur b
b = np.zeros(N * N)
for i in range(N):
ii = i*N + N -1
b[ii] = np.sin(pi * x[i])
for i in range(N):
ii = (N - 1) * N + i
b[ii] = np.sin(pi * y[i])
return A,b
#### fonction résolvant le systà¨me ####
def direct(ksolve, ksch , N, Te): #
A,b = systeme(ksolve,ksch,N)
# Résolution du systà¨me par méthode LU
#matrice ordinaire:
if ksolve == 4:
resu = np.linalg.solve(A, b)
#matrice creuse:
else:
resu = scipy.sparse.linalg.spsolve (scipy.sparse.csr_matrix(A), b)
# Remise en forme de la solution
T = np.zeros([N, N])
ii = 0
for i in range(0,N):
for j in range(0,N):
T[i, j] = resu[ii]
ii = ii + 1
err = np.linalg.norm(T-Te)/N
return (err,T)