Transcript ASI 3 Méthodes numériques pour l’ingénieur Résolution de systèmes linéaires par des méthodes itératives : Jacobi, Gauss Seidel, relaxation.

ASI 3
Méthodes numériques
pour l’ingénieur
Résolution de systèmes linéaires
par des méthodes itératives :
Jacobi, Gauss Seidel, relaxation
Résoudre un système linéaire
 6
2 x1  4 x2  2 x3
 x  3x
 x4  0
 1
2
Ax  b  
3 x1  x2  x3  2 x4  8

 x2  2 x3  x4  6
soit  ~
x1 , ~
x2 , ~
x3 , ~
x4  une " presque" solution
~
~
~
imaginons que A~
x b
avec par exemple dist  A, A   
 6  4~
x2  2 ~
x3

x

 1
2

0~
x1  ~
x4
x

 2
3
si on essayait : 
~
8  3 x1  ~
x2  2 ~
x4
 x3 
1

~
~
 x  6  x2  2 x3
 4
1
n
bi   aij ~
xj
soit
xi 
j 1, j i
aii
Résoudre un système linéaire en itérant
 6  4~
x2  2 ~
x3

x

 1
2
~

0  x1  ~
x4
 x2 
3

~
8  3 x1  ~
x2  2 ~
x4
 x3 
1

~
~
 x  6  x2  2 x3
 4
1
n
bi   aij ~
xj
soit
xi 
j 1, j  i
aii
Si Ax n’est pas encore égale à b, on recommence !
tant que dist ( Ax new , b)  
n
(i.e.10-12 )
bi   aij x old
j
xinew 
fin
j 1, j i
aii
Osons itérer ! méthode de Jacobi
tant que dist ( Ax new , b)  
(i.e.10-12 )
n
bi   aij x old
j
xinew 
j 1, j i
aii
fin
Soit D la diagonale de la matrice A, et G le reste :
A = D+G
0 
 a11 0 0 0


0 
 0  0 0
D   0 0 aii 0
0 


 0 0 0  0 
 0 0 0 0 a 

nn 
a12
a1i

a1n 
 0


0
ai 1,i

 
 a21
0
ai 1,i
ai ,n 
G   ai1 ai ,i 1


 ai ,i 1
0
an 1,n 
 
a
ani an,n1
0 
 n1 

Dx new  b  Gx old  x new  D 1 b  Gx old

Gauss Seidel
2 x1  4 x2  2 x3
 x  3x
 x4
 1
2
Ax  b  
3 x1  x2  x3  2 x4

 x2  2 x3  x4
~
~
 6  4 x2  2 x3

x

 1
2
i 1
n

~
0  x1  ~
x4
b

a
x

a


i
ij j
ij x j
 x2 
j 1
j i 1
3
soit
x


i
8  3 x1  x2  2 ~
x4
aii
 x3 
1

 x  6  x2  2 x3
 4
1
 6
 0
 8
 6
méthode de Gauss-Seidel
tant que dist ( Ax new , b)  
i 1
bi  
xinew 
j 1
(i.e.10-12 )
aij x new
j
n
  aij x old
j
j i 1
aii
fin
Soit E la triangulaire inférieure et F la supérieure de la matrice A :
A = D+E+F
0
0

 0

0
0

 a21
0
0
E   ai1 ai ,i 1

 ai ,i 1
0
 
a
ani an,n1
 n1 
0
 0 a12



0 0
0 F  0 0


0
 
0 
0 

a1i

ai 1,i
0

ai 1,i
0
0
0
0
 D  E x new  b  Fx old  x new   D  E 1 b  Fx old 
a1n 

 
ai ,n 

an1,n 
0 
La relaxation
 
x new   x old
xrnew  x new  (1   ) x old
  0 : paramètre de relaxation
Jacobi :
1
x new   D  b  ( E  F ) x old


   0,1 : interpollation

  1 : status quo
  1,2 : extrapollation



xrnew    D 1 b  Gx old  1   x old
en multipliant par D
Dxrnew  b  1   D  G x old
Gauss  Seidel :
x new   D  E 1 b  Fx old




xrnew    D  E  b  Fx old  1   x old
1
 D  E xrnew  b  1   D  F x old
Résumé « algorithmique »
Mx new  Nx old  b
Jacobi :
Dx
new

 b  (E  F )x
Gauss  Seidel :
 D  E x
new

 b  Fx
old
old


 M D
:
N  E  F
M  D  E
:
 N  F
Relaxation :


new
D

E
x
b

1

1


D  F x old
 M  1 D  E

:
1
 N   D  F
Convergence
Principes généraux
 x 0 donné,
 k 1
k
 x  Cx  d
Théorème : S' il existe une norme matricielle subordonné e telle que
C 1
Alors l' algorithme ci  dessus converge pour tout x 0
vers x* la solution de :
 I  C x*  d
Éléments de démonstration
- x* est un point fixe de l’algorithme
k
k
-e  x  x*
 Cx k 1  d  Cx * d
 e k  Ce k 1  C k e 0
 C x k 1  x *
k
k
e k  C k e0  C e0
si C  1, C k
 0



Normes matricielles
 n( x)  0 positivité
 n( x )  0  x  0

vérif iant 
 n(x)   n( x)
n( x  y )  n( x)  n( y )
 n : E  R
norme : 
 x  n( x )
exemples E  R n ;
n
x 
2
2
i 1
xi2
; x
p
p
n
n
  xi ( p  1) ; x 1   xi ; x   sup xi
i 1
p
i 1
i 1, n
Définition
Soit A une matrice nxm, étant donnée une norme vectorielle,
on appelle norme matricielle subordonnée, la norme matricielle
définie par :
Ax
A  max
xR n , x  0 x
Conséquence :
~
x est ce max : A~
x  A ~
x
Exemples
norme de Frobénius : A
n
2
F
m
n
   aij2  tr  A'A ;
j 1 i 1
m
A 1  max  aij ;
A   max  aij ;
1 j  m i 1
1i  n j 1
A 1  A' 
à utiliser pour le calcul !
Illustration 2d
x 1
2
A  max
n
Ax
xR , x 1
x2
x2
2
A
2
Ax
1.5
1
Ax
0.5
x
0
x1
x1
-0.5
-1
-1.5
-2
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Calculez les valeurs propres de
1 0 2


 0 1  1
 1 1 1 


1  1, 2  1  i 3, 3  1  i 3,
Et ses vecteurs propres ?
2i 3 
2i 3 



1
 3 
 3 
 
v1   1  , v2   i 3  , v3    i 3 
 3 
 3 
0
 1 
 1 
 




Rayon spectrale d’une matrice
Définition : on appelle rayon spectrale d’une matrice carrée A,
le nombre réel r(A) tel que :
r  A  max i ( A)
i 1,n
Théorème : soit A une matrice nxm, alors :
2
A 2  r  A' A
Corollaire : si A est une matrice carrée symétrique nxn, alors :
A 2  r  A
Remarque : en général, le rayon spectrale n’est pas une norme :
0 1
A
 , A 2  1  0, r  A  0
0 0
Convergence : le retour
 x donné,
 k 1
k
 x  Cx  d
0
Principes généraux
Théorème : S' il existe une norme matricielle subordonné e telle que
C 1
Alors l' algorithme ci  dessus converge pour tout x 0
vers x* la solution de :
 I  C x*  d
Théorème : les points suivants sont équivalents :
• C est une matrice convergente, (i.e. Ck tend vers 0)
• r(C)<1
• C 1
Résumé « algorithmique »
Mx new  Nx old  b
Jacobi :
Dx
new

 b  (E  F )x
Gauss  Seidel :
 D  E x
new

 b  Fx
old
old


 M D
:
N  E  F
M  D  E
:
 N  F
Relaxation :


new
D

E
x
b

1

1


D  F x old
 M  1 D  E

:
1
 N   D  F
Convergence
Théorème : Si A est une matrice à diagonale strictement dominante,
alors la méthode de Jacobi converge.
Démonstration Jacobi :

x k 1   D  b  ( E  F ) x k
1
x k 1  Cx k  d
C
avec C  D 1 ( E  F ) , d  D 1b
1


 D (E  F )

 1 n
 max 
aij

i 1, n aii j 1


  1

Remarque : il en est de même pour la méthode de Gauss-Seidel
Convergence
Théorème : soit une méthode itérative : Mxk 1  Nx k  b
Si A est une matrice symétrique définie positive
telle que si A = M-N alors M+N’ est définie positive
Alors la méthode itérative est convergente
Démonstration A symétrique définie positive  x 2A  x' Ax
2
2
2
1
1
1
M Nx  M M  Ax  x  M Ax
A
A
A
posons y  M 1 Ax
1
x  M Ax
2
A
 My  Ax
 x y
2
A
  x  y ' A x  y ' A  x A  x' Ay  y ' Ax  y ' Ay
2
 x A  y ' M ' y  y ' My  y ' Ay  x A  y ' M  N ' y
2
1
on a donc : x  M Ax
2
2
A
 x
2
A
 M 1 N
A
 max M 1 Nx
x , x A 1
A
1
Théorème : Si A est une matrice symétrique définie positive,
la méthode de la relaxation converge pour : 0    2
Influence de w
rayon spectral de la matrice M-1N
rayon spectral
1
0

2
Remarques
pratique :
• pas de preuve de convergence généralisée,
• on préfère la relaxation avec différents tests pour ,
• on préfère les méthodes directes,
• voir les méthodes semi directes pour les problèmes de
grande taille (cf les méthodes « multigrilles »),
Conditionnement d’un système linéaire
2  x1   3 
 1
Ax  b  
    
1.0001 2  x2   3 
examinons deux solutions
 x1  1
x  
 x2  1
 y1   3 
y    
 y2   0 
rx  ry  0.0002 et
 rx  Ax  b  0 0T

T
ry  Ay  b  0 0.0002
x y 2
Deux vecteurs très différents donnent des solutions très proches
x2
1
1
3
x1
Conditionnement : influence du second membre
 8 10 10 10 
 38 
 1


 
 
1 5 6 1 
 13 
 1
A
, b   , admet comme solution x   

7 10 4 7
28
1


 
 
7 8 1 7 
 23 
 1
 38.1 
  34.5 




12.9 
11.0 
b


 3 x
b  b  
, x  x  
:
 3.7 10
 22.6


28.1
 5. 6
b
x




 22.9 
 26.0 
1  24.1, 4  0.007 : cond ( A)  7363
r  b  b  b 
 Ax  A~
x
avec ~
x  ( x  x )
 A x  ~
x   x  ~
x   A1r

x~
x  A1 r
A non singulière
Conditionnement
Définition : on appelle conditionnement d’une matrice carrée A, relatif
à une norme subordonnée, le nombre réel c(A) :
c  A  A
Remarque :
 AA 1  A
I
A1
A1  c  A  1
Théorème : Si A est une matrice carrée, non singulière (régulière)
Ax  b,
x
x
A x  x   b  b
 c  A
b
b
Perturbation du second membre
Ax  b,
 A  A x  x   b
x
A
 c  A
x  x
A
Perturbation de la matrice
Un problème est dit « bien conditionné » si c(A) est proche de 1,
il est dit « mal conditionné » si c(A) est grand (et mal posé si c(A) est infini)
Conditionnement
1
1
1) Ax  b  b  A x   A
x
b
2) A x  x   b  b  b  A x  x  A1 b
1)  2) 
x
x
 A
A
1
b
c  A  A
b
A1
Remarque : si A est symétrique, si on note ses valeurs propres i
1
1 , 2 ,..., i ,..., n ,
 A1 n  A
1
2
2
n
et donc c 2  A 
1
Dans l ’exemple, c 2  A  7363, b  0.0037
la borne de l' erreur est de l' ordre de c 2  A b  27.2 (on a trouvé 22.5)
Comment améliorer le
conditionnement ?
Ajouter un « chouia » sur la diagonale
n  
c 2  A  I  
1  
Itérations !
Ax  b
~
x  gauss( A, b)
err  A~
x b
A ~
x  x   b
~
Ax  b  b  A~
x
A = randn(n);
b= ones(n,1);
x = A\b;
err = A*x-b;
norm(err)
ans = 2.8246e-013
dx = A\err;
err2 = A*(x-dx)-b;
norm(err2)
ans = 6.4789e-014
TP - la relaxation
Le but du TP est d’écrire un programme matlab résolvant
un système linéaire par la méthode de la relaxation
x = relax(A,b,w,nite,err)
Pour ce faire, il faut étudier l’évolution du rayon spectal
r M 1 N en fonction de 
- mettez vous par binôme
- rédigez une page :
recto : ce que vous avez fait
verso ce que vous en pensez
- a rendre pour le 8 décembre à 17h30 (publication du corrigé)


Indices : créer un problème test,
les fonction cputime et flops
tril et triu pourraient vous simplifier la vie
et diag(diag()) et eig aussi
Propriétés
Définition : on appelle quotient de Rayleigh la fonction qA(x)
x ' Ax
x' x
Soit A une matrice carrée, on appelle
polynôme caractéristique de A le polynôme défini par :
p( )  det A  I 
q A ( x) 
Les n racines i de ce polynôme sont les valeurs propres de A,
vi est un vecteur propre de A. Il existe n vecteurs vi tels que :


Avi  i vi
i est une valeur propre de A,
vi est un vecteur propre de A.
Théorème : si A est symétrique,
max i  max q A ( x)
i
x