TD MATLAB

Stefanie Hahmann, ENSIMAG

1 Méthode QR, recherche de valeurs propres

Rappel sur la factorisation QR :
Soit A une matrice m × n, sa factorisation QR est une factorisation A = QR avec Q
matrice m × ·m orthogonale
¸ (ie Q ∗ Q0 = Id) et R matrice m × n triangulaire supérieure
R1
de la forme où R1 est une matrice triangulaire supérieure n × n et la matrice
0
nulle au-dessous est de taille n × m − n.
Soit A une matrice carrée n × n.

1. Algorithme QR pour les valeurs propres

Partant de A = A0 , on crée une suite de matrices A1 , A2 , . . . , Ar , Ar+1 ainsi définie :

Méthode I
On factorise A en QR : Ar = Qr Rr
on retourne la factorisation : Ar+1 = Rr Qr
(r = 0, 1, 2, . . .)
On sait que :
- toutes les matrices Ar sont semblables à A0 .
- si A0 est symétrique, toutes les Ar le sont aussi.
- si on pose :
Qr = Q0 Q1 . . . Qr
Rr = Rr Rr−1 . . . R0
alors on a :
Qr Ar+1 = A0 Qr
Qr Rr = Ar+1 0

2. Ecrire un programme générant la suite Ar . L’expérimenter sur une matrice A0

(10 × 10) symétrique. On observera l’évolution de Ar (r jusqu’à 100 ou 150 si
nécessaire) : Visualiser la convergence de Ar vers la diagonale des valeurs
propres de A0 .
Recommencer pour quelques matrices A0 symétriques 5 × 5 ou 10 × 10.

2 Résolution de systèmes linéaires Ax = b

1. Méthode LU
Il s’agit ici d’utiliser la commande de factorisation LU, puis de résoudre successive-
ment les deux systèmes triangulaires.
2. Méthode de Jacobi
C’est une méthode itérative, il est donc intéressant de noter le nombre d’itérations.
On décompose la matrice A en L + D + U où L et U sont respectivement les parties
triangulaires strictement inférieure et strictement supérieure, et D est la diagonale
de A. On définit les matrices MJ = D et NJ = −(L + U ).
x0 = 0
xk+1 = MJ−1 ∗ (b + NJ xk )
On peut aussi l’écrire à l’aide des composantes :

i−1
X n
X
xk+1
i = (bi − Ai,j xkj − Ai,j xkj )/Ai,i
j=1 j=i+1

Traiter quelques exemples en générant des systèmes à votre guise.

3. Gauss-Seidel
On utilise la décomposition A = L + D + U , on définit les matrices MG = L + D et
NG = −U . L’itération s’écrit (en prenant x0 = 0) :
−1
xk+1 = MG ∗ (b + NG xk )

c’est-à-dire
i−1
X n
X
xk+1
i = (bi − Ai,j xk+1
j − Ai,j xkj )/Ai,i
j=1 j=i+1

Traiter les mêmes exemples pour comparer avec Jacobi.

4. SOR (Successive Over-Relaxation)

Il s’avère que la méthode de Gauss-Seidel converge d’autant plus lentement que le
−1 −1
rayon spectral de MG NG est proche de 1, en effet on a kxk − x∗ k ≤ k|MG NG k|k ∗
0 ∗ ∗
kx − x k, si Ax = b. Pour pallier cet inconvénient, on équilibre la répartition de
D entre M et N . On utilise un paramètre ω ∈]0; 1[ pour définir Mω = ω1 D + L et
Nω = (1−ω) 0
ω D − U , et l’itération devient (en prenant x = 0) :

xk+1 = Mω−1 ∗ (b + Nω xk )

c’est-à-dire
i−1
X n
X
xk+1
i = ω(bi − Ai,j xk+1
j − Ai,j xkj )/Ai,i + (1 − ω)xki
j=1 j=i+1

Vous expérimenterez cette méthode sur les mêmes exemples (pour pouvoir comparer
toutes ces méthodes), en essayant différentes valeurs de ω.
3 Calcul et visualisation du Laplacien

Cette troisième partie du TP a pour but de visualiser le comportement du Laplacien

d’une fonction u : Ω̄ ⊂ IR2 → IR
½
−∆u(x, y) = f (x, y) pour (x, y) ∈ Ω
u(x, y) = 0 pour (x, y) ∈ ∂Ω

pour plusieurs fonctions f : Ω → IR differentes† . Il s’agit d’appliquer plusieurs com-

mandes de visualisation de MATLABT M , comme mesh, surf, etc.
Le calcul du Laplacien par les différences finies est rappelé en Annexe.
Il y a en effet des commandes MATLABT M déjà prédéfinies pour ce problème. Vous
trouverez dans la suite certaines commandes qui vous seront utiles: Pour plus d’infor-
mation utilisez la commande help.
G = numgrid(’S’, n) Numérote les points de la grille de taille (n,n) d’une région
2-dimensionnelle. La forme de la région est à spécifier.
Ici on choisit le rectangle par ’S’.
D = delsq(G) Construit le Laplacien avec le schéma des différences finies
à 5 points. D sera une matrice sous forme sparse.
spy(D) visualise la densité d’une matrice
Attention, si vous créez une grille G(i, j), i, j = 1, . . . , n pour la discrétisation de taille
n, alors la matrice D du Laplacien est de la taille [(n − 2)2 × (n − 2)2 ]. Il faut donc
calculer la fonction u à l’intérieur de la grille, i.e., uij pour i, j = 2, . . . , n−1, et rajouter
la valeur de u sur le bord.
Après avoir calculé u sur la grille pour une fonction f que vous choisirez, il est demandé
de visualiser le résultat de différentes façons.

†
Comme vous allez certainement le découvrir tout seul on vous indique içi qu’il existe
dans les demos de MATLABT M déjà un exemple pour f ≡ 1.
Annexe: Rappel des différences finies
Problème:
Soit Ω =]x1 , xn [×]y1 , yn [ et une fonction f : Ω → IR de classe C 1 .
On cherche la fonction u : Ω̄ → IR, telle que
½
−∆u(x, y) = f (x, y) pour (x, y) ∈ Ω
(∗)
u(x, y) = 0 pour (x, y) ∈ ∂Ω.

Solution: Différences finies à 5-points

Discrétisation du problème (*) par différences finies:
y

yn n*n points
uij = u(xi,yj)
i,j=1,...,n
uij
yj
h

y1

x
x1 xi xn

On connait u déjà sur le bord de Ω:

u1,j = un,j = ui,1 = ui,n = 0 ∀ i, j

Reste à calculer u à l’intérieure de la grille:

ui−1,j + ui,j−1 − 4ui,j + ui,j+1 + ui+1,j
− = fi,j , i, j = 2, . . . , n − 1
h2
Ce qui donne un système linéaire de taille [(n − 2)2 , (n − 2)2 ] à résoudre:
   u2,2  
f2,2

T I
I T I  ..   .. 
  
. 



. 

 I T I    
    
  u   fi,j 
 ..   i,j
 = −h2  
 .    
    
    
    
 I T I   .
..   .
.. 
I T un−1,n−1 fn−1,n−1
Il s’agit d’une matrice tridiagonal par blocs, où T et I sont des blocs de taille [n−2, n−2]:
   
−4 1 1
 1 −4 1   1 
 ..   .. 
T =
 . 
 I=
 . 

 1 −4 1   1 
1 −4 1

Pour n = 7 on obtient donc une matrice du Laplacien de la forme suivante:

0

10

15

20

25
0 5 10 15 20 25
nz = 105