ENS de Rennes Année 2020 – 2021

Magistère de mathématiques
Algorithme de la factorisation LU.
TP 3 – Méthodes directes de résolution de Ax = b et méthodes Forme compacte
Pour k = 1, . . . , n − 1
itératives. Pour i = k + 1, . . . , n
Ai,k ← Ai,k /Ak,k
Les résolutions de systèmes linéaires n’utiliseront pas la commande [Link] Pour j = k + 1, . . . , n
durant ce TP. Ai,j ← Ai,j − Ai,k Ak,j
Fin pour j
Fin pour i
Théorème (Décomposition LU). Soit A ∈ Mn (K), dont tous les mineurs princi- Fin pour k.
paux sont non nuls, c’est-à-dire que pour tout k ∈ {1, . . . , n} la sous-matrice extraite
(Ai,j )1≤i,j≤k est inversible. • Programmer cet algorithme (on essayera si possible d’obtenir une programmation
Alors il existe un unique couple (L, U ) de matrices carrées d’ordre n tel que : concise, idéalement avec seulement une boucle sur k) et le tester sur la matrice de
· L est triangulaire inférieure avec uniquement des 1 sur la diagonale ; Pascal suivante  
1 1 1 1 1
· U est triangulaire supérieure avec des coefficients diagonaux tous non nuls ; 1 2 3 4 5 
· A = LU .  
P =1 3 6 10 15 .

Une fois cette décomposition connue, la résolution d’un système linéaire Ax = b est 1 4 10 20 35
équivalente à la résolution de deux systèmes linéaires triangulaires, élémentaires à ré- 1 5 15 35 70
soudre :
(Ax = b) ⇐⇒ (Ly = b et U x = y). Le but de la suite de ce TP est d’introduire et de comparer quelques méthodes itéra-
tives pour résoudre le système Ax = b.

Exercice 1. Résolution de systèmes triangulaires Principe : La matrice A est décomposée sous la forme A = M − N avec M facile et
• Programmer les algorithmes de descente et de remontée permettant de résoudre des peu coûteuse à inverser et la solution x est approchée par une suite (xk )k∈N d’éléments
systèmes triangulaires inférieurs (Ly = b) et supérieurs (U x = y) respectivement. de Rn :
x0 ∈ R n donné,
M xk+1 = N xk + b.
Exercice 2. Algorithme de décomposition LU
La matrice M −1 N est fondamentale dans l’analyse de convergence. En effet la solution
L’algorithme de construction des matrices L et U s’appuie sur les opérations d’élimi-
x est l’unique point fixe d’une application affine sur Rn dont la différentielle (constante)
nation du pivot de Gauss. Pour réduire de moitié la complexité en espace (forme com-
a pour matrice M −1 N . En particulier la convergence itérative de la méthode, requise
pacte), l’algorithme agit directement sur la matrice A, la matrice U étant construite
pour tout choix de l’initialisation x0 , est équivalente à la condition suivante sur le rayon
au fil des itérations dans la partie supérieure de A, et la matrice L correspondant aux
spectral :
opérations du pivot étant stockée dans la partie inférieure de A. On a ainsi à l’issue de
ρ(M −1 N ) < 1.
l’algorithme Ui,j = Ai,j si i ≤ j, et Li,j = Ai,j si i > j, les autres coefficients de L et U
étant nuls ou égaux à 1. Quelques possibilités pour construire M et N : écrivons A = D − E − F avec D la partie
diagonale de A, −E sa partie triangulaire inférieure stricte et −F sa partie triangulaire
supérieure stricte. Autrement dit

Dii = Aii , i = 1, . . . , n,
Eij = −Aij , 1 ≤ j < i ≤ n,
Fij = −Aij , 1 ≤ i < j ≤ n,
Année 2020 – 2021 ENS de Rennes

tous les autres termes étant nuls. Pour la construction de ces trois parties utilisera la Pour ce faire, on approche l’équation par différences finies. Soit N un entier naturel,
commande [Link]. on pose h = 1/(N + 1) le pas de discrétisation pour chacune des directions. Pour tout
0 ≤ i, j ≤ N + 1, ui,j désigne une approximation de la valeur u(xi , yj ), les points de la
Méthode de Jacobi : M = D et N = E + F .
grille étant xi = ih, yj = jh. L’EDP est approchée par des différences centrées :
Méthode de Gauss-Seidel : M = D − E et N = F .
ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1
1 1−ω + = f (xi , yj ), 1 ≤ i, j ≤ N,
Méthode de relaxation : M = ωD − E et N = ω D + F avec ω ∈ R+
∗. h2 h2
et les conditions de bord sont données : ui,j = 0 lorsque ij(N + 1 − i)(N + 1 − j) = 0.
a. Construire la matrice A de taille N 2 × N 2 permettant de définir le vecteur des
Exercice 3. Comparaison des méthodes de Jacobi, Gauss-Seidel et de relaxation inconnues U = (ui,j )1≤i,j≤N de RN ×N comme solution du problème
Soit la matrice carrée A de taille n, tridiagonale, intervenant dans le TP précédent,
AU = F.
définie par
−1
 
2 0 ... 0
 .. ..  b. En utilisant les méthodes itératives précédemment mises en œuvre, résoudre ce
−1
 2 −1 . . 
 problème pour le second membre f donné dans l’exemple de code ci-dessous.
A= 0
 .. .. .. 
c. Représenter la solution. On pourra adapter l’exemple suivant pour tracer une sur-
. . .0
face en 3 dimensions :
 
 . ..
 ..

. −1 2 −1
0 ... 0 −1 2 from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
a. Programmer les trois méthodes précedentes.
b. Application : utiliser ces trois méthodes sur la matrice A et le second membre fig = [Link]()
B = (1, 0, . . . , 0, 1)T , pour n = 10. On effectuera 100 itérations. Dans le cas de la ax = [Link](projection = ’3d’)
méthode de relaxation, on choisira ω = 23 . n = 20
X = [Link](0,1,n+2)
c. Dans le cas où n = 20, déterminer le paramètre optimal dans la méthode de
Y = X
relaxation, en représentant le rayon spectral de la matrice d’itération obtenue en fonction
X, Y = [Link](X, Y)
de ω.
def f(x,y):
d. Pour différentes valeurs de n, comparer le nombre d’itérations nécessaires pour return -(x-0.5)**2 - (y-0.5)**2
avoir une précision de 10−12 entre la solution approchée et la solution exacte. Quelle Z = f(X,Y)
relation y a t-il entre le nombre d’itérations et le rayon spectral ? surf = ax.plot_surface(X, Y, Z,cmap = [Link],linewidth = 0,
antialiased = False)
[Link]("La surface f(x,y)=-(x-0.5)**2-(y-0.5)**2")
Exercice 4. Résolution d’un problème de Poisson 2d [Link]()
On souhaite résoudre le problème de Poisson en dimension 2 d’espace où X et Y sont les vecteurs avec les coordonnées des points du plan, et Z est la matrice
2 2 contenant les valeurs f(X,Y).
∂ u ∂ u
+ 2 = f (x, y), (x, y) ∈ [0, 1]2 ,
∂x2 ∂y

assorti des conditions de Dirichlet au bord

u(x, y) = 0, (x, y) ∈ ∂[0, 1]2 .