Calcul Scientifique TP1 -Corrig?

September 29, 2019

1 TP1 : Résolution numérique de systèmes d’équations linéaires

1.1 Objectif
L’objectif de ce TP est de pouvoir programmer les algorithmes du pivot de Gauss, Jacobi, et
Gauss-Seidel pour résoudre des Systèmes d’Equations Linéaires (SEL) et ceci en utilisant la bib-
liothèque numpy.

1.2 Rappel
On considère le SEL 

 a11 x1 + a12 x2 + . . . + a1n xn = b1

 a21 x1 + a22 x2 + . . . + a2n xn = b2
(S) : ..

 .


an1 x1 + an2 x2 + . . . + ann xn = bn
L’écriture matricielle de (S) est donnée par AX = b, avec :
     
a11 a12 a13 · · · a1n x1 b1
 a21 a22 a23 · · · a2n   x2   b2 
     
     
A =  a31 a32 a33 · · · a3n  , X =  x3  , et b =  b3  .
 .. .. .. .. ..   ..   .. 
 . . . . .  . .
an1 an2 an3 · · · ann xn bn
Il existe deux catégories de méthodes pour la résolution du système (S) :

• les méthodes directes : la solution X est exacte, et donnée après un nombre d’opérations
algorithmiques fixé en fonction de n, la taille de la matrice A.

Exemples : la méthode du pivot de Gauss, la décomposition LU, . . .

• les méthodes itératives (ou indirectes) : la solution X est approchée et donnée une fois qu’une
condition d’arrêt, imposée à l’algorithme, est satisfaite.

Exemples : Jacobi, Gauss-Seidel, . . .

1
 2 Méthodes directes
2.1 Description de l’algorithme du pivot de Gauss
On suppose que le système linéaire AX = b, avec A = ( aij )1≤i,j≤n , X = ( xi )1≤i≤n et b = (bi )1≤i≤n ,
est de Cramer, c’est-à-dire il admet une unique solution. La matrice A est donc inversible.
L’algorithme du pivot de Gauss, consiste à triangulariser la matrice A à l’aide d’opérations élé-
mentaires.
En notant Li la ime ligne de A, 1 ≤ i ≤ n, on appelle opération élémentaire sur A toute opération
de l’un des types suivants:

Li ←− Li + λL j , λ ∈ R∗ ,
Li ←− λLi ,
Li ←→ L j .

2.2 Algorithme de pivot de Gauss

Pour résoudre le système AX = b, on suit les trois étapes suivates :

• Etape 1 : on construit la matrice élargie ( A|b) de type (n, n + 1), en concaténant la matrice
A et le second membre b comme suit :
 
a11 · · · a1n b1
 .. 
( A|b) =  ... ..
. . 
an1 · · · ann bn

• Etape 2 : on effectue des opérations élementaires sur les lignes de ( A|b) jusqu’à arriver à un
système triangulaire.
• Etape 3 : on remonte par substitutions pour obtenir la solution X du système.

2.2.1 Etape 1 : Matrice élargie

Ecrire une fonction concatener(A,b) prenant en argument une matrice carrée A de taille n et
une matrice colonne b de taille n également et renvoyant la matrice élargie ( A|b) correspondante.
Compléter le code suivant :
[1]: import numpy as np
def concatener(A,b):
n=len(b)
Ab=np.zeros((n,n+1))
Ab[0:n,0:n]=A
Ab[0:n,n]=b[:,0]
return Ab
Tester la fonction concatener(A,b) sur le SEL du TP0 (S0 ) : AX = b, avec :
     
4 −1 −1 0 x1 50
 −1 4 0 − 1   x2  30
A=  −1 0
 , X =   , et b =   ,
4 −1  x3  70
0 −1 −1 4 x4 50

2
 :
[2]: A=5*np.eye(4)-np.ones((4,4))
[A[0,3],A[1,2],A[2,1],A[3,0]]=np.zeros(4)
b=np.array([[50,30,70,50]]).transpose()

2.2.2 Etape 2 : Triangulariser le système AX = b

a- Opération élémentaire : Permutation de deux lignes Ecrire une fonction permuter(A,i,j)
qui permute la ième et la jème colonnes de la matrice A.
[3]: def permuter(A,i,j):
'''Échange des lignes i et j de la matrice A.'''
P=A.copy()
L=P[i,:].copy() # sans .copy() L est un alias de P[i,:]
P[i,:]=P[j,:]
P[j,:]=L
return P
Tester la fonction permuter(A,i,j) sur A. Le résultat trouvé est-il celui attendu?
[ ]:

b- Opération élémentaire : Combinaison linéaire de lignes Ecrire une fonction

ajouter(A,i,j,alpha) qui ajoute, à la ligne Li , une multiple α de la ligne L j :

Li ←− Li + αL j

[4]: def ajouter(A,i,j,alpha):

'''Ajouter à la ligne i, alpha * la ligne j'''
C=A.copy()
C[i,:]+=alpha*C[j,:]
return C
Tester la fonction ajouter(A,i,j,alpha) sur A.
[ ]:

c- Recherche d’un pivot non nul Ecrire une fonction choix_pivot(A,j) qui, si le pivot de l’étape
j est nul, 1 ≤ j ≤ n − 1, cherche le premier élément non nul de a j+1,j , a j+2,j ,. . . ,an,j et renvoie l’indice
de la ligne correspondante. Sinon, elle retourne j. Compléter le code ci-dessous en utilisant la
fonction np.where. Utiliser le help de python pour se documenter sur l’utilité de cette fonction.
[5]: def choix_pivot(A,j):
n=A.shape[0]
if A[j,j]!=0:
Ind_max=j
else:
Ind=np.where(A[j+1:n,j]!=0)
Ind_max=Ind[0][0]+j+1
return Ind_max

3
 Tester la fonction choix_pivot(A,j) sur B=A.copy() avec B[1,1]=0 et pour j=1.
[ ]:

d- Recherche du pivot partiel Ecrire une fonction choix_pivot_partiel(A,j) qui, pour l’étape
j, 1 ≤ j ≤ n − 1, cherche le maximum des éléments a j,j , a j+1,j ,. . . ,an,j en valeur absolue et renvoie
l’indice de la ligne correspondante.
[6]: def choix_pivot_partiel(A,j):
n=A.shape[0]
return np.argmax(np.abs(A[j:n,j]))+j
Tester la fonction choix_pivot_partiel(A,j).
[ ]:

e- Triangularisation du système Ecrire une fonction triangulariser qui triangularise une ma-
trice A par la méthode du pivot de Gauss. Compléter le code suivant :
[7]: def triangulariser(A):
n=A.shape[0]
T=A.copy()
for j in np.arange(0,n):
ind_pivot=choix_pivot(T,j)
if ind_pivot>j:
T=permuter(T,j,ind_pivot)
for i in np.arange(j+1,n):
T=ajouter(T,i,j,-T[i,j]/T[j,j])
return T
Tester la fonction triangulariser(A) sur A.
[ ]:

2.2.3 Etape 3 - Résolution d’un système triangulaire par remontée

Cas où A est triangulaire supérieure.

 bn

 xn = a ;
nn
n
 1

 x i =
aii
( b i − ∑ aij x j ), ∀i = n − 1, n − 2, . . . , 1.
j = i +1

Ecrire une fonction remonter()qui résout un SEL AX = b par remotée, avec A une matrice
triangulaire supérieure.
[8]: def remonter(A):
n=A.shape[0]
x=np.zeros((n,1))
x[n-1]=A[n-1,n]/A[n-1,n-1]
for i in np.arange(n-1,-1,-1):
S=A[i,i+1:n].dot(x[i+1:n])

4
 x[i]=(A[i,n]-S)/A[i,i]
return x

2.2.4 Synthèse
Ecrire une fonction PivotDeGauss(A,b)qui résout un SEL AX = b par la méthode du pivot de
Gauss. Completer le code suivant :
[9]: def PivotDeGauss(A,b):
if np.linalg.det(A)==0:
print('A n est pas inversible')
else :
Ab=concatener(A,b)
T=triangulariser(Ab)
X=remonter(T)
return X
Appliquer la fonction PivotDeGauss(A,b) pour résoudre le SEL (S0 ).
[ ]:

2.2.5 Comparaison PivotDeGauss(A,b) et solve(A,b) en terme de temps de calcul

[10]: %timeit PivotDeGauss(A,b)

%timeit np.linalg.solve(A,b)
# 1microsecond=10-^6 seconds

159 ţs ś 10.9 ţs per loop (mean ś std. dev. of 7 runs, 10000 loops each)
16 ţs ś 3.46 ţs per loop (mean ś std. dev. of 7 runs, 100000 loops each)

3 Méthodes itératives
Contrairement aux méthodes directes (pivot de Gauss, LU,..) qui consistent à un gros calcul exacte,
les méthodes itératives ( Jacobi, Gauss-Seidel,. . . ) consistent en une succession de petit calcul
approché de la solution.

3.0.1 Principe
On considère le système linéaire AX = b de Cramer. Les méthodes itératives se basent sur une dé-
composition de A sous la forme A = M − N où M est une matrice inversible. Une suite récurrente
de solutions X (k) , k ≥ 0, est ensuite génére comme suit :

X (k+1) = M1 NX (k) + M−1 b

Cette suite converge, sous certaines conditions, vers la solution exacte X.

• Pour un choix de M = D, où D est une matrice diagonale dont les éléments diagonaux
correspondent à ceux de A, et N = D − A on se trouve avec la méthode de Jacobi.

5
 • Pour un choix de M = L, où L est une matrice triangulaire inférieure dont les éléments sont :
lij = aij si i ≥ j et lij = 0 sinon, et N = L − A on se trouve avec la méthode de Gauss-Seidel.

3.0.2 Convergence
Théorème: Si la matrice A est à diagonale strictement dominante alors les méthodes de Jacobi et
de Gauss_Seidel convergent quelque soit le choix du vecteur initial X (0) .
Ecrire une fonction qui vérifie si une matrice est à diagonale strictement dominante (etat 0) ou
non (etat 1).
[11]: def Matrice_diag_dominante(A):
for i in np.arange(0,A.shape[0]):
if np.abs(A[i,i])<=np.sum(np.abs(A[i,:]))-np.abs(A[i,i]):
#print('A n\'est pas à diagonale strictement dominante')
etat=0
break
else:
#print('A est à diagonale strictement dominante')
etat=1
return etat
Tester la fonction Matrice_diag_dominante(A) sur A.
[ ]:

3.0.3 Algorithme de Jacobi

Donner un vecteur initial X (0) et un réel ε > Tant que | X (k+1) − X (k) | > ε faire > > Pour i de 1 à n
faire > > >
1( n
(k) )
bi − ∑ aij x j
( k +1)
xi =
aii j =1
j ̸ =i

> > > > Fin pour > > Fin Tant que
Exercie
Ecrire une fonction jacobi(A, b, X0, epsilon) qui résout le SEL AX = b par la méthode de
Jacobi avec une précision ε et un vecteur initial X0. Proposer deux façons pour le faire.
[12]: # Méthode 1
def jacobi(A, b, X0, epsilon):
etat=int(Matrice_diag_dominante(A))
if etat==0:
print('A n\'est pas à diagonale strictement dominante')
return
else:
D = np.diagflat(np.diag(A))
N = D-A
X1=np.linalg.inv(D).dot(N.dot(X0)+b)
k=1
while np.linalg.norm(X1-X0,1)>epsilon:
X0=X1

6
 X1=np.linalg.inv(D).dot(N.dot(X0)+b)
k+=1
return X0,k
 
1
 1
Tester la méthode de Jacobi pour X (0) =   −6
1 et ε = 10 . Interpréter le résultat.
1
[13]: X0=np.ones((4,1))
epsilon=10**(-6)
jacobi(A, b, X0, epsilon)
[13]: (array([[24.99999964],
[19.99999964],
[29.99999964],
[24.99999964]]), 27)

3.1 Algorithme de Gauss-Seidel

Donner un vecteur initial X (0) et un réel ε > Tant que | X (k+1) − X (k) | > ε faire > >

1 ( n
(k) )
b1 − ∑ a1j x j
( k +1)
x1 =
a11 j =2

> > Pour i de 2 à n − 1 faire > > >

i −1
1( n
(k) )
bi − ∑ aij x j − ∑ aij x j
( k +1) ( k +1)
xi =
aii j =1 j = i +1

> > > > Fin pour > >

n −1
1 ( ( k +1) )
bn − ∑ a1j x j
( k +1)
xn =
ann j =1

> > Fin Tant que

Ecrire une fonction Gauss_Seidel(A, b, X0, epsilon)qui résout le SEL AX = b par la méth-
ode de Gauss-Seidel avec une précision ε et un vecteur initial X0. Proposer deux façons pour le
faire.
[14]: # Méthode 1
def Gauss_Seidel(A, b, X0, epsilon):
etat=int(Matrice_diag_dominante(A))
if etat==0:
print('A n\'est pas à diagonale strictement dominante')
else:
D = np.tril(A)
N = D-A
X1=np.linalg.inv(D).dot(N.dot(X0)+b)
k=1
while np.linalg.norm(X1-X0,1)>epsilon:
X0=X1

7
 X1=np.linalg.inv(D).dot(N.dot(X0)+b)
k+=1
return X0,k
 
1
 1
Tester la méthode de Gauss_Seidel pour X (0) =   −6
1 et ε = 10 . Interpréter le résultat.
1
[15]: X0=np.ones((4,1))
epsilon=10**(-6)
Gauss_Seidel(A, b, X0, epsilon)
[15]: (array([[24.99999973],
[19.99999987],
[29.99999987],
[24.99999993]]), 15)
Comparer le temps de calcul de chacune des méthodes de résolution de SEL appliquées sur
( S0 ) .
[16]: %timeit jacobi(A,b,X0,epsilon)
%timeit Gauss_Seidel(A,b,X0,epsilon)
%timeit np.linalg.solve(A,b)
%timeit PivotDeGauss(A,b)

2.07 ms ś 333 ţs per loop (mean ś std. dev. of 7 runs, 100 loops each)
720 ţs ś 136 ţs per loop (mean ś std. dev. of 7 runs, 1000 loops each)
14.4 ţs ś 1.6 ţs per loop (mean ś std. dev. of 7 runs, 100000 loops each)
164 ţs ś 29.7 ţs per loop (mean ś std. dev. of 7 runs, 1000 loops each)

Exercice
On considère le SEL (Sn ) :AX = b, avec
  

8 −1 −1 0 0 ··· 0 x1  
  6
 −1 8 −1 −1 0 ···  0 x2  5
   x3   
 −1 −1 8 −1 −1 ···  0  4
     
  .. 
A=
..
.
..
.
..
. 
,X =  . , b =  .
 ..  , et n ≥ 4.
     
 .. .. ..   ...  4
 . . .     
   ..  5
 −1 −1 8 −1 .
0 ··· ··· 0 −1 −1 8 xn
6

1- Pour ε = 10−6 , représenter sur un même graphe, le nombre d’itérations effectuées pour
atteindre la convergence par les deux méthodes itératives : Jacobi et Gauss-Seidel, en fonction de
n la taille de la matrice. On considère n ∈ {5, 10, 15, 20, 25, 30}.
2- Pour n = 20, représenter sur un même graphe, le nombre d’itérations effectuées pour at-
teindre la convergence par les deux méthodes itératives : Jacobi et Gauss-Seidel, en fonction de la
précision ε. On considère ε ∈ {10−3 , 10−6 , 10−9 , 10−12 }.
On considère le vecteur nul comme vecteur initial.
Interpréter les résultats.

8
[24]: def S_n(n):
A=8*np.eye(n)+np.diagflat(-np.ones(n-1),1)+np.diagflat(-np.ones(n-2),2)+np.
,→diagflat(-np.ones(n-1),-1)+np.diagflat(-np.ones(n-2),-2)

b=4*np.ones((n,1))
[b[0,0], b[-1,0]]=[6,6]
[b[1,0], b[n-2,0]]=[5,5]
return A,b
[25]: # question 1
epsilon=10**(-6)
N=np.arange(5,31,5)
Iter_J=[]
Iter_GS=[]
for n in N:
S=S_n(n)
X0=np.zeros((n,1))
J=jacobi(S[0],S[1],X0,epsilon)
GS=Gauss_Seidel(S[0],S[1],X0,epsilon)
Iter_J.append(J[1])
Iter_GS.append(GS[1])
[26]: import matplotlib.pyplot as plt
plt.figure(figsize=(20,10))
plt.plot(N,Iter_J,'ro--',N,Iter_GS,'b^--', linewidth=3,markersize=12)
plt.xlabel('Taille de la matrice',fontsize=30)
plt.ylabel('Nombre d\'itérations',fontsize=30)
plt.legend(('Méthode de Jacobi','Méthode de Gauss-Seidel'),fontsize=20, loc = 0)
plt.grid(True)
plt.savefig('Comparaison_n.png',format='png')

9
[27]: # question 2
n=20
epsilon=[10**-3,10**-6,10**-9,10**-12]
S=S_n(n)
X0=np.zeros((n,1))
Iter_J=[]
Iter_GS=[]
for eps in epsilon:
J=jacobi(S[0],S[1],X0,eps)
GS=Gauss_Seidel(S[0],S[1],X0,eps)
Iter_J.append(J[1])
Iter_GS.append(GS[1])
[28]: plt.figure(figsize=(20,10))
plt.plot(epsilon,Iter_J,'ro--',epsilon,Iter_GS,'b^--',␣
,→linewidth=3,markersize=12)

plt.xscale('log')
plt.xlabel('La précision',fontsize=30)
plt.ylabel('Nombre d\'itérations',fontsize=30)
plt.legend(('Méthode de Jacobi','Méthode de Gauss-Seidel'),fontsize=20, loc = 0)
plt.grid(True)
plt.savefig('Comparaison_n.png',format='png')

4 Application
On considère le circuit, donné ci-dessous, pour lequel on donne :
U1 = 110V, U2 = 105V, U3 = 90V, R1 = 0.5W, R2 = 0.25W, et R3 = 0.5W.
V et W désignent respectivement Volt et Watt.

10
[22]: plt.figure(figsize=(20,10))
plt.imshow(plt.imread('circuit.png'))
[22]: <matplotlib.image.AxesImage at 0x1bcbe0df630>

 
i1
Soit I = i2  , on souhaite résoudre un système de la forme RI = U, où R est une matrice
i3
dont les éléments se définissent avec R1 , R2 , et R3 et U est un vecteur dont les éléments se
définissent avec U1 , U2 , et U3 .
A l’aide des lois des nuds, des mailles et d’Ohm, on va déterminer R et U.
Pour calculer les courants i1 , i2 et i3 , on leur choisit un sens, ce qui permet de flécher les tensions
aux bornes des résistances. Le circuit ayant deux noeuds, il n’y a q’une équation de noeuds, en A
par exemple:

i1 + i2 = i3
Il y a 3 mailles, donc on choisit 2 mailles indépendantes afin de compléter le système à résoudre :
Maille 1 : A, R1 , U1 , B, U2 , R2 , A, avec le sens de parcours indiqué :

R1 .i1 − U1 + U2 − R2 .i2 = 0

Maille 2 : A, R2 , U2 , B, U3 , R3 , A avec le sens de parcours indiqué :

R2 .i2 U2 + U3 + R3 .i3 = 0

1. Déterminer R et U. 1. Détermine I en utilisant solve() et PivotDeGauss.

Solution

11
[23]: R=np.array([[1, 1, -1],[0.5, -0.25, 0],[0, 0.25, 0.5]])
U=np.array([[0],[5],[15]])
print(np.linalg.solve(R,U))
print(PivotDeGauss(R,U))

[[15.]
[10.]
[25.]]
[[15.]
[10.]
[25.]]

Références [1] Kiusalaas, J. (2013). Numerical methods in engineering with Python 3. Cam-
bridge university press.
[2] Numpy Package
[3] Mathplotlib Package
[4] Jupyter markdowns

12