Solveurs Linéaires pour les problèmes industriels

TP 1 : Résolution de l’équation de Poisson

On s’intéresse à l’équation de la chaleur stationnaire (connue aussi comme l’équation de
Poisson) :
−∆u = f
avec des conditions aux limites de type Dirichlet. f représente la source de chaleur alors que
u représente la température dans la pièce. On considère une discrétisation par différences
finies sur un domaine carré [0, 1]2 :
4ui,j − ui+1,j − ui−1,j − ui,j+1 − ui,j−1
= fi,j
∆x2
où ui,j est la température au point (xi , yj ) = (i∆x, j∆x). avec

1
∆x =
n+1
où n est le nombre de points suivant x ou y. On obtient ainsi le système linéaire

AU = F

avec la matrice A tridiagonale par blocs

 
H −I ··· 0
0
. .
 −I H −I . . ..
 


A = (n + 1)2  .. .. .. 
 0 . . . 0 

 . .
 .. . . −I H −I


0 · · · 0 −I H

avec I la matrice identité de taille n et H la matrice tridiagonale

 
4 −1 0
 −1 4 −1 
 
H=
 ... ... ... 

 
 −1 4 −1 
0 −1 4

Les développements du TP se feront en Python. on pourra utiliser Visual Studio ou le duo

emacs/ipython. Pour cette dernière possibilité, emacs sert à éditer les fichiers sources alors
que ipython vous donne une console pour tester les fonctions. Cette console s’ouvre en tapant
sur un terminal :

1
 conda_pytorch
ipython --pylab
L’appel à conda pytorch est nécessaire pour les machines de l’enseirb. Ce TP est évalué, il
vous faudra déposer un rapport ainsi que les sources sur Moodle au plus tard le jeudi 28
novembre. Vous pouvez réaliser ce travail seul ou à deux personnes. Il faudra indiquer les
noms des participants sur le rapport.

Question 1 Ecrire une fonction générant la matrice A pour un stockage plein. On rappelle
qu’en python, on peut utiliser la fonction diag, par exemple
1 B = diag (2* ones ( n )) - diag ( ones (n -1) , 1) - diag ( ones (n -1) , -1)

construit la matrice  
2 −1 0
 −1 2 −1 
 
B=
 .. .. .. 
. . . 
 
 −1 2 −1 
0 −1 2
Il est à noter que la matrice A est une matrice de taille n2 avec cinq diagonales non nulles.
Faites attention que la sous-diagonale de A est le vecteur suivant

(n + 1)2 × (−1, −1, · · · , −1, 0, −1, · · · , −1, 0, · · · )

On a (n-1) coefficients égaux à -1 suivis d’un 0, puis (n-1) coefficients à -1, etc. En python,
on peut sélectionner certains indices avec les deux points
1 x = rand (10)
2 x [0:4] = 1.0; # 0:4 prend les indices 0 , 1 , 2 et 3
3 x [6:] = -2.0; # 6: prend les indices 6 , 7 , etc
4 # jusqu ’ au dernier element du vecteur
5 x [2:7:2] = 3.0; # 2:7:2 : on part de l ’ indice 2
6 # jusqu ’a 7 par pas de 2 (7 etant exclu )
7 x [3::2] = 4.0; # on part de l ’ indice 3 jusqu ’a la fin
8 # par pas de 2 ( donc 3 , 5 , 7 , etc )

Vérifier que vous avez la bonne matrice pour des petites valeurs de n (3 et 4 par exemple).

Question 2 Construire un second membre qui correspond à une source de chaleur au

centre. Une solution est de prendre un vecteur F de taille n2 rempli de 0 excepté 1 à l’indice
n2 n
+ . Appelez ensuite solve (du package numpy.linalg) suivi de imshow pour afficher la
2 2
solution
1 u = solve (A , F )
2 imshow ( reshape (x , (n , n )) , cmap = cm . hot )

2
 Vous devez obtenir une image similaire à la figure ci-dessous (pour n = 50)

Question 3 Comme vous l’avez noté, la matrice A est particulièrement creuse. on va

maintenant utiliser des fonctions du package scipy.sparse afin de construire directement une
matrice creuse et utiliser des fonctions des matrices creuses afin d’avoir des temps de calcul
réduits. Ecrire une fonction générant la matrice A pour un stockage creux. On pourra utiliser
la fonction spdiags. Voici un exemple d’utilisation de cette fonction
1 data = zeros ([3 , 5]) # syntaxe : zeros ([ nb_diag , n ])
2 data [0 , :] = 4.0* ones (5) # diagonale principale
3 data [1 , :] = array ([ -2.0 , -3.0 , 1.0 , -1.0 , 0.0]) # sous - diagonale
4 data [2 , :] = array ([0.0 , -2.0 , -3.0 , 1.0 , -1.0]) # sur - diagonale
5 B = spdiags ( data , [0 , -1 , 1] , 5 , 5). tocsr ()
6 # derniers arguments : nombre de lignes et colonnes

qui construit la matrice  

4 −2 0 0 0

 −2 4 −3 0 0 

B=
 0 −3 4 1 0 

 0 0 1 4 −1 
0 0 0 −1 4
L’appel à la methode tocsr() permet de convertir la matrice au format CSR (Compressed
Sparse Row). On utilisera plus loin des fonctions qui ne marcheront que pour ce type de
stockage.

Question 4 Calculer la solution en appelant la fonction spsolve (au lieu de solve pour les
matrices pleines). Observez le gain en temps de calcul en utilisant une matrice creuse.

Question 5 Implémenter l’algorithme du gradient conjugué (en prenant en argument la

matrice creuse). Pour réaliser le produit matrice vecteur, il faut appeler la méthode dot :
1 y = A . dot ( x ) # produit matrice vecteur y = A x

3
 La fonction dot est à utiliser pour le produit scalaire entre deux vecteurs
1 alpha = dot (x , y ) # produit scalaire alpha = <x , y >

On rappelle l’algorithme du gradient conjugué précondtionné

if ||b|| = 0 then
Retourner x = 0
end if
x = x0
r = b − Ax0
k=0
while ||r||/||b|| > ε et k ≤ Nmax do
Appliquer preconditionneur z = M −1 r
ρ = < r, z >
if ρ == 0 then
Echec de l’algo
end if
if k == 0 then
p=z
else
γ = ρ/ρo
p = γp + z
end if
Produit matrice-vecteur : q = Ap
δ = < p, q >
if δ == 0 then
Echec de l’algo
end if
α = ρ/δ
x = x + αp
r = r − αq
ρo = ρ
k =k+1
end while
Pour le préconditionneur, on vous propose de faire une classe contenant la méthode apply
et de passer la classe en argument de la fonction implémentant le gradient conjugué. On
implémentera dans un premier temps le préconditionneur identité :
1 # preconditioneur identite
2 class id_precond :
3 # fonction principale appliquant le preconditionneur
4 def apply ( self , x ):
5 return x
6

4
 7 # constructeur , mettre ici les calculs prealables du preconditionneur
8 # de telle sorte qu ’ apply soit le plus rapide possible
9 def __init__ ( self , A ):
10 # si on a besoin par exemple de la diagonale de A
11 # les variables precedees de self . sont des attributs de la classe ( acce
12 self . D = A . diagonal ();
13

14 # pour initialiser un objet de la classe :

15 # ici on donne la matrice en argument du constructeur
16 prec_id = id_precond ( A )
17

18 # ensuite lancer le gradient conjugue avec cet objet

19 epsilon = 1e -6; Nmax = 1000
20 x , residu = c on ju ga te _g ra di en t (A , b , x0 , prec_id , epsilon , Nmax )
21 # dans co nj ug at e_ gr ad ie nt
22 # la ligne z = M ^{ -1} r s ’ implementera comme z = prec . apply ( r )

L’avantage de faire une classe est que vous pourrez stocker des attributs dedans (utiles
pour les autres préconditionneurs). La méthode implémentée renverra la solution ainsi que
l’historique des résidus relatifs (coefficients ||r||/||b||). On pourra afficher dans un graphe le
résidu en fonction du nombre d’itérations pour vérifier la convergence de la méthode.

Question 6 Implémenter le précondionneur Gauss-Seidel symétrique vu en cours. On rap-

pelle qu’il s’écrit :

 Résoudre le système triangulaire (D-E) q = r
Résoudre M z = r ⇔ Multiplier par la diagonale p = D q
Résoudre le système triangulaire (D-F) z = p


En python, on peut récupérer la partie triangulaire inférieure (ou supérieure) de A en écrivant

1 L = tril ( A ). tocsr ();
2 U = triu ( A ). tocsr ();

Les fonctions tril et triu font partie du package scipy.sparse.

Pour la résolution avec un système triangulaires, la fonction spsolve triangular fournie par
scipy est très lente. On vous propose de télécharger deux fichiers en tapant sur un terminal :
wget https://www.math.u-bordeaux.fr/~durufle/module_solve.py
wget https://www.math.u-bordeaux.fr/~durufle/mod_solve_triangular.f90
On compilera le fichier mod solve triangular.f90 avec f2py :
f2py -c mod_solve_triangular.f90 -m mod_solve_triangular
C’est une commande à taper dans un terminal à part, il ne faut pas la taper dans Visual
Studio par exemple. Bien sûr, il faut avoir placé le fichier mod solve triangular.f90 dans
le même répertoire que les autres fichiers Python que vous avez créés. Dans votre fichier
Python, il suffira d’inclure le module pour l’utiliser :

5
1 from module_solve import *
2 x = solve_tri (L , b , lower = True ); # pour resoudre L x = b
3 # en supposant L triangulaire inferieure
4

5 x = solve_tri (U , b , lower = False ); # pour resoudre U x = b

6 # en supposant U triangulaire superieure

La fonction solve tri n’est implémentée que pour des matrices stockées au format CSR.
Si vous n’arrivez pas à utiliser f2py sur votre machine, vous pouvez utiliser la fonction
spsolve triangular ou écrire une fonction Python adaptée à des matrices tridiagonales. Si
vous utilisez spsolve triangular, il faudra plutôt comparer les complexités (calculées à partir
du nombre d’itérations) et ne pas prêter attention aux temps de calcul Python.
Tester le gradient conjugué préconditionné par Gauss-Seidel symétrique.

Question 7 La factorisation de Cholesky incomplète consiste à calculer les coefficients de

la matrice L uniquement pour les éléments non-nuls de A. On vous propose d’implémenter
cette factorisation pour le cas particulier de la matrice A (matrice de taille n2 ), l’algorithme
de factorisation s’écrit
L=0
for i = 1, n*n doq
Calculer Li,i = ai,i − L2i,i−1 − L2i,i−n
ai+1,i
Calculer Li+1,i =
Li,i
ai+n,i
Calculer Li+n,i =
Li,i
end for
On fera attention aux indices pour éviter de dépasser les limites des tableaux. Par exemple,
pour le calcul de Li,i , on fera attention que le coefficient Li,i−1 n’existe que si i est plus
grand que 1, et le coefficient Li,i−n n’existe que si i est plus grand que n. Implémenter la
factorisation de Cholesky incomplète pour la matrice A. On retournera la matrice L sous la
forme d’une matrice creuse (on peut calculer les diagonales et utiliser spdiags).

Question 7 On veut utiliser la factorisation de Cholesky incomplète pour préconditionner

le système. L’application du préconditionneur M −1 revient à résoudre deux systèmes trian-
gulaires :
z = M −1 r ⇔ Résoudre LLT z = r
(
Résoudre Ly = r
⇔
Résoudre LT x = y
Pour la transposée de L, on pourra utiliser L.T.tocsr() pour avoir une matrice creuse
de type CSR. Implémenter la classe applicant le préconditionneur de Cholesky incomplet.

6
Tester l’algorithme du gradient conjugué avec ce préconditionneur. Comparer les différents
préconditionneurs (identité, Gauss-Seidel symétrique et Cholesky incomplet).