Cours d’analyse numérique Semestre 1 - 2022/2023

[Link] Première année GSInfotroniques - Mécatronique

Résolution des systèmes linéaires

Exemple de motivation et notations

L’objet de ce chapitre est le développement des méthodes numériques pour la résolution de systèmes linéaires.
Comme dans de nombreux phénomènes physiques on est amené à résoudre des systèmes linéaires de grandes tailles,
il faudrait développer des méthodes qui convergent rapidement et qui se rapprochent le plus possible de la solution.
Problème de réseaux : Soit un réseau fermé de
noeuds Pi et d’arêtes Eij reliant les noeuds Pi et
Pj . C’est le cas de canalisations d’eau, de lignes de
télécommunication, circuits électrique, etc.
A chaque noeud est associé un potentiel Ui et
dans chaque arête circule un fluide dont le débit
qij est proportionnel à la différence de potentiels
qij = kij (Ui − Uj ).
Le réseau est alimenté à partir de noeuds sources P ∗ et P ∗∗ de potentiels associés U ∗ et U ∗∗ donnés. XLa loi de
conservation impose que le débit total est nul si le noeud Pi est isolé du réseau avec ses voisins c.a.d qij = 0
j∈V (i)
X
avec V (i) les voisins de Pi (exemple V (4) = {3, 5, P ∗ }). Donc, on ∀i kij (Ui − Uj ) = 0. Pour simplifier on
j∈V (i)
prend kij = 1. Pour déterminer les potentiels Ui , nous avons à résoudre le système linéaire


 (U1 − U2 ) + (U1 − U5 ) + (U1 − U9 ) = 0
 3 · U2 − U3 − U1 − U5

= 0
.. .

 . = ..
= U ∗∗

3 · U9 − U1 − U6


Ensuite il est facile d’écrire le système sous la forme A · U = b pour ensuite le résoudre c.a.d trouver U .
D’une manière générale, on note Mn (R) l’ensemble des matrices carrés d’ordre n. Soit A ∈ Mn (R) une matrice
inversible et b ∈ Rn . On cherche à résoudre le système linéaire ci-dessous c.a.d trouver x ∈ Rn tel que Ax = b avec
x = t (x1 , . . . , xn ) ; b = t (b1 , . . . , bn ) ; A = (aij )1≤i,j≤n

 a11 x1 + · · · + a1n xn = b1

.. .
 . = ..
an1 x1 + · · · + ann xn = bn


Remarque 1 : Comme A est inversible c.a.d aucune des équations n’est combinaison linéaire des n − 1 autres, il
existe un unique vecteur x solution de Ax = b.
Remarque 2 : Si A matrice de taille (n, m) avec n ̸= m • si n > m si les équations sont indépendantes entre
elles alors il n’y a pas de solution sinon on élimine les équations dépendantes une par une. • si n < m il existe une
infinité de solution. Pour résoudre Ax = b, on distingue deux types de méthodes :
Les méthodes directes : On obtient la solution en un nombre fini d’opérations. L’avantage c’est que nous avons
l’assurance de trouver une solution en un nombre fini d’opérations élémentaires par contre cette solution peut être
numériquement "mauvaise" à cause des erreurs d’arrondi.
Les méthodes itératives : On construit une suite qui tend vers la solution. Si on n’est pas satisfait du résultat,
on peut continuer les itérations pour l’améliorer et dans ce cas plus la précision de la solution est grande plus le
nombre d’itérations est important.

Ecole Nationale d’Ingénieurs de Carthage 1

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique

Résolution des systèmes linéaires

Méthodes Directes

1 Introduction
Dans ce chapitre nous étudions certaines méthodes directes de résolution d’un système de n équations linéaires
à n inconnues Ax = b c.a.d. des méthodes qui mènent à la solution en un nombre fini d’opérations élémentaires
(X,+,-,/).
La méthode de Cramer en est une mais elle est numériquement inacceptable. Elle demande le calcul de (n + 1)
det(Ai )
déterminants puisque xi = ; ∀i = 1, . . . , n. Ai : matrice obtenue en remplaçant la ième colonne de A par
det(A)
le vecteur b.
Chaque déterminant requiert plus que n! opérations élémentaires (matrices de tailles n) donc la méthode requiert
plus que (n + 1)! opérations élémentaires. Pour un exemple de 50 équations et avec un ordianateur exécutant 108
opérations élémentaires par seconde cette méthode nécessitera un temps d’exécution énorme et aberrant ! !
Si on utilisait la méthode de Gauss pour calculer le déterminant, il en faudrait n + 1 application de cette méthode
autant l’appliquer une seule fois et trouver directement la solution.

2 Méthode de Gauss
L’idée étant de transformer la matrice A en   
a11 . . . . . . . . . a1n
  
une matrice triangulaire supérieure ( ou infé- x1 b1
 .. ..  ..   .. 
rieure ) c.a.d avoir un système équivalent à ce 
 0 . . 

  .   . 
système ( ou un système avec une matrice tri-
 .. .. 
 
x
  
 =  bi 
 . . aii ain   i
   
angulaire inférieure). Dans ce cas on trouve les

.. .. .. ..   ..   .. 
.   .   . 

. . .
xi par substitution arrière c.a.d, si on suppose

0 ... . . . 0 ann xn bn
les aii ̸= 0,
n
bn 1 X

on calcule xn puis xn−1 , . . . , x1 avec xn = et xi = bi − aij · xj ∀i = n − 1, . . . , 1.
ann aii
j=i+1
Dans le cas d’une matrice triangulaire inférieure, on retrouve les xi par substitution avant c.a.d on calcule x1 puis
x2 , . . . , xn . Pour transformer le système Ax = b en ce système équivalent, on peut faire des combinaisons linéaires
de lignes et⧸ou permuter des lignes ou des colonnes.

2.1 Description
On pose A(1) = A et b(1) = b et on construit une suite de matrices A(1) , . . . , A(k) , . . . , A(n) et de vecteurs
b(1) , . . . , b(k) , . . . , b(n) .
Première étape :
   
(1) (1) (1) (1) (1) (1)
a11 ... . . . a1j . . . a1n a11 ... . . . a1j . . . a1n
 (1) (1)
 a21 a22 . . . a(1) . . . a(1) (2) (2) (2)
  
2n
  0 a22 . . . a2j . . . a2n 
 2j   
 .. .. ..   .. .. .. 
 . . .
Soit A = A(1)  si a(1) ̸= 0 ; A(2) =  . . .
  
=
 a(1) a(1) . . . a(1) . . . a(1) 11  0 (2) (2) (2)

 i1 i2 ij in

  ai2 . . . aij . . . ain 

 . .. ..   . .. .. 
 .. . .   .. . . 
   
(1) (1) (1) (1) (2) (2) (2)
an1 an2 . . . anj . . . ann 0 an2 . . . anj . . . ann

Ecole Nationale d’Ingénieurs de Carthage 2

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique

avec
(1)
(2) (1) ai1 (1)
aij = aij − (1)
· a1j pour i = 2, . . . , n et j = 1, . . . , n
a11
(1)
(1) (2) (1) ai1 (1)
Les mêmes opérations doivent être effectués sur bi avec bi = bi − b .
(1) 1
a11
A la k ème étape :
   
(1) (1) (1) (1) (1) (1)
a11 . . . . . . a1k ... a1n a ... ... ... a1 k+1 ... a1n
 11 (2) (2) (2)
 0 a(2) (2) (2)
  
22 ... a2k ... a2n   0 a22 . . . . . . a2 k+1 ... a2n 
.. ..
   
.. ..
. .
   
 0 0 .   0 0 . 
. .  .. ..
   
 . .. (k) (k)  (k) (k) (k) 
 . akk ... akn   . . akk ak k+1 ... akn 
 . ..  . ..
   
A(k) =  .. .
(k)
ak+1 k . . . ak+1
(k)
 ; A(k+1) =  .. (k+1)
0 ak+1 k+1
(k+1)
. . . ak+1
 
n . n

 . .. .. ..   . .. .. .. .. 
 .  .
 . . . .  . . . . .
 
 
 . ..   . .. 
 . (k) (k)  . (k+1) (k+1)
 . . aik ... ain  . . 0 ai k+1 ... ain
 
 
 . .. .. ..   . .. .. .. .. 
 .  .
 . . . .  . . . . .
 
 
(k) (k) (k+1) (k+1)
0 0 ank ... ann 0 0 0 an k+1 ... ann
(k)
pour akk ̸= 0 et k = 1, . . . n − 1 avec
(k)
(k+1) (k) aik (k)
aij = aij − (k)
· akj pour i = k + 1, . . . , n et j = k, . . . , n
akk
(k)
(k) (k+1) (k) aik (k)
Les mêmes opérations doivent être effectués sur bi avec bi = bi − b .
(k) k
Tant que le pivot ne s’annule pas
akk
on continu, à la fin on obtient un système triangulaires supérieur.

Nombre d’opérations élémentaires :

Pour passer de A(1) , . . . , A(n) , il faudrait (n − 1)2 + (n − 2)2 + . . . + 1 additions, le même nombre en multiplications
n−1
X 3n(n − 1)(2n − 1)
et le même en divisions. Donc 3 · i2 = opérations élémentaires.
6
i=1

3n(n − 1)
Pour passer de b(1) , . . . , b(n) , il faudrait 3 · (n − 1) + (n − 2) + . . . + 1 = opérations élémentaires.
2
Donc la méthode de Gauss nécessite un nombre d’opérations élémentaires de l’ordre de n3 incomparable avec la
méthode de Cramer. Exemple pour n = 100 la solution est donnée en moins d’une seconde.

2.2 Méthode du pivot partiel

(k)
Dans la pratique, il se peut qu’à l’étape k, le pivot akk = 0. Dans ce cas, on cherche dans la même colonne s’il
(k)
y’a une ligne i0 > k où il n’est pas nul ai0 k ̸= 0 et on fait une permutation entre les 2 lignes (En permutant
les lignes de la matrice A, on permute aussi celles du second membre b). Par contre on ne choisit pas n’importe
quelle ligne. En effet, dans un ordinateur plus on divise par des nombres petits plus on risque de faire des erreurs
d’arrondi importantes. On choisit donc le nouveau pivot de manière à minimiser ces erreurs. On choisit i0 tel que
| ai0 k |= maxi=k,...,n | aik | c.a.d celui dont la valeur absolue est la plus grande.

Ecole Nationale d’Ingénieurs de Carthage 3

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique
     
ϵ 1 x 1
Exemple : Soit 0 < ϵ < 1 et le système linéaire = . On applique Gauss directement
1 1 y 2
     
ϵ 1 x 1
L2 ← L2 − 1ϵ L1 donc = .
0 1 − 1ϵ y 2 − 1ϵ
2ϵ − 1 1−y
Par substitution arrière, on calcule y = et ensuite x = .
ϵ−1 ϵ
−0, 999 . . . 9
A.N : Si la précision de l’ordinateur est 10−15 et ϵ = 10−20 , y = = 1 et donc x = 0. Donc la solution
−0, 999 . . . 9
donnée en appliquant Gauss directement est (0, 1) qui représente une solution non acceptable vu qu’elle ne vérifie
pas le système ! Ceci est dû au fait d’avoir pris un pivot très petit plus petit que la précision de l’ordinateur.

Appliquons maintenant la stratégie du pivot partiel. Dans ce cas on fera tout
 d’abord
  une permutation des deux
1 1 x 2
lignes et on prendra comme pivot 1. Le système est alors = . L2 ← L2 − ϵL1 donne
ϵ 1 y 1
     
1 1 x 2 1 − 2ϵ
= . Par substitution arrière, on calcule y = ce qui donne en A.N. y = 1 et
0 1−ϵ y 1 − 2ϵ 1−ϵ
ensuite x = 2 − y = 1. Donc la stratégie du pivot partiel donne (1, 1) comme solution du système ce qui représente
une solution acceptable.

2.3 Méthode du pivot total

On peut aussi envisager de chercher le nouveau pivot dans tout le carré situé sous le pivot. On intervertit dans ce
cas les lignes et les colonnes si nécessaires (En intervertissant les colonnes, on renumérote les inconnus xi de la
solution). Donc on choisit i0 et j0 tels que | ai0 j0 |= maxi,j=k,...,n | aij |.

Remarque
Tout échange ou modification de lignes doit se faire sur A(k) et b(k) . L’intérêt de la méthode suivante est qu’elle ne
modifie pas le second membre b. Elle est donc très avantageuse lorsqu’on doit résoudre plusieurs systèmes ayant
tous la même matrice A et des seconds membres différents. C’est la factorisation L · U

3 Factorisation L · U
Supposons que nous sommes capables d’écrire la matrice A sous la forme d’un produit de deux matrices : L(Lower)
triangulaire
 inférieure et U (U pper) triangulaire supérieure, le système Ax = b peut s’écrire L · U x = L · (U x) = b
(1) Ly = b
⇔ .
(2) U x = y
Donc si on connaît L et U , (1) se résout facilement par substitution avant et (2) par substitution arrière.

Théorème 1 (Décomposition LU d’une matrice)

Soit A ∈ Mn (R). Si tous les pivots de Gauss ne s’annulent pas, il existe un unique couple de matrices (L, U ) : L
triangulaire inférieure avec des 1 sur la diagonale et U triangulaire supérieure tels que A = L · U .
Démo : Cette décomposition peut se calculer facilement à partir de la méthode de Gauss. En effet, on peut montrer
que le passage de A(k) à la matrice A(k+1) , ∀k = 1, . . . , n − 1 revient à multiplier à gauche la matrice A(k) par la
matrice E (k) de taille (n, n) suivante :

Ecole Nationale d’Ingénieurs de Carthage 4

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique

1 0 ... ... ... ... 0

 
..
0 . 0 ... ... ... ... 0
 
 
 .. .. 
. 1 .
 
 
 .. .. 
(k)
(k)

. −ℓk+1;k .
 aik
Soit E = ; avec ℓi;k = pour i = k + 1, . . . , n
 
.. .. .. ..  (k)

 . . 0 . . 
 akk
 .. .. .. .. .. 

 . −ℓi;k . . . . 

 .. .. .. .. .. 
 . . . . . 0 
0 −ℓn;k 0 ... ... 0 1
Donc, A(k+1) = E (k) · A(k) et par récurrence A(n) = E (n−1) × E (n−2) × · · · × E (1) · A(1) .
En posant M = E (n−1) × · · · × E (1) et puisque A(1) = A, nous avons A(n) = M · A avec A(n) la dernière matrice
triangulaire supérieure obtenue en appliquant l’élimination de Gauss et M triangulaire inférieure puisque c’est le
produit de matrices triangulaires inférieures.
n−1
Y
La matrice M étant inversible puisque det(M ) = detE (k) = 1. L’inverse d’une triangulaire inférieure étant une
k=1
triangulaire inférieure, en posant L = M −1 et U = A
(n) , on aura L · U

= A. Puisqu’on

connait les termes de U ,
−1 (1) −1 (2) −1 (n−1) −1
cherchant la structure de la matrice L ; M = E × ··· × E ··· E avec

1 0 ... ... ... ... 0

 
..
0 . 0 ... ... ... ... 0
 
 
 .. .. 
. 1 .
 
 
 .. .. 
. ℓk+1;k .

−1  
−1 (k)
Soit E (k) =  ; ce qui vérifie E (k) ·E = In pour k = 1, . . . , n−1
 
.. .. .. ..

 . . 0 . . 

 .. .. .. .. .. 

 . ℓi;k . . . . 

 .. .. .. .. .. 
 . . . . . 0 
0 ℓn;k 0 ... ... 0 1

−1
En calculant le produit des E (k) , on aura
 
1 0 ... ... ... 0
 .. 
 ℓ2;1
 1 . 0  
 .. .. .. ..  (k)
 . ℓ3;2 . . .  aik
L=  .. .. .. .. .. 
 ; avec ℓi;k = (k)
pour k = 1, . . . n − 1 et i = k + 1, . . . , n
 . . . . .  akk
 .. ..
 
.. ..
. .

 . . 0 
ℓn;1 ℓn;2 ℓn;n−1 1

Nous avons donc démontrer l’existence de L triangulaire inférieure avec des 1 sur la diagonale et U triangulaire
supérieure tels que A = L · U . Pour l’unicité (voir TD avec corrigé)
Maintenant regardons quand est ce que tous les pivots ne s’annulent pas et donc quand est ce qu’on peut utiliser
la factorisation LU pour une matrice.

Ecole Nationale d’Ingénieurs de Carthage 5

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique

Théorème 2
Tous les pivots de la méthode de Gauss ne s’annulent pas si et seulement si toutes les sous-matrices
 principales
 Ãk
a11 . . . a1k
de A ∀k = 1, . . . n sont inversibles. La k ème sous-matrice de A est de la forme Ãk =  ... .. .. .

. . 
ak1 . . . akk

En effet, à l’étape k de la méthode de Gauss, on a déjà construit L̃(k) et Ũ (k) tel que Ãk = L̃k · Ũk avec
   (1) (1)

1 0 ... 0 a11 . . . . . . a1k
. .. .. 
 ℓ2;1 . .  0 a(2) (2)  k
 
. .  22 . . . a2k  Y (i)
L̃k =  .

. .
 et Ũk =  .
 . .. .. ..  ; det(Ãk ) = det(L̃k ) · det(Ũk ) = 1 ·
 aii .
 .. .. ..  . . . . 

0  i=1
(k)
ℓk;1 . . . ℓk;k−1 1 0 . . . 0 akk
(k)
Donc Ãk inversible ⇔ det(Ãk ) ̸= 0 ⇔ tous les pivots akk des étapes précédentes sont non nuls.

4 Méthode de Cholesky
4.1 Matrice symétrique définie positive
• Une matrice A ∈ Mn (R) est symétrique si t A = A.
• Une matrice A ∈ Mn (R) est définie positive ssi < A · x, x >= t x · A · x > 0, ∀x ∈ R∗ n.
• Une matrice A ∈ Mn (R) est définie positive ssi les valeurs propres de A sont > 0.
• Les mineurs principaux dominants de A qui représentes les déterminants des sous-matrices principales Ãk sont
strictement positifs.

4.2 Décomposition de Cholesky

Dans le cas où la matrice A est symétrique définie positive, la condition du théorème précédent de la décomposition
LU est satisfaite (et donc les déterminants
 des sous-matrices
 principales Ãk sont ̸= 0) et la matrice A admet une
1 0 ... 0
 
U11 . . . . . . U1n
.. 
 ℓ2;1 . . . ..  0 U22 . . . U2n 

. . 
factorisation LU = A avec L =  .  et U =  .. .


. .  .. . . . .
 .. .. ..  . . . . 

0 
ℓn;1 . . . ℓn;n−1 1 0 . . . 0 U nn
k
Y
Comme A est définie positive, on a det(Ãk ) = det(L̃k )·det(Ũk ) = 1· Uii , ∀k = 1, . . . n et les Uii > 0, ∀i = 1, . . . n.
i=1
 √
√1 0 ... 0
  
U11 0 ... 0 U11
 .. .. ..   .. .. .. 
 0 . . .  qui est inversible avec D−1 = 
  0 . . . 
On pose D = 
 .. .. ..   .. .. .. ..
.

. . . .
. √0 . .
   
0 ... 0 Unn 0 ... 0 √1
Unn
√
A = L · U = L · D · D−1 · U = R · t B avec R = L · D qui est
√ triangulaire inférieure avec Uii en diagonale donc
t −1
inversible et B = D · U est triangulaire supérieure avec Uii en diagonale donc inversible.
Comme A est symétrique donc A = t A ⇔ R · t B = B · t R ⇔ B −1 · R · t B · ( t B)−1 = B −1 · B · t R · ( t B)−1 donc
B −1 · R = t R · ( t B)−1 .
De plus, B −1 · R est triangulaire inférieure puisque c’est le produit de deux matrices triangulaires inférieures

Ecole Nationale d’Ingénieurs de Carthage 6

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]
Cours d’analyse numérique Semestre 1 - 2022/2023
[Link] Première année GSInfotroniques - Mécatronique

et t R · ( t B)−1 triangulaire supérieur puisque c’est le produit de deux

√ matrices triangulaires supérieures donc
B −1 · R = ∆ : matrice diagonale d’éléments diagonaux δii = √U 1
· Uii = 1. Donc B −1 · R = In et B = R et
ii
donc il existe une matrice R triangulaire inférieure avec des éléments de la diagonale strictement positifs tel que
A = R · t R.
De plus cette matrice est unique. En effet supposons qu’il existe R1 triangulaire inférieure avec des éléments de la
′ > 0 telles que A = R · t R =
diagonale rii > 0 et R2 triangulaire inférieure avec des éléments de la diagonale rii 1 1
t −1 t t −1 −1 t t −1
R2 · R2 ⇔ R1 · R1 · R1 · ( R2 ) = R1 · R2 · R2 · ( R2 ) .
Donc t R1 · ( t R2 )−1 triang. sup. = R1−1 · R2 triang. inf. = ∆ : matrice diagonale d’éléments diagonaux δii = rii
′
rii
=
′
rii ′ )2 et r = r ′ . Donc δ = 1 et R−1 · R = I et R = R .
⇔ (rii )2 = (rii
rii ii ii ii 1 2 n 1 2
On a donc démontré le théorème suivant :

Théorème 3 (Décomposition de Cholesky)

Soit A ∈ Mn (R) une matrice symétrique définie positive, il existe une unique matrice R triangulaire inférieure
avec des éléments de la diagonale strictement positifs telle que A = R · t R.

4.3 Calcul de la matrice R

D’après la construction de la matrice R et si une décomposition LU a été effectuée pour la matrice A, la matrice
R = L · D avec L la matrice triangulaire inférieure avec des 1 √ sur la diagonale regroupant les coefficients de la
méthode d’élimination de Gauss et D la matrice diagonale avec Uii c’est à dire les racines des pivots de Gauss.
Sinon pour calculer R à partir de A, on peut poser des éléments inconnus pour R et les déterminer par identification
des éléments de R · t R à ceux de A.

t (1) R·y = b
Une fois la factorisation de Cholesky réalisée, le système Ax = b ⇔ R · Rx = b ⇔ tR · x = y .
(2)
Donc si on connaît R, (1) se résout facilement par substitution avant et (2) par substitution arrière.

Ecole Nationale d’Ingénieurs de Carthage 7

Charguia II - 2035 Tunis-Carthage • Contact : [Link]@[Link]