142
4 M´ ethodes it´ eratives pour la r´ esolution des syst` emes lin´ eaires
Exemple 4.5 R´ esolvons par la m´ ethode du gradient le syst` eme lin´ eaire associ´ e ` a la matrice Am ∈ R
m×m construite dans MATLAB par les commandes
G=numgrid(’S’,n); A=delsq(G) avec m = (n − 2)
2 . Cette matrice provient de la
discr´ etisation de l’op´ erateur de Laplace sur le domaine [−1, 1]
2 . Le second membre
bm est tel que la solution exacte soit le vecteur 1
T ∈ R
m . La matrice Am est sym´ etrique d´ efinie positive pour tout m et elle est mal conditionn´ ee pour m grand. On
ex´ ecute le Programme 19 pour m = 16 et m = 400, avec x
(0) = 0, tol=10
−10 et
nmax=200. Si m = 400, la m´ ethode ne parvient pas ` a satisfaire le crit` ere d’arrˆ et dans
le nombre d’it´ erations imparti et le r´ esidu d´ ecroˆ ıt tr` es lentement (voir Figure 4.6).
Dans ce cas en effet K2(A400) 258. En revanche, si on pr´ econditionne le syst` eme
avec la matrice P = R
T
in Rin, o` u Rin est la matrice triangulaire inf´ erieure de la factorisation de Cholesky incompl` ete de A, l’algorithme satisfait le test d’arrˆ et dans le
nombre d’it´ erations fix´ e (` a pr´ esent K2(P
−1 A400) 38).
•
0
50
100
150
200
250
10
−14
10
−12
10
−10
10
−8
10
−6
10
−4
10
−2
10
0
(a)
(b)
(c)
(d)
Fig. 4.6. R´ esidu normalis´ e par le r´ esidu initial en fonction du nombre d’it´ erations
pour la m´ ethode du gradient appliqu´ ee aux syst` emes de l’Exemple 4.5. Les courbes
(a) et (b) concernent le cas m = 16 sans et avec pr´ econditionnement, les courbes
(c) et (d) concernent le cas m = 400 sans et avec pr´ econditionnement
4.3.4 La m´ ethode du gradient conjugu´ e
La m´ ethode du gradient consiste essentiellement en deux phases : choisir une
direction de descente (celle du r´ esidu) et un point o` u Φ atteint un minimum
local le long de cette direction. La seconde phase est ind´ ependante de la premi` ere. En effet, ´ etant donn´ e une direction p
(k) , on peut choisir α k comme
´ etant la valeur du param` etre α telle que Φ(x
(k) + αp
(k) ) est minimum. En
´ ecrivant que la d´ eriv´ ee par rapport `
a α s’annule au point o` u la fonction admet
un minimum local, on obtient
α k =
p
(k) T r
(k)
p (k) T Ap (k)
,
(4.38)
(ce qui redonne (4.36) quand p
(k) = r
(k) ). On peut se demander si un autre
choix de direction de recherche p
(k) ne pourrait pas donner une convergence
4 M´ ethodes it´ eratives pour la r´ esolution des syst` emes lin´ eaires
Exemple 4.5 R´ esolvons par la m´ ethode du gradient le syst` eme lin´ eaire associ´ e ` a la matrice Am ∈ R
m×m construite dans MATLAB par les commandes
G=numgrid(’S’,n); A=delsq(G) avec m = (n − 2)
2 . Cette matrice provient de la
discr´ etisation de l’op´ erateur de Laplace sur le domaine [−1, 1]
2 . Le second membre
bm est tel que la solution exacte soit le vecteur 1
T ∈ R
m . La matrice Am est sym´ etrique d´ efinie positive pour tout m et elle est mal conditionn´ ee pour m grand. On
ex´ ecute le Programme 19 pour m = 16 et m = 400, avec x
(0) = 0, tol=10
−10 et
nmax=200. Si m = 400, la m´ ethode ne parvient pas ` a satisfaire le crit` ere d’arrˆ et dans
le nombre d’it´ erations imparti et le r´ esidu d´ ecroˆ ıt tr` es lentement (voir Figure 4.6).
Dans ce cas en effet K2(A400) 258. En revanche, si on pr´ econditionne le syst` eme
avec la matrice P = R
T
in Rin, o` u Rin est la matrice triangulaire inf´ erieure de la factorisation de Cholesky incompl` ete de A, l’algorithme satisfait le test d’arrˆ et dans le
nombre d’it´ erations fix´ e (` a pr´ esent K2(P
−1 A400) 38).
•
0
50
100
150
200
250
10
−14
10
−12
10
−10
10
−8
10
−6
10
−4
10
−2
10
0
(a)
(b)
(c)
(d)
Fig. 4.6. R´ esidu normalis´ e par le r´ esidu initial en fonction du nombre d’it´ erations
pour la m´ ethode du gradient appliqu´ ee aux syst` emes de l’Exemple 4.5. Les courbes
(a) et (b) concernent le cas m = 16 sans et avec pr´ econditionnement, les courbes
(c) et (d) concernent le cas m = 400 sans et avec pr´ econditionnement
4.3.4 La m´ ethode du gradient conjugu´ e
La m´ ethode du gradient consiste essentiellement en deux phases : choisir une
direction de descente (celle du r´ esidu) et un point o` u Φ atteint un minimum
local le long de cette direction. La seconde phase est ind´ ependante de la premi` ere. En effet, ´ etant donn´ e une direction p
(k) , on peut choisir α k comme
´ etant la valeur du param` etre α telle que Φ(x
(k) + αp
(k) ) est minimum. En
´ ecrivant que la d´ eriv´ ee par rapport `
a α s’annule au point o` u la fonction admet
un minimum local, on obtient
α k =
p
(k) T r
(k)
p (k) T Ap (k)
,
(4.38)
(ce qui redonne (4.36) quand p
(k) = r
(k) ). On peut se demander si un autre
choix de direction de recherche p
(k) ne pourrait pas donner une convergence