Université de Picardie Jules Verne

Master 2 AAM
Calcul Scientifique
2019-2020

Méthodes Itératives 1 - Méthodes classiques

1 Commandes Scilab
• Partie triangulaire supérieure : triu(A)

• Partie triangulaire inférieure : tril(A)

• Calcul du temps de calcul

tic
.
.
.
instructions
.
.
.
toc

• Comparaison de l’historique des résidus

plot(it1,log(R1)/log(10))
plot(it2,log(R2)/log(10))
drawnow

• Matrice creuse
Convertir une matrice pleine en matrice creuse
A=sparse(X)
Convertir une matrice creuse en matrice pleine
A=full(X)

• Opérations sur les matrices creuses

– sprand : matrice creuse aléatroire

– speye : matrice identité creuse
– spones : matrice de 1 creuse
– spzeros : matrice nulle
– nnz : nombre déléments non nuls
– lufact: factorisation LU d’une matrice creuse
– lusolve : solveur LU de système linéaire creux
– spchol: factorisation de Cholesky d’une matrice creuse
– chsolve : solveur de Cholesky de système linéaire creux

1
 – spget : entrées non nulles d’une matrice creuse

Exercice 1
Effectuer et chronométrer le produit A ∗ B de deux matrices creuses, déclarées comme pleines, puis comme
creuses. Que constatez-vous ?

2
2 Méthodes classiques
On rappelle que les méthodes itératives classiques peuvent se programmer comme suit :

Algorithm 1 Méthode itérative classique

1: x(0) donné dans Rn
2: k = 0, r(0) = b − Ax(0)
3: while k < Kmax & kr(k) k >  do
4: Résoudre : M z (k) = r(k)
5: Poser : x(k+1) = x(k) + z (k)
6: Poser : r(k+1) = r(k) − Az (k)
7: Poser : k = k + 1
8: end while

Exercice 2 (Programmation générale)

1. Construire un seul programme principal Scilab pour les méthodes itératives classiques.

2. Soit la matrice  
2 −1
 −1 2 −1 
1  .

A= 2
 .. ... ... .

h  
 −1 2 −1 
−1 2
Résoudre le système linéaire Ax = b avec b=ones(n,1) avec les méthodes de Jacobi, Gauss-Seidel et
relaxation ; on partira à chaque fois de la donné initiae x=zeros(n,1). On rangera dans un vecteur
la suite des normes des résidus et on calculera le temps de calcul à chaque fois.

3. Comparer l’historique des résidus pour ces méthodes.

4. Conclusion ?
Exercice 3 (Méthode de Jacobi pour les systèmes triangulaires)
1. Soit A = D−F n×n, triangulaire supérieure, inversible, on note D = diag(A). Montrer que J = D−1 F
est nilpotente.

2. En déduire que la méthode de Jacobi converge en au plus n itérations.

3. Illustration sur un exemple de votre choix. On représentera graphiquement l’historique des résidus.

Exercice 4 (Méthode SSOR (Symmetric Successive Over Relaxation))

Soit A = D + L + LT une matrice symétrique. On considère la matrice M :

M = (D + L)D−1 (D + L)T

et plus généralement
   T
ω 1 −1 1
M (ω) = D + L (D) D+L
2−ω ω ω

3
 1. Programmer cette méthode pour différentes valeurs de ω ∈ [1, 2[. On calculera le nombre d’it́erations
nécessaires à la convergence.

2. Que constatez-vous ?

Exercice 5 (Directions alternées)N N N

1. On rappelle la propriété (A B)(C D) = AC BD. N N N
Soit A une admettant une décomposition LU . Montrer que (Id A) = (Id L)(Id U ).
N N
2. Soit A une matrice SDP et M = Id A + A Id = K1 + K2 . On considère la méthode itérative

Algorithm 2 Méthode des directions alternées

1: x(0) donné dans Rn
2: k = 0, r(0) = b − Ax(0)
3: while k < Kmax & kr(k) k >  do
4: Résoudre : (Id + α
2 K1 )u
(k+1/2) = (Id − α K )u(k) + α b
2 2 2
α
Résoudre : (Id + 2 K2 )u(k+1) = (Id − 2 K1 )u(k+1/2) + α
α
5:
2b
6: (k)
Poser : r = b − M u (k+1)

7: Poser : k = k + 1
8: end while

Ici α > 0. Exprimer u(k+1) en fonction de u(k) sous la forme u(k+1) = Gu(k) + c.

3. En déduire une CNS de convergence et montrer que si u(k) converge, c’est vers la solution de M u = b.

4. Montrer que K1 et K2 commutent. Montrer que ρ(G) < 1 et conclure.

5. Programmation de la méthode en Scilab

(a) Pour α > 0 fixé, calculer les décompositions de Cholesky de Id + α α

2 K1 et de Id + 2 K2 .
(b) Programmer la méthode des directions alternée pour la matrice obtenue avec
A=2*eye(N,N)-diag(ones(N-1,1),1)-diag(ones(N-1,1),-1)

Exercice 6 (Méthode de Multi-décompositions)

Soit A une matrice inversible. On considère m couples de matrices (Mi , Ni ) avec Mi inversible telles que

A = Mi − Ni , i = 1, · · · , m.

Soit alors m matrices diagonales à coefficients dans [0, 1] Ei telles que

m
X
Ei = Id
i=1

où Id désigne la matrice identité.

On introduit maintenant la méthode itérative

m
X
u(k+1)
= Ei Mi−1 (Ni u(k) + b) = Su(k) + T b
i=1

4
 1. On se place dans le cas m = 1. Expliciter la méthode et donner une CNS de convergence.

2. Maintenant m > 1.

(a) Expliciter S et T .
(b) Montrer que
Id − T A = S.

(c) En déduire une CNS de convergence.

3. On se place dans le cas m = 2 et pour simplifier on impose

E1 = αId, E2 = (1 − α)Id

(a) Montrer que l’algorithme peut s’écrire sous la forme

u(k+1) = αv (k+1) + (1 − α)w(k+)

où v (k+1) et w(k+1) sont à préciser.

(k+1) (k+1)
(b) Exprimer le résidu r(k+1) = b−Au(k+1) en fonction de rv = b−Av (k+1) et de rw b−Aw(k+1) .
Quelle valeur de α permettrait d’obtenir une convergence plus rapide ?
(c) On modifie maintenant l’algorithme précédent en changeant α = αk à chaque itération. Calculer
αk de sorte à minimiser kr(k+1) k2 .
(d) Montrer qu’avec ce choix

kr(k+1) k2 ≤ krv(k+1) k2 et kr(k+1) k2 ≤ krw

(k+1)
k2 .

(e) Conclure.

Remarque : l’intérêt de cette méthode provient de ce que les calculs Ei Mi−1 (Ni u(k) + b) peuvent être
exécutés indépendamment les uns des autres. Cela se réalise avantageusement sur des ordinateurs ”par-
allèles”, c’est à dire avec plusieurs processeurs, chacun ayant un calcul à executer ”en même temps”
indépendemment des autres. Le résultat final est assemblé une fois tous les calculs terminés.

5
3 Méthodes particulières
Les Méthodes suivantes reposent sur l’observation que la solution d’un système linéaire m × n Ax = b , si
elle existe, est dans l’intersection des hyperplans affines :

Hi : ai,1 x1 + ai,2 x2 + · · · ai,n xn = bi , i = 1, · · · , m.

Bien sûr l’intersection est réduite à un point lorsque A est inversible. Une idée naturelle consiste donc à
construire une suite de projetée orthogonales successives sur ces hyperplans affines. Il n’est pas difficile de
voir que
P rojHi (x) = Pi (x)y
avec < ai , y >= bi et x − y = αai . Du coup

< x, ai >=< y, ai > +αkai k22 = bi + αkai k22 .

< x, ai > −bi b − < x, ai >

Il en découle que α = 2 et que Pi (x) = x + i ai .
kai k2 kai k22

Exercice 7 (Méthode de Kaczmarz)

Soit A une matrice m × n. Soit b ∈ Rm . On considère le système linéaire

Ax = b

On désigne par ai la ième ligne de A.

Algorithm 3 Méthode de Kaczmarz

1: x(0) donné dans Rn
2: while k < Kmax & kr (k) k >  do
bi − hai , xk i
3: Poser : xk+1 = xk + ai
kai k2
4: (avec i = k mod m, i = 1, 2, ..., m)
5: Poser : k = k + 1
6: end while

Cette méthode s’interprète de la manière suivante : xk+1 est la projection orthogonale de xk dans l’hyperplan
< ai , x >= bi , soit
xk+1 = Pi (xk ).
Le taux de convergence est 1 − 12 , pour κ = kAkkA−1 k assez grand.
κ
On peut généraliser cette méthode en introduisant un paramètre de relaxation λk comme suit :

bi − hai , xk i
xk+1 = xk + λk ai
kai k2

1. Programmer la méthode pour A m × n de rang maximal (rang(A) = m ≤ n).

2. Comparer la solution obtenue avec celle au sens des moindres carrés.

6
Exercice 8 (Méthode de Cimino)
La méthode de Karzmarz utilise les lignes de A séquentiellement, c’est à dire les unes après les autres. L a
méthde de Cimmino quant à elle utilise les lignes de A simultanément et calcule l’itération courante comme
moyenne de toutes les projection de l’itérée précédente, soit
m m m
1X 1X 1 X bi − hai , xk i
xk+1 = Pi (xk ) = xk + Pi (xk ) = xk + ai
m m m kai k2
i=1 i=1 i=1

1. Programmer la méthode pour A m × n de rang maximal (rang(A) = m ≤ n).

2. Comparer la solution obtenue avec celle au sens des moindres carrés.

Remarque 1
La méthode de Cimino peut se réécrire comme suit :
m
1
X bi − hai , xk i
xk+1 = xk + m Pi (xk ) ai
kai k2
i=1
b1 − < a1 , xk >
 
 
1
= xk + m a1 , · · · , am  .. 
ka1 k22 kam k22
 . 
bm − < am , x >k

 T
 1 
   
a1 ka k2 a1
 1 2 
1 
= xk + m .   .. ..  xk 
 ..    b − 
 
. .  
1
 
am 2
am
kam k2
= xk + AT M −1 b − Axk