Analyse numérique

Thomas Harbreteau

13 avril 2020

Notes du cours de Benjamin Boutin.

Université de Rennes 1, année 2019/2020.

Document sous license Art Libre ([Link]

ANU SOMMAIRE

Sommaire

1 Résolution de systèmes linéaires 2

1 Rappels et notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Réduction en dimension finie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
3 Propriétés spectrales des matrices hermitiennes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
4 Normes subordonnées et rayon spectral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
5 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
6 Exemple important : matrice du laplacien 1d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Méthodes directes 14
1 Remarques préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Décomposition de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3 Méthodes itératives 18
1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2 Méthodes usuelles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Critères explicites de convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4 Résolution de systèmes linéaires au sens des moindres carrés . . . . . . . . . . . . . . . . . . . . . . . . 21
5 Méthodes variationelles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4 Approximation spectrale 27
1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2 Méthodes numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

1
Chapitre 1

Résolution de systèmes linéaires

1 Rappels et notations
1.1 Rappels
Soit A ∈ Mn (K) et λ ∈ σ(A).
• Polynôme caractéristique : χA = det(A − XIn ).
• Multiplicité algébrique de λ : multiplicité en tant que racine du polynôme caractéristique.
• Multiplicité géométrique de λ : dimension du sous-espace propre associé à λ.
• Valeur propre défective : multiplicité géométriques et algébriques différentes.

1.2 Notations
• Tnu (K) l’ensemble des matrices triangulaires supérieures.
• Un (K) l’ensemble des matrices unitaires, i. e. dont les colonnes forment une base orthonormée pour le produit
scalaire hermitien (analogue aux matrices orthogonales réelles). Caractérisation matricielle : U ∗ U = U U ∗ = In .

2 Réduction en dimension finie

Théorème 1.2.1 (Réduction de Jordan) Soit A ∈ Mn (K), supposée trigonalisable sur K,
 
Jk1 (λ1 ) 0 ... 0
. ..
Jk2 (λ2 ) . .
 
 0 .
P −1 AP = 

∃P ∈ GLn (K),  .
,
 .. .. .. 
. . 0 
0 ... 0 Jkm (λm )
où
• les (λi )1≤i≤m non-nécéssairement distincts deux à deux,
• pour tout i ∈ {1, . . . , m}, Jki (λi ) désigne un bloc de Jordan de taille ki ≥ 1 :
 
λ 1 0 ... 0
. .. .. .
. .. 

0 λ
 
Jk (λ) =  ... . . . . . . .. .
 
 . 0 
. .. ..
 ..

. . 1
0 ... ... 0 λ

2.1 Réduction en base orthonormée

Proposition 1.2.2 (Produit scalaire hermitien) ∀x, y, z ∈ C∗ , ∀λ, µ ∈ C, ∀A ∈ Mn (C),
• hx, y + zi = hx, yi + hx, zi,
• hx, yi = hy, xi,
• hλx, µyi = λ̄µhx, yi,
• hAx, yi = hx, A∗ xi, où A∗ = ĀT (définition de l’adjoint),
• |hx, yi| ≤ kxk2 kyk2 .

2
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Proposition 1.2.3 (Projection orthogonale) On considère F un sous-espace de Kn et (e1 , . . . , en ) une base

orthonormée de F . Alors la projection orthogonale sur F parallèlement à F ⊥ , notée pF , est obtenue comme :
m m
!
X X
∀x ∈ Kn , pF (x) = hei , xi = ei e∗i x.
i=1 i=1

Preuve. "m #
m
X m
X m
X X
∗ ∗ ∗
hei , xiei = (ei x)ei = ei (ei x) = (ei ei ) x.
i=1 i=1 i=1 i=1


n
Proposition 1.2.4 (Orthonormalisation de Gram-Schmidt) Soit (vi )1≤i≤n une base de K . On construit une
base orthonormée (ei )1≤i≤n de Kn par projections successives sur les directions précédemment construites :

u1 := v1 , e1 := u1 /kuk,

et pour tout k ∈ {2, . . . , n},

k−1
X hui , vk i uk
uk := vk − pVect(e1 ,...,ek−1 ) (vk ) = vk − ui , ek := .
i=1
hui , ui i kuk k

Proposition 1.2.5 (Factorisation QR) Soit A ∈ GLn (C), alors il existe Q ∈ Un (C) et R ∈ Tnu (C) telles que
A = QR. Cette décomposition est de plus unique si on suppose les termes diagonaux de R réels positifs.
Preuve.
k−1
X hui , vk i k−1
X
∀k ∈ {1, . . . , n}, vk = uk + ui = kun kek + hei , vk iei .
i=1
hui , ui i i=1

Ainsi,
.. .. . .. ku k
    
. .   .. . 1 (hei , vj i)
..
 v1
A= ... vn   =  e1 ... en   .
 
  .
.. .. .. .. (0) kun k
. . . .

∗
Théorème 1.2.6 (Schur) Soit A ∈ Mn (C), alors il existe U ∈ Un (C) et T ∈ Tnu (C) telles que U AU = T .
Preuve. P admet une factorisation QR, notée P = QR, où Q ∈ Un (C) et R ∈ Tnu (C). Alors

A = QRT (QR)−1 = QRT R−1 Q−1 = Q(RT R−1})Q∗ .

| {z
∈Tnu (C)

2.2 Réduction des matrices normales

Définition 1.2.7 (Matrice normale) Une matrice A ∈ Mn (C) est dite normale si AA∗ = A∗ A.
Exemple 1.2.8 Sn (R) ; Hn (C) hermitiennes complexes (A∗ = A) ; On (R) ; Un (C).
Lemme 1.2.9 Toute matrice triangulaire (supérieure ou inférieure) et normale est diagonale.
Preuve. Soit A ∈ Tnu (C).
 i i
X X
 (A∗ A) |a2i,j |

i,i = Ā A
j,i j,i =
∀i ∈ {1, . . . , n},
 j=1 j=1
(AA∗ )i,i = |ai,i |2 .


3
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Théorème 1.2.10 Une matrice A ∈ Mn (C) est normale si et seulement si elle est diagonalisable en base ortho-
normée.
Preuve. ( ⇐= ) : Supposons A diagonalisable en base orthonormée. On note A = U DU ∗ , où U ∈ Un (C) et D
diagonale. Alors  ∗
A A = (U DU ∗ )∗ (U DU ∗ ) = U D∗ (U ∗ U )DU ∗ = U D∗ DU ∗
AA∗ = U DD∗ U ∗ .
Ainsi,
AA∗ = A∗ A ⇐⇒ DD∗ = D∗ D ⇐⇒ ∀i ∈ {1, . . . , n}, Di,i D̄i,i = D̄i,i Di,i .
( =⇒ ) : Soit A ∈ Mn (C) normale, AA∗ = A∗ A. Par théorème de Schur, il existe U ∈ Un (C) et T ∈ Tnu (C)
telles que A = U T U ∗ . Du fait de l’égalité AA∗ = A∗ A, on trouve T T ∗ = T ∗ T après calculs. Ainsi, T est normale et
triangulaire supérieure, donc est diagonale par le lemme précédent. 
Proposition 1.2.11 Les valeurs propres d’une matrice normale sont
• réelles si A ∈ Sn (R), ou plus généralement si A ∈ Hn (C),
• réelles positives si A ∈ Sn+ (R) ou si A ∈ Hn+ (R),
• réelles strictement positives si A ∈ Sn++ (R) ou si A ∈ Hn++ (R),
• imaginaires pures si A est anti-symétrique réelle ou anti-hermitienne,
• complexes de module 1 si A ∈ On (R) ou plus généralement si A ∈ Un (C).
Preuve.
• Valeurs propres de A ∈ Un (C) ? Soit A normale, alors A = U DU ∗ avec U unitaire et D vérifie D∗ D = DD∗ = In
(car A∗ A = AA∗ = In ). On déduit de cette dernière égalité que

∀i ∈ {1, . . . , n}, |Di,i |2 = 1,

donc σ(A) = σ(D) = {z ∈ C | |z| = 1}.



3 Propriétés spectrales des matrices hermitiennes

Définition 1.3.1 (Quotient de Rayleigh) Soit A ∈ Hn (C), on définit l’application

(Cn )∗ −→ R
RA : hAx, xi .
x 7−→
hx, xi

Preuve. Si A ∈ Hn (C),

∀x ∈ Cn , hAx, xi = (Ax)∗ x = x∗ (A∗ x) = hx, A∗ xi = hx, Axi,

car A = A∗ . Par caractère hermitien du produit scalaire, hAx, xi = hx, Axi, donc ici, hx, Axi = hx, Axi ∈ R. Ainsi,
RA (x) ∈ R. 
Théorème 1.3.2 Soit A ∈ Hn (C) de valeurs propres (non nécéssairement positives) λmin = λ1 ≤ λ2 ≤ · · · ≤ λn =
λmax , alors
max RA (x) = λmax , et min RA (x) = λmin .
x6=0 x6=0

Preuve. Soient x ∈ Cn \ {0} et z ∈ C \ {0},

hAzx, zxi z̄zhAx, xi

RA (zx) = = = RA (x).
hzx, zxi z̄zhx, xi

Étudier les extremums de RA sur Cn \ {0} se ramène à les étudier sur S := {x ∈ Cn | kxk2 = 1}.
S étant compacte et RA continue sur S, RA est bornée et atteint ses bornes.
A étant hermitienne, elle est diagonalisable en base orthonormée. Notons A = U DU ∗ , avec U ∈ Un (C) et D
diagonale réelle (car A est hermitienne), D = diag(λ1 , . . . , λn ) avec λ1 ≤ · · · ≤ λn . Alors

hU DU ∗ x, xi hDU ∗ x, U ∗ xi
RA (x) = = = RD (U ∗ x).
U U ∗ x, xi hU ∗ x, U ∗ xi
L’application
S −→ S
σ:
x 7−→ U ∗x

4
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

définit une bijection sur S. Les extremums de RA sur S sont ceux de RD . Soit y ∈ S,
n n
hDy, yi X X
RD (y) = = hDy, yi = y ∗ D∗ y = y ∗ Dy = ȳi λi yi = λi |yi |2 .
hy, yi i=1 i=1

Comme y est de norme 1,

n
X n
X
λmin |yi |2 ≤ RD (y) ≤ λmax |yi |2 .
i=1 i=1
| {z } | {z }
=1 =1

Les égalités sont atteintes à gauche pour y = (1, 0, . . . , 0) et à droite pour y = (0, . . . , 0, 1)T .
T

Remarque 1.3.3 On peut bien permuter les vecteurs colonnes de D. Soit P une matrice de permutation, D̃ =
P DP −1 = P DP ∗ , alors
A = U DU ∗ = U P ∗ D̃P U ∗ = ( U P ∗} )D̃(U P ∗ )∗ .
| {z
unitaire

4 Normes subordonnées et rayon spectral

4.1 Définitions
Définition 1.4.1 (Normes usuelles) On considère x ∈ Cn et A ∈ Mn (C).
n
X n
X
kxk1 = |xi |, kAk1 = max |Ai,j |.
1≤j≤n
j=1 i=1

n
X
kxk∞ = max |xj |, kAk∞ = max |Ai,j |.
1≤j≤n 1≤i≤n
j=1

4.2 Rayon spectral

Définition 1.4.2 (Rayon spectral) Soit A ∈ Mn (K). On appelle rayon spectral de A la quantité

ρ(A) = max{|λ| | λ ∈ σ(A)}.

Lemmep1.4.3 Si A est normale, alors kAk2 = ρ(A). Plus généralement, pour toute matrice A ∈ Mn (C), on a
kAk2 = ρ(A∗ A).
Théorème 1.4.4 Soit A ∈ Mn (K).
1. Pour toute norme subordonnée k · k, on a ρ(A) ≤ kAk.
2. Réciproquement, pour tout ε > 0, il existe une norme subordonnée k · kε telle que kAkε ≤ ρ(A) + ε.
Ainsi,
ρ(A) = inf kAk.
k·k subordonnée

Preuve.
1. Soit λ ∈ σ(A) telle que ρ(A) = |λ|, et x ∈ Cn \ {0}, Ax = λx.

kAxk = kλxk = |λ|kxk = ρ(A)kxk,

et
kAxk kAyk
= ρ(A) = max = kAkMn (C) .
kxk y∈C \{0} kyk
n

2. Soit A ∈ Mn (C). Par théorème de Schur, il existe U ∈ Un (C) et T ∈ Tnu (C) telles que A = U T U ∗ . Alors
ρ(A) = ρ(T ), mais a-t-on une relation entre kAk et kT k, pour une norme quelconque ? On calcule
 
 j−1
X 
kT1 k = max |λj | + |Ti,j | .
 
1≤j≤n  
i=1
| {z }
à « éliminer »

5
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Soit η > 0, on note  

η (0)
 η2 
D(η) :=  ,
 
..
 . 
n
(0) η
et on définit
(η −i+j Ti,j )
 
λ1
T (η) := 
 .. .

.
(0) λn
On remarque que T (η) →(η→0) diag(λ1 , . . . , λn ), donc

j−1
!
X
kT (η)k1 = max |λj | + η η j−1−i |Ti,j | −−−→ max |λj | = ρ(A).
j η→0 j
i=1

Soient ε > 0 et η suffisamment petit de sorte que kT (η)k1 ≤ ρ(A) + ε. De plus, on a

A = (U D(η))T (η)(D(η)−1 U ∗ ).

On considère sur Cn la norme

∀x ∈ Cn , kxkε = kD(η)−1 U ∗ xk1 .
Ainsi,

kAxkε = kD(η)−1 U ∗ Axk1

= k (D(η)−1 U ∗ )(U D(η)) T (η)(D(η)−1 U ∗ )xk1
| {z }
=In

= kT (η)D(η)−1 U ∗ xk1
≤ kT (η)k1 kD(η)U ∗ xk1
≤ (ρ(A) + ε)kxk1 .

Par conséquent,
ρ(A) = inf kAk.
k·k subordonnée

4.3 Utilisation du rayon spectral

Théorème 1.4.5 Soit A ∈ Mn (K), il y a équivalence entre les propositions suivantes.
1. lim(k→+∞) Ak = 0.
2. Il existe une norme subordonnée k · k telle que kAk < 1.
3. ρ(A) < 1.
Preuve. 2) =⇒ 1) : Ak → 0 ⇐⇒ kAk k → 0, mais kAk k ≤ kAkk → 0.
2) ⇐⇒ 3) : Conséquence du théorème précédent.
1) =⇒ 3) : Par contraposée, supposons que ρ(A) ≥ 1. Alors il existe λ ∈ C et x ∈ Cn tels que |λ| ≥ 1, x 6= 0 et
Ax = λx. Par récurrence, pour tout k ∈ N, Ak x = λk x, donc

kAk xk = |λk |kxk ≥ kxk =

6 0.

Ainsi,
kAk yk
kAk k = sup ≥ 1,
y6=0 kyk
donc (Ak )k ne peut tendre vers 0. 
Corollaire 1.4.6 Soit A ∈ Mn (K), pour toute norme subordonnée k · k, on a

lim kAk k1/k = ρ(A).

k→+∞

6
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Preuve. Soient λ ∈ C et x ∈ Cn tels que x 6= 0, Ax = λx et |λ| = ρ(A). Alors

kAk xk = |λ|k kxk = ρ(A)k kxk,

donc kAk k ≥ ρ(A)k , d’où kAk k1/k ≥ ρ(A), et ce pour tout k ∈ N.

Soit ε > 0, On souhaiterait obtenir pour k assez grand kAk k1/k ≤ ρ(A) + ε. Posons
A
Aε = .
ρ(A) + ε

On a alors ρ(Aε ) = ρ(A)/(ρ(A) + ε) < 1, donc Akε →(k→+∞) 0. À partir d’un certain rang k0 , kAkε k ≤ 1, or
 k
1
kAkε k = kAkk .
ρ(A) + ε

Par conséquent, pour k ≥ k0 , kAk k ≤ (ρ(A) + ε)k , d’où kAk k1/k ≤ ρ(A) + ε. 
Théorème 1.4.7 Soit A ∈ Mn (K), il y a équivalence entre les propositions suivantes.
1. La suite (Ak )k≥0 est bornée.
2. ρ(A) ≤ 1 et les valeurs propres de A de module 1 sont toutes semi-simples.
Preuve. Idée essentielle : T la réduite de Jordan de A, de blocs de Jordan les Jkj (λj ), alors

(Ak )k bornée ⇐⇒ ∀j ∈ {1, . . . , n}, (Jkj (λj )k )k est bornée.

or   
k k−i

k
λj kλk−1
j ... λ
i j 
..
 
k
 .. .. 
Jkj (λj ) =  . . . .
 
 .. k−1

 . kλj 
k
λj
Ainsi,
• si |λj | < 1, Jkj (λj )k → 0,
• si |λj | > 1, |λkj | → +∞, donc (Jkj (λj )k ) n’est pas bornée.

Exemple 1.4.8
• Soit  
2 1 3
A := 0 1 5 ,
0 0 1/2
alors ρ(A) = 1 > 1, donc (Ak )k n’est pas bornée.
• Soit  
i 1 0
B := 0 1 1 ,
0 0 1
alors ρ(B) = 1 et σ(B) = {1, i}. La valeur propre 1 est défective (associée à un bloc de Jordan de taille 2 > 1),
donc (B k )k n’est pas bornée.
• Soit  iπ 
e 0 0
C :=  0 eiπ/6 0 .
0 0 eiπ/6
Alors ρ(C) = 1 et (C k )k est bornée.

5 Conditionnement
On se place dans Kn muni d’une norme vectorielle k · k et Mn (K) muni de la norme induite k · k.
Définition 1.5.1 (Conditionnement) Étant donnée A ∈ GLn (K), on appelle conditionnement de A dans la norme
k · k la quantité
cond(A) = kAkkA−1 k.

7
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Proposition 1.5.2
1. cond(A) ≥ 1.
2. ∀α ∈ K∗ , cond(αA) = cond(A).
3. ∀A, B ∈ GLn (K), cond(AB) ≤ cond(A) cond(B).
Proposition 1.5.3 (Conditionnement en norme 2)
1. Soit A ∈ GLn (K), s
max{|λ| | λ ∈ σ(A∗ A)}
cond2 (A) = .
min{|λ| | λ ∈ σ(A∗ A)}
2. Si de plus A est normale, alors
max{|λ| | λ ∈ σ(A)}
cond2 (A) = .
min{|λ| | λ ∈ σ(A)}
3. En particulier, si A ∈ Un (C), ou si A ∈ On (R), alors cond2 (A) = 1.
Preuve. Si A est normale, on sait que ρ(A) = kAk2 = max{|λ| | λ ∈ σ(A)}. De plus, si A ∈ GLn (C),
ρ(A−1 ) = max{|λ−1 | | λ ∈ σ(A)} = min{|λ| | λ ∈ σ(A)}−1 .
Mais A−1 est normale ((AA∗ )−1 = (A∗ A)−1 ), donc ρ(A−1 ) = kA−1 k2 . Ainsi,
max |λ|
cond2 (A) = .
min |λ|

n n
Théorème 1.5.4 Soient A ∈ GLn (K), b ∈ K et x ∈ K tels que Ax = b.
1. Soient δx ∈ Kn et δb ∈ Kn tels que A(x + δx) = b + δb, alors
kδxk kδbk
≤ cond(A) .
kxk kbk
2. Soient δx ∈ Kn et δA ∈ Mn (K) tels que A + δA ∈ GLn (K) et (A + δA)(x + δx) = b, alors
kδxk kδAk
≤ cond(A) .
kx + δxk kAk
Preuve. Soient A ∈ GLn (C), b ∈ Cn et x ∈ Cn tels que Ax = b. Soit δb ∈ Cn une perturbation. Notons y := x + δx
l’unique solution de Ay = b + δb. Donc si Ax + Aδx = b + δb, alors δx = A−1 δb. Par conséquent,
kδxk = kA−1 δbk ≤ kA−1 kkδbk,
donc
kbk = kAxk ≤ kAkkxk,
d’où
kδxk kδbk
≤ kAkkA−1 k .
kxk kbk

Exemple 1.5.5 Probèmes mal conditionnés :
• Interpolation polynomiale : Soient (xi )1≤i≤n et (yi )1≤i≤n des complexes. On suppose les (xi )i deux à deux
distincts. On cherche le polynôme P ∈ C[X] de degré minimal vérifaint pour tout i ∈ {1, . . . , n}, P (xi ) = yi .
Il existe un unique P ∈ Cn−1 [X] solution. On peut le chercher dans la base canonique de Cn−1 [X],
n−1
X
P = aj+1 X j .
j=0

Le vecteur A = (a1 , . . . , an )T ∈ Cn est alors solution de

 
y1
 .. 
V A =  .  , avec V = (xij−1 )1≤i,j≤n la matrice de Vandermonde (inversible).
y1
cond(V ) est très grand pour n assez grand. Conséquence (en présence d’une arithmétique approchée) : la solution
obtenue numériquement est loin d’être bonne. Sur le test, cond(V ) ' 1024 , avec une erreur machine de 10−16 .
On peut donc s’attendre à kδxk/kxk ' 108 .
Pour contourner cette difficulté, on reformule le problème au niveau algébrique. Cela revient à changer de base
dans la formulation d’origine. Par exemple :

8
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

— Base de Lagrange associée aux points (xi )1≤i≤n ,

n n
Y X − xj X
Li = , P = yi Li .
x − xj
j=1 i i=1
j6=i

Dans ce cas, la formulation « algèbre linéaire »revient à inverser la matrice identité, de conditionnement 1.
— Base de Newton :
P = b0 + b1 (X − x1 ) + b2 (X − x1 )(X − x2 ) + · · · + bn−1 (X − x1 ) . . . (X − xn−1 ).
L’algorithme des différences divisées (mémo Scilab) permet de trouver les (bi )1≤i≤n−1 (de type « taux
d’accroissement »).

6 Exemple important : matrice du laplacien 1d

De nombreux problèmes sur des applications des mathématiques font intervenir l’opérateur laplacien
∂2u ∂2u ∂2u
∆u = + 2 + 2 en dimension 3,
∂x2 ∂y ∂z
comme
∂u
• l’équation de la chaleur : = K∆u,
∂t
∂2u
• l’équation des ondes : − c2 ∆u = 0.
∂t2
Exemple 1.6.1 f ∈ C 0 ([0, 1]) est donnée, on considère le système
−u00 (x) = f (x) ,

x ∈]0, 1[
(S) : .
u(0) = u(1) = 0
Si f = 0, l’unique u ∈ C 2 ([0, 1]) est u = 0. Sinon, on peut trouver une formule par intégrations successives (en
dimension plus grande, cette méthode échoue).

6.1 Discrétisation aux différences finies

Soit n ∈ N, n ≥ 1, on pose h := 1/(n + 1) ainsi que xi := ih pour 0 ≤ i ≤ n + 1,
0 = x0 ≤ x1 ≤ · · · ≤ xn+1 = 1.
Supposons qu’il existe une unique solution u ∈ C 4 ([0, 1]) au problème précédent. Par formules de Taylor,
u(xi+1 ) − u(xi ) u(xi ) − u(xi−1 )
− h2
00
∀i ∈ {1, . . . , n}, u (xi ) = −f (xi ) = h h + u(4) (xi + θi h), θi ∈ [−1, 1].
h 12
On simplifie le problème (en réalité on l’approche) en dimension finie, sous forme linéaire en le problème suivant.
Trouver U ∈ Rn tel que, ayant posé U0 = Un+1 = 0 (conditions aux limites), on ait
Ui+1 − 2Ui + Ui−1
∀i ∈ {1, . . . , n}, = −f (xi ),
h2
i. e. AU = F , où F = (f (xi ))i et  
2 −1 (0)
 .. .. 
1 −1 . . 
A=  .
h2  .. .. 
 . . −1
(0) −1 2
Définition 1.6.2 (Laplacien) La matrice
 
2 −1 (0)
 .. .. 
1 −1 . . 
An = 2  .
h  .. .. 
 . . −1
(0) −1 2
est égale à l’opposé de la matrice du laplacien.
Le problème à résoudre est (n + 1)2 An U = F , avec U ∈ Rn et F ∈ Rn donnée, pour approcher (S).

9
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

6.2 Propriétés hermitiennes de An

Proposition 1.6.3 La matrice An est symétrique définie positive.
Preuve. Soit U ∈ Rn , on impose U0 = Un+1 = 0. Par sommation d’Abel,
n
X n
X n
X n
X
hAn U, U i = Ui (−Ui−1 + 2Ui − Ui+1 ) = (Ui − Ui−1 )Ui−1 + (Ui − Ui−1 )Ui + Un2 = (Ui − Ui−1 )2 + U12 + Un2 .
| {z }
i=1 i=2 i=1 i=2
=−(Ui−1 +Ui )+(Ui −Ui−1 )

Si hAN U, U i = 0, alors U1 = Un = 0 et pour tout n ∈ {2, . . . , n}, Ui − Ui−1 = 0, donc U = 0. Par conséquent,
σ(An ) ⊂ R∗+ . En utilisant le rayon spectral, ρ(An ) ≤ kAn k∞ ≤ 4 (égalité si n ≥ 3). Ainsi, σ(An ) ⊂]0, 4]. 
Remarque 1.6.4 On peut montrer que An − 4In est symétrique défini négative, et donc σ(An ) ⊂]0, 4].
Proposition 1.6.5 (Spectre de An )
   
kπ
σ(An ) = 4 sin2 |1≤k≤n .
2(n + 1)

Preuve. 1ère méthode : Voir TD, par coïncidence algébrique (avec le problème continue).
2ème méthode : Cherchons U ∈ Rn et λ ∈]0, 4] tels que U 6= 0 et An U = λU . On résout

∀i ∈ {1, . . . , n}, −Ui−1 + 2Ui − Ui+1 − Ui+1 = λUi ,

avec U0 = Un+1 = 0. (Ui )1≤i≤n est solution d’une récurrence linéaire d’ordre 2. Le polynôme caractéristique de cette
récurrence est X 2 − (2 − λ)X + 1 = 0, de racines
1 p
r := [(2 − λ) + i λ(4 − λ)], et r̄.
2
Ainsi, pour tout j, Uj = αrj + βr̄j , avec α, β ∈ C. Comme |r|2 = rr̄ = 1, on a r = eiθ , avec θ ∈ R. Les conditions
U0 = Un+1 = 0 permettent d’obtenir α et β. En l’occurence, en (α, β),

U0 = α+β = 0
Un+1 = αr n+1
+ βr̄ n+1
= α[ei(n+1)θ − e−i(n+1)θ ].

Si la différence d’exponentielles est non nulle, alors α = β = 0, et U = 0, et λ ∈ / σ(An ). Si elle est nulle, on obtient
une droite vectorielle de solutions de la forme Uj = α(rj − r̄j ). Il suffit de trouver λ ∈]0, 4[ tel que r = eiθ vérifie

ei(n+1)θ = e−i(n+1)θ ,

ce qui est équivalent à θ = kπ/(n + 1), où k ∈ Z. Pour résoudre λ ∈]0, 4[ tel que r(λ) = eikπ/(n+1) , où k ∈ Z, on
observe que  
kπ 1
R(r(λ)) = cos = (2 − λ),
n+1 2
d’où    
kπ kπ
λ = 2 − 2 cos = 4 sin2 .
n+1 2(n + 1)
Parmis ces valeurs pour k ∈ Z, il y a exactement n valeur deux à deux distinctes. Ainsi,
       
2 kπ 2 kπ
σ(An ) ⊂ 4 sin | k ∈ Z = 4 sin |1≤k≤n .
2(n + 1) 2(n + 1)

L’inclusion réciproque est dans la synthèse de la preuve. 

On peut ainsi identifier explicitement
 
nπ
kAn k2 = ρ(An ) = 4 sin2 ,
2(n + 1)

et également
  −1
nπ
kA−1
n k2 = ρ(A−1
n ) = 4 sin2
.
2(n + 1)
De plus, cond2 (An ) = kAn k2 kA−1
n k2 .

10
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

6.3 Monotonie de la matrice An

Définition 1.6.6 (Positivité et monotonie)
• Un vecteur v ∈ Rn est dit positif si ses composantes le sont. On note alors v ≥ 0.
• Une matrice A ∈ Mn (R) est dite positive, noté A ≥ 0, si ses coefficients le sont.
• Une matrice A ∈ Mn (R) est dite monotone si A ∈ GLn (R) et A−1 ≥ 0.
Remarque 1.6.7 Il n’y a aucun lien entre la positivité des matrices symétriques et cette positivité ci. Par exemple, la
matrice  
1 2
2 1
est positive mais de spectre {−1, 3}.
Proposition 1.6.8 Soit A ∈ Mn (R). Il y a équivalence entre
• A est monotone,
• ∀v ∈ Rn , Av ≥ 0 =⇒ v ≥ 0.
Preuve. ( =⇒ ) : On suppose A ∈ GLn (R) et A−1 ≥ 0. Soit v ∈ Rn tel que Av ≥ 0, alors v = A−1 (Av), d’où
n
X
∀j ∈ {1, . . . , n}, vj = (A−1 )j,i (Av)i ≥ 0,
i=1
| {z } | {z }
≥0 ≥0

donc v ≥ 0.
( ⇐= ) : Soit v ∈ Rn tel que Av = 0, alors Av ≥ 0, donc v ≥ 0, et A(−v) ≥ 0, donc v ≤ 0. Donc v = 0, et
Ker(A) = {0} et A ∈ GLn (R). De plus, on pose ei le ième vecteur de la base canonique de Rn . ei ≥ 0 donc A−1 ei ≥ 0
par la propriété. Ce vecteur étant la ième colonne de A−1 , on a bien A−1 ≥ 0. 
Si A est monotone et x, y ∈ Rn sont tels que A(y − x) ≥ 0, alors y − x ≥ 0. On obtient un principe de comparaison
des solutions de Ax = f .
Proposition 1.6.9 Pour tout n ∈ N∗ , An est monotone.
Preuve. Soit v ∈ Rn tel que An v ≥ 0. Montrons que v ≥ 0. Les coefficients de v vérifient

∀i ∈ {1, . . . , n}, −vi−1 + 2vi − vi+1 ≥ 0,

après avoir posé v0 = vn+1 = 0. La première équation est 2v1 − v2 ≥ 0. Posons k ∈ {1, . . . , n} tel que vk = min1≤i≤n vi .
Il suffit de montrer que vk ≥ 0.
Si k = 1, alors
2v1 − v2 ≥ 0 =⇒ v1 = vk ≥ v2 − v1 ≥ 0.
|{z} |{z}
équation définition de k

Si k = n, alors
2vn − vn−1 ≥ 0 =⇒ vn = vk ≥ vn−1 − vn ≥ 0.
Si k ∈ {2, . . . , n − 1}, alors la k ème
équation (Avk )k = −vk+1 + 2vk − vk−1 ≥ 0, de sorte que

(−vk+1 + vk ) + (vk − vk−1 ) ≥ 0.

| {z } | {z }
≤0 ≤0

Ainsi, vk = vk−1 = vk+1 . On peut, par récurrence sur l’indice, montrer que les vi sont tous égaux 
Conséquence : Pour n ≥ 3, on peut montrer que kAn k∞ vaut (n + 1)2 /8 si n est impair et n(n + 2)/8 sinon.
Preuve. Observons que la fonction
[0, 1] −→ R
ϕ: 1
x 7−→ x(1 − x)
2
est solution de −ϕ00 = 1 sur ]0, 1[, et que ϕ(0) = ϕ(1) = 0. Posons pour n ≥ 3, φ ∈ Rn tel que
 
i 1
φi = ϕ = i(n + 1 − i).
n+1 2(n + 1)2

On peut calculer
1
(n + 1)2 An φ = [2i(n + 1 − i) − (i − 1)(n + 1 − i + 1) − (i + 1)(n + 1 − i − 1)]1≤i≤n = (1)1≤i≤n .
2

11
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Soit F ∈ Rn tel que

∀i ∈ {1, . . . , n}, kF k∞ ≤ Fi ≤ kF k∞ ,
alors    
1 1
 ..   .. 
0 ≤ v1 := kF k∞  .  − F, et v2 := F + kF k∞  .  ≥ 0.
1 1
Or (n + 1)2 An est monotone, donc

(n + 1)2 An v ≥ 0 =⇒ v ≥ 0 ou v ≥ 0 =⇒ ((n + 1)2 An )−1 V ≥ 0.

On en déduit que
≤ ((n + 1)2 An )−1 v1

0
,
0 ≤ ((n + 1)2 An )−1 v2
mais  
1
2 −1  .. 
((n + 1) An )  .  = φ,
1
d’où, en comparant composante par composante,

−kF k∞ φ ≤ ((n + 1)2 An )−1 F ≤ kF k∞ φ.

| {z }
=U

Ainsi,
kU k∞ ≤ kF k∞ kφk∞ .
On peut interpréter cela comme

A−1
n F k∞ (n + 1)2 k((n + 1)2 An )−1 F k∞
kA−1
n k = sup = sup ≤ (n + 1)2 kφk∞ .
F 6=0 kF k∞ F 6=0 kF k ∞

Il y a égalité pour F = (1, . . . , 1)T . Reste à calculer kφk∞ . 

6.4 Convergence de la méthode (en norme uniforme)

Soit f ∈ C 2 ([0, 1]) et u l’unique solution du système

−u00

= f
(S) :
u(0) = u(1) = 0

En fait, u ∈ C 4 ([0, 1]). Pour n ∈ N∗ , on considère les objets suivants.

• U ex ∈ Rn le vecteur de composantes Uiex = u(i/(n+1)), pour 1 ≤ i ≤ n. En particulier, il vérifie (n+1)2 An U ex =
F + ε, avec F = (f (i/(n + 1)))1≤i≤n , et ε ∈ Rn , appelé erreur de coïncidence. Par développement de Taylor,
 
1 1 (4) i θi
∀i ∈ {1, . . . , N }, ∃θi ∈ [−1, 1], εi = u + .
(n + 1)2 12 n+1 n+1
Donc
1
kεk∞ ≤ max |f (2) (x)|.
12(n + 1)2 x∈[0,1]
• U (n) ∈ Rn la solution du problème approché (n + 1)2 An U (n) = F .
Question : kU ex − U (n) k∞ →(n→+∞) 0 ? On a

(n + 1)2 An (U ex − U (n) ) = F + ε − F = ε,

donc
kU ex − U (n) k∞ ≤ (n + 1)−2 kA−1 k∞ kεk∞ .
| n{z } | {z }
∼n2 /8 ∼n−2 kf (2) k ∞ /12

Ainsi,
kU ex − U (n) k∞ ≤ Cn−2 kf (2) k∞ −−−−→ 0,
n→∞
n
pour C ∈ R une constante indépendante de n et de f .

12
ANU CHAPITRE 1. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Remarque 1.6.10
n
(n)
X
kU ex − U (n) k22 = |Uiex − Ui |2 ≤ nkuex − u(n) k1∞ ,
i=1

donc √
kU ex − U (n) k2 ≤ nCn−2 kf (2) k∞ ≤ Cn−3/2 kf (2) k∞ −−−−−→ 0.
n→+∞

13
Chapitre 2

Méthodes directes

1 Remarques préliminaires
Soient A ∈ GLn (K) et b ∈ Kn . Pour résoudre Ax = b par formule de Cramer, et calculs de déterminants de taille
n, la complexité de calcul est en O(nn!). Sur un supercalculateur à 200 petaflops (= 200.1015 opérations flottants
par seconde), pour n = 50, il faut 2.1061 années. Pour n = 12, il faut environ une année. En réalité, on peut trouver
des algorithmes en O(n3 ) opérations. Pour n = 10000, votre machine peut résoudre en quelques secondes. Le cas
particulier des matrices triangulaires mène à deux algorithmes :
• l’algorithme de descente, pour A ∈ Tnl (K), en O(n2 ) opérations,
• l’algorithme de remontée, pour A ∈ Tnu (K), de même complexité.

2 Factorisation LU
Rappel : Opérations du pivot de Gauss.
• Opérer sur les lignes équivaut à multiplier à gauche par A.
• Opérer sur les colonnes équivaut à multiplier à droite par A.
• La matrice D(λ, i) = In + (λ − 1)Ei,i , où λ ∈ / {0, 1}, est la matrice de dilatation de la ième ligne par λ, Li ← λLi .
• La matrice P (i, j) = In − Ei,i − Ej,j + Ei,j + Ej,i est la matrice de permutation, qui permute les ième et j ème
lignes, (Li , Lj ) ← (Lj , Li ).
• La matrice T (λ, i, j) = In + λEi,j , où i 6= j et λ 6= 0, est la matrice de transvection, Li ← Li + λLj .
Remarque 2.2.1 Le calcul du rang, déterminant, noyau, image peuvent se faire à partir de ces opérations en examinant
les invariants de ces opérations. Par exemple,

Ker A = Ker(T (λ, i, j)A), Im A = Im(AT (λ, i, j)).

Pour résoudre Ax = b, les opérations du pivot de Gauss peuvent être « stockées »dans une factorisation LU de A.
De même, inverser A revient à résoudre n systèmes linéaires Ax = ei , où ei est le ième vecteur de la base canonique.
Définition 2.2.2 (Mineurs principaux dominants) Soit A ∈ Mn (K), les mineurs principaux dominants de A
sont les det((Ai,j )1≤i,j≤k ), pour 1 ≤ k ≤ n.
Remarque 2.2.3
• Toute matrice dont les mineurs principaux dominants sont non nuls est inversible.
• La réciproque est fausse. Par exemple,  
0 1
1 1
est inversible, mais le mineurs principaux dominants de taille 1 est nul.
• Toute matrice symétrique définie positive réelle a ses mineurs principaux dominants tous non nuls. C’est une
conséquence de la réduction des formes bilinéaires. En réalité, si A est symétrique, il y a équivalence entre A
symétrique définie positive et ses mineurs principaux dominants sont strictement positifs.
Lemme 2.2.4 Soit A ∈ Mn (K) dont tous les mineurs principaux dominants sont non nuls. Alors pour tous λ ∈ K,
pour tous i < j, T (λ, i, j)A a aussi ses mineurs principaux dominants non nuls.

14
ANU CHAPITRE 2. MÉTHODES DIRECTES

Preuve. Soit k ∈ {1, . . . , n},

    
Tk 0 Ak Bk Tk Ak ∗
T (λ, i, j)A = = .
T̃ Tn−k Ck Dk ∗ ∗

Le mineurs principaux dominants de taille k est det(Tk Ak ), avec


Ik si k ≥ i ou k ≥ j
Tk = ,
Ik + λEi,j sinon

donc det Tk = 1, d’où det(Tk Ak ) = det Ak 6= 0. 

Théorème 2.2.5 Soit A ∈ Mn (K) dont les mineurs principaux dominants sont tous non nuls. Alors il existe un
unique couple (L, U ) ∈ Tnl (K) × Tnu (K) tel que

∀i ∈ {1, . . . , n}, Li,i = 1
.
A = LU

Preuve. Unicité :
A = L1 U1 = L2 U2 ⇐⇒ L−1 L1 = U2 U −1 .
| 2{z } | {z1 }
∈Tnl (K) ∈Tnu (K)

Les coefficients diagonaux de L−1 −1

2 L1 étant tous égaux à 1, on a L2 L1 = In , d’où L1 = L2 et U1 = U2 .
Existence : Par l’algorithme effectif de construction. Soit k ∈ {1, . . . , n − 1}, supposons déterminées des matrices
T1 , . . . , Tk−1 chacune produit de transvections inférieures telles que
 (k) (k) 
a1,1 (ai,j
 .. 
 0
 . 

 . .. 
 .. . (k)
ak−1,k−1 
Tk−1 . . . T1 A =  . .
 
 . (k)
 . 0 ak,k


 .. .. ..
 

 . 0 . . 
(k)
0 ... 0 an,k

(k) (k)
Par le lemme, pour tous j ∈ {1, . . . , k}, aj,j 6= 0. En particulier, le pivot ak,k 6= 0. On peut faire l’opération Li ←
(k+1) (k) (k) (k)
Li − xi,k Lk , avec xi,k := ai,k /ak,k pour tous i ≥ k + 1. On pose
n
(k)
Y
Tk := T (−xi,k , i, k),
i=k+1

alors  (k) (k) 

a1,1 (ai,j
 0 
 
 . ..
 ..

 . 

 . 
Tk . . . T1 =  ..
(k)
 ak,k .

 . 
 .
 . 0


 . .. 
 . 
 . . 
(k+1)
0 ... 0 (ai, j )
On conclut par récurrence sur k. Comment obtenir A = LU ? On pose U := Tn−1 . . . T1 A ∈ Tnu (K), ainsi que
L := (Tn−1 . . . T1 )−1 ∈ Tnl (K) de diagonale 1. On en déduit l’existence. U est construite dans la récurrence. L s’obtient
simplement :
1 0
 
..
.
 
 
L= 
..
.

 . 
(j)
(xi,j ) 1

15
ANU CHAPITRE 2. MÉTHODES DIRECTES

Pourquoi ?
n
(k)
Y
Tk = (In − xi,n Ei,k ),
i=k+1
n n
(k) (k)
Y X
Tk−1 = (Ik + xi,n Ei,k ) = In + xi,k Ei,k .
i=k+1 i=k+1

Par exemple, Tk−1 Tk+1

−1
fait intervenir des produits Ei,k Ej,k+1 avec i ≥ k + 1 et j ≥ k + 2,
  
X (k) X (k+1)
Tk−1 Tk+1
−1
= I + xi,k Ei,k  I + xj,k+1 Ej,k+1 
i≥k+1 j≥k+2
(k) (k+1)
X X
=I+ xi,k Ei,k + xj,k+1 Ej,k+1 + 0.
i≥k+1 j≥k+1

et ce produit est nul. On obtient ainsi

n−1 n
(j)
X X
L = T1−1 . . . Tn−1
−1
= In + xi,j Ei,j .
j=1 i=j+1


Conséquence : Pour résoudre un ou plusieurs systèmes Ax = b (pour différents b, par exemple pour A−1 ), on peut
• précalculer A = LU , en O(2n3 /2) opérations,
• résoudre Ly = b et U x = y, chacun en O(n2 ) opérations.
On a « stocké »les opérations du pivot du second memebre b dans la matrice L.
Concernant le choix des pivots, si l’un des mineurs principaux dominants est nul, supposons le k ème le plus petit
(k) (k)
non nul, alors l’algorithme s’arrêt à la k ème étape car ak,k = 0. Mais si A est inversible, il existe ai,k 6= 0 pour i ≥ k + 1
(k)
(sans que rg A ≤ n − 1). On pourrait permuter les lignes k et la i pour utiliser ai,k comme pivot. Cela donne une
factorisation de A = P LU , où P est une matrice de permutation.
(k)
Si un pivot est très petit, une difficulté numérique se pose sur l’arithmétique approchée. L’erreur de calcul sur ai,k
(k)
par |ak,k | >> 1. On peut alors établir une stratégie de pivot en permutant pour limiter cet effet.
(k)
• Pivot partiel : Soit avec la ligne i de sorte que maxi≥k |ai,k | est réalisé sur la diagonale.
(k)
• Pivot total : On suppose |ai,j | maximal pour i, j ≥ k. On permute lignes et colonnes pour utiliser cette valeur
comme pivot, A = P LU Q.
On peut citer l’exemple des matrices bandes / creuses, dans lesquelles seule un groupe de sous-diagonales est non
nul. On appelle demi-largeur de bande l’entier p tel que p sous-diagonales soient non nulles de chaque côté de la
diagonale. La complexité de la factorisation LU est alors en O(np2 ) opérations, et les matrices L et U ont le même
profil (seulement p sous-diagonales non nulles).

3 Décomposition de Cholesky
Théorème 2.3.1 (Décomposition de Cholesky) Soit A ∈ Sn++ (R), il existe une unique B ∈ Tnl (R) telle que

∀i ∈ {1, . . . , n}, Bi,i > 0,
.
A = BB T .

Preuve. Existence : A admet une unique factorisation LU = A, avec Li,i = 1 pour tous 1 ≤ i ≤ n.
 
∗ ∗
U =
 ..  = DŨ ,

.
0 ∗

avec D = diag((ui,i )1≤i≤n ), Ũ ∈ Tnu (R), Ũi,i = 1 pour tous 1 ≤ i ≤ n. Comme A = LU = AT , et que AT admet une
factorisation LU, qui est U T LT . On en déduit que Ũ T DT LT = LDŨ . Par unicité, L = Ũ T , Ũ = LT et D = DT , et
donc A = LDLT . Pour x 6= 0, comme A est symétrique définie positive, on a

0 < hAx, xi = hLDLT x, xi = hDLT x, LT xi,

16
ANU CHAPITRE 2. MÉTHODES DIRECTES

donc
∀y 6= 0, hDy, yi > 0,
√ √
ce qui montre que D ∈ Dn++ (R). On peut donc écrire A = LDŨ = L D DŨ . On pose B := L diag(( Di,i )1≤i≤n ) ∈
p

Tnl (R), avec Bi,i = Li,i Di,i > 0. Ainsi, A = BB T .

p

Unicité : laissé en exercice. 

Algorithme : A = BB T , équations d’inconnues (Bi,j ).
j−1
X
∀j ∈ {1, . . . , n}, b2j,j + b2j,k = aj,j .
k=1

On pose v
u n
X
u
bj,j := taj,j + b2j,k ,
k=1

et donc
j+1
X
ai,j − bj,k bi,k
k=1
∀i ≥ j + 1, bi,j = .
bj,j
Cette méthode est en O(n3 /3) opérations.

17
Chapitre 3

Méthodes itératives

1 Généralités
Soit A ∈ GLn (K), on cherche à résoudre Ax = b. On choisit deux matrices M, N ∈ Mn (K) telles que

A=M −N
.
M est facile à inverser (algorithmiquement peu coûteuse à inverser).

On définit alors pour un x0 ∈ Kn donné la suite récurrente

∀k ≥ 1, M xk+1 = N xk + b.

Définition 3.1.1 (Méthode itérative convergente) On dit qu’une méthode itérative est convergente si pour tout
x0 ∈ Kn , la suite (xk )k converge.
Remarque 3.1.2 A étant inversible, la seule limite possible est la solution x = A−1 b.
Définition 3.1.3 (Vocabulaire)
• Résidu à l’étape k : rk := b − Axk .
• Erreur à l’étape k : ek := x − xk .
On a
xk −−−−−→ x ⇐⇒ ek −−−−−→ 0 ⇐⇒ rk −−−−−→ 0.
k→+∞ k→+∞ k→+∞

Remarque 3.1.4 ek = x − xk = A(b − Axk ) = A−1 rk , donc kek k ≤ kA−1 kkrk k. Cela peut servir à définir un critère
d’arrêt à partir de krk k et de kA−1 k.
Théorème 3.1.5 La méthode itérative pour A = M − N est convergente si et seulement si ρ(M −1 N ) < 1.
Preuve. On a
∀k ≥ 0, xk+1 = M −1 N xk + M −1 b

,
x = M −1 N x + M −1 b
donc pour tout k ∈ N, ek+1 = M −1 N ek , soit ek = (M −1 N )k e0 . On souhaite que pour tout e0 ∈ Kn , ek →(k→+∞) 0,
donc que (M −1 N )k →(k→+∞) 0. C’est équivalent au fait que ρ(M −1 N ) < 1. 
Remarque 3.1.6 La vitesse de convergence est au plus géométrique, et est liée à la quantité ρ(M −1 N ) :

∀ε > 0, ∃k · k norme sur Kn , kM −1 N k ≤ ρ(M −1 N ) + ε.

Ainsi,
∀k ∈ N, k(M −1 N )k kke0 k ≤ (ρ(M −1 N ) + ε)k ke0 k,
avec pour ε petit, ρ(M −1 N ) + ε < 1. Pour une norme donnée, par exemple k · k2 sur Kn , par équivalence des normes
en dimension finie, il existe Cε > 0 tel que

∀k ∈ N, kek k2 ≤ C(ρ(M −1 N ) + ε)k ke0 k2 .

On appelle taux de convergence la quantité − log(ρ(M −1 N )). La convergence est d’autant plus rapide (pour le x0 le
pire) que le taux est grand.

18
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

2 Méthodes usuelles
2.1 Méthode de Richardson
Soit α ∈ R∗ , on pose M := α−1 In et N := α−1 In − A. La suite récurrence est

xk+1 = M −1 N xk + M −1 b = xk + α(b − Axk ),

et M −1 N = In − αA. Par conséquent,

σ(M −1 N ) = {1 − αλ, λ ∈ σ(A)}.

Supposons que A ∈ Sn++ (R), alors σ(A) ⊂ R∗+ . On note 0 < λ1 ≤ · · · ≤ λn les valeurs propres de A. On a

−1 |1 − αλn | si α ∈]0,
/ +α∗ [
ρ(M N ) = ,
|1 − αλ1 | si α ∈]0, +α∗ ],

avec α∗ := 2/(λ1 + λn ). En particulier,

   
−1 2 2
ρ(M N ) < 1 ⇐⇒ α ∈ 0, = 0, .
λn ρ(A)
Mieux, le rayon spectral est minimal lorsque α = α∗ . On a alors
λn − λ1 cond2 (A) − 1
ρ(M −1 N ) = = .
λn + λ1 cond2 (A) + 1
Si cond2 (A) >> 1, cette quantité approche 1 par valeurs inférieures. Si cond2 (A) ' 1, elle approche 0.

2.2 Méthode de Jacobi

On choisit M := D := diag(A). Il est nécéssaire de supposer D ∈ GLn (R). La matrice d’itération est M −1 N =
D (D − A) = In − D−1 A. En pratique, l’algorithme « en place » prend la forme :
−1

Soit x(0) ∈ Rn ,  
(k+1) 1  X (k) 
∀k ≥ 0, ∀i ∈ {1, . . . , n}, xi = bi − ai,j xj .
ai,i
j6=i

Idée : x vérifie
n
X
∀i ∈ {1, . . . , n}, ai,i xi + ai,j xj = bi .
|{z} |{z}
inconnu j=1
supposé connu
j6=i

2.3 Méthode de Gauss-Seidel

Soit A ∈ Mn (R), on note D ∈ Mn (R) sa diagonale, −E ∈ Tnl (R) sa sous-diagonale et −F ∈ Tnu (R) sa sur-
diagonale, de sorte que A = D − E − F . On pose

M := D − E
N := F.

M est inversible si et seulement si D ∈ GLn (R).

Algorithme : On pose pour tous k ≥ 0, i ∈ {1, . . . , n},
 
i+1 n
(k+1) 1  X (k+1)
X (k) 
xi = bi − ai,j xj − ai,j xj .
ai,i j=1 j=i+1

2.4 Méthode de relaxation (SOR)

Soit ω ∈ R∗+ , avec les notations précédentes, on pose

1
 M := D−E

ω
1−ω
 N
 := D + F.
ω
M est inversible si et seulement si D ∈ GLn (R).

19
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

3 Critères explicites de convergence

3.1 Matrices à diagonale dominante
Définition 3.3.1 (Matrice à diagonale strictement dominante) Soit A ∈ Mn (C), on dit qu’elle est à diagonale
strictement dominante si
n
X
∀i ∈ {1, . . . , n}, |ai,i | > |ai,j |.
j=1
j6=i

Remarque 3.3.2 Toute telle matrice a une diagonale inversible. Les méthodes précédentes peuvent donc s’appliquer,
i. e. on peut bien définir la suite récurrente (xk )k .
Proposition 3.3.3 Si A ∈ Mn (C) est à diagonale strictement dominante, alors A ∈ GLn (C).
Preuve. Soit x ∈ Ker A, alors
n
X
∀i ∈ {1, . . . , n}, ai,j xj = 0,
j=1

donc X
∀i ∈ {1, . . . , n}, ai,i xi = − ai,j xj .
j6=i

Soit i ∈ {1, . . . , n}, par inégalité triangulaire,

 
X X
|ai,i ||xi | ≤ |ai,j ||xj | ≤  |ai,j | kxk∞ .
j6=i j6=i

Or il existe i ∈ {1, . . . , n} tel que |xi | = kxk∞ , et donc

 
X
 |ai,j | − |ai,i | kxk∞ ≥ 0,
j6=i
| {z }
<0

d’où x = 0. 
A-t-on convergence des précédentes méthodes pour de telles matrices ?
Théorème 3.3.4 Soit A ∈ Mn (C) à diagonale strictement dominante. Le méthodes de Jacobi, Gauss-Seidel et
relaxation convergent pour ω ∈]0, 1].
Preuve. On se contentera de donner la preuve pour Jacobi. On a avec les notations précédentes, M −1 N = In −D−1 A.
Soit λ ∈ σ(M −1 N ) tel que |λ| = ρ(M −1 N ) et x ∈ Cn un vecteur propre associé. Alors M −1 N x = λx, donc N x = λM x,
soit (D − A)x = λDx. On considère l’indice i ∈ {1, . . . , n} tel que |xi | = kxk∞ 6= 0. On a
X
ai,j xj
X j6=i
− ai,j xj = λai,i xi , donc λ = − .
ai,i xi
j6=i

Ainsi, X X
|ai,j ||xj | |ai,j |
j6=i j6=i kxk∞
|λ| ≤ ≤ .
|ai,i ||xi | ai,i |xi |
| {z } |{z}
<1 par hypothèse =kxk∞

Donc ρ(M −1 N ) < 1 et la méthode converge. 

3.2 Matrices Hn++ (C)

Théorème 3.3.5 Soit A ∈ Hn++ (C), on écrit A = M − N avec M ∈ GLn (C). Alors M ∗ + N est hermitienne.
Si de plus M ∗ + N ∈ Hn++ (C), alors ρ(M −1 N ) < 1 (et la méthode itérative converge).
Preuve. (M ∗ + N )∗ = M + N ∗ = (A + N ) + N ∗ = (A∗ + N ∗ ) + N = M ∗ + N , donc cette matrice est hermitienne.

20
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

Supposons que M ∗ + N ∈ Hn++ (C). Il suffit de trouver une norme subordonnée sur Mn (C) tele que kM −1 N k < 1.
On va utiliser la norme induite par A, induite par le produit scalaire

∀x, y ∈ Cn , hx, yiA = hAx, yi.

Considérons x ∈ Cn tel que kxkA = 1 et kM −1 N xk1 = kM −1 N kA . On a

kM −1 N xk21 = hAM −1 N x, M −1 N xi
= hAM −1 (M − A)x, M −1 (M − A)xi
= hAx − Ax, x − yi (y := M −1 Ax)
= hAx, xi − hAy, xi −h |{z}
Ax , yi + hAy, yi
| {z } | {z }
=kxkA =1 =hy,Axi =M y
=hy,M yi
| {z }
=hM ∗ y,yi =−hN y,yi

∗
= 1 − h(M + N )y, yi (y 6= 0 car M −1 A ∈ GLn (C), x 6= 0)
< 1.


Exemple 3.3.6 On considère la matrice An de l’opposé du Laplacien. On sait que An ∈ Sn++ (R). On lui applique la
méthode de Jacobi. Posons M := 2In , alors M ∗ + N = −An + 4In , mais σ(An ) ⊂]0, 4[, d’où σ(M ∗ + N ) ⊂]0, 4[⊂ R∗+ .
Par conséquent, M ∗ + N ∈ Sn++ (R) et la méthode de Jacobi converge.
Exemple 3.3.7 Pour la méthode de Gauss-Seidel, M ∗ + N = 2In et la méthode converge aussi.

4 Résolution de systèmes linéaires au sens des moindres carrés

SOient n, p ∈ N∗ , A ∈ Mn,p (K) et b ∈ Kn . On cherche un vecteur x ∈ Kp tel que Ax = b. On dit que le système
linéaire Ax = b est
1. carré si n = p,
2. surdéterminé si n > p,
3. sous-déterminé si n < p,
4. de rang maximal si rg A = min(, np),
5. incompatible si b ∈/ Im A,
6. compatible si b ∈ Im A.
Si le système est compatible, alors il existe des solutions. Si rg A = p, il existe une unique solution et si rg A < p, il
existe une infinité de solutions. En particulier, s’il est surdéterminé, de rang maximal et compatible, alors il existe une
unique solution. S’il est sous-déterminé de rang maximal, alors il est compatible et il existe une infinité de solutions.
Mais que se passe-t-il si le système est incompatible ?
Définition 3.4.1 (Résolution au sens des moindres carrés) Résoudre le système Ax = b au sens des moindres
carrés, c’est trouver un vecteur x ∈ Kp minimisant la quantité kAx − bk22 .

4.1 Existence et unicité

Théorème 3.4.2 L’ensemble des solutions du système Ax = b au sens des moindres carrés est

{x ∈ Kp | Ax = p(b)},

où l’application p est la projection orthogonale sur Im A.

Preuve. Comme p(b) ∈ Im A, la caractérisation du projeté donne

∀w ∈ Im A, hp(b) − b, wi = 0.

Soit v ∈ Im A, d’après le théorème de Pythagore,

kv − bk2 = kv − p(b)k2 + kp(b) − bk2 ≥ kp(b) − bk2 ,

avec égalité si et seulement si v = p(b). De plus, il existe x0 ∈ Kp tel que Ax0 = p(b). L’unicité est garantie si et
seulement si Ker A = {0}. 

21
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

4.2 Équation normale

Lemme 3.4.3 Les solutions du système Ax = b au sens des moindres carrés sont exactement les solutions de
l’équation A∗ Ax = A∗ b.
Preuve. D’après le lemme précédent, il faut montrer que Γ1 := {x ∈ Kp | Ax = p(b)} est égal à Γ2 := {x ∈ Kp |
A∗ Ax = A∗ b}.
Soit x ∈ Γ1 , alors A∗ Ax = A∗ p(b). Montrons que A∗ p(b) = A∗ b, i. e. que p(b) − b ∈ Ker A∗ . On sait que
Ker A∗ = (Im A)⊥ , or p(b) − b ∈ (Im A)⊥ , d’où A∗ p(b) = A∗ b, et x ∈ Γ2 .
Soit x ∈ Γ2 . Comme p(b) − b ∈ Ker A∗ , on a

A∗ Ax = A∗ (b − p(b) + p(b)) = A∗ p(b),

d’où Ax − p(b) ∈ Ker A∗ . Mais Ker A∗ = (Im A)⊥ , et Ax − p(b) ∈ Im A. Par conséquent, Ax − p(b) = 0, donc x ∈ Γ1 .
Ainsi, Γ1 = Γ2 . 
Remarque 3.4.4 L’équation normale A∗ Ax = A∗ b fait intervenir une matrice A∗ A qui est hermitienne et positive. De
plus, son noyau est égal à celui de A. Ainsi, le système carré AA∗ x = b est inversible et il existe une unique solution
pour les moindres carrés.

4.3 Exemple de la régression linéaire

Soit (xi , yi )1≤i≤n une famille de points de R2 . On cherche un couple (a, b) ∈ R2 tel que

∀i ∈ {1, . . . , n}, axi + b = yi .

Matriciellement, le problème s’écrit sous la forme A(a, b) = b, avec

   
x1 1 y1
 .. .
.  .. 
A := .  , et b :=  . .

 .
xn 1 yn

L’équation normale fait intervenir la matrice et le vecteur

n
X n
X  n
X 
 x2i xi   xi yi 
∗  i=1 i=1 ∗ i=1
A A = X , et A b =  X .
 
n X n n
   
xi 1 yi
i=1 i=1 i=1

De plus, on a Ker A = Ker A∗ A = {0} dès que n ≥ 2, et qu’au moins deux des réels xi sont distincts. Après résolution
de l’équation normale, on obtient que
 ! n !
Xn Xn X Xn

n x y − x y (xi − x̄)(yi − ȳ)

i i i i





 a = i=1 i=1 i=1 i=1
!2 = n
n n X ,
(xi − x̄)2
X X
2 2



 n x i − x i

 i=1 i=1 i=1

b = ȳ − ax̄


avec
n n
1X 2 1X 2
x̄ := x , et ȳ = y .
n i=1 i n i=1 i

4.4 Factorisation QR par les matrices de Householder

On veut utiliser de symétries orthogonales pour produire une factorisation QR.
Remarque 3.4.5 Si Q ∈ Hn (C) ∩ Un (C), alors Q−1 = Q∗ = Q.
Définition 3.4.6 (Matrice de Householder) On appelle matrice de Householder associée à un vecteur v ∈ Rn la
matrice 
 In si v=0
H(v) := vv T
 In − 2 sinon .
vT v

22
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

Proposition 3.4.7 Soit v ∈ Rn , alors

1. la matrice H(v) est symétrique et orthogonale,
2. si v 6= 0, alors H(v) est la matrice de la symétrie orthogonale sur v ⊥ parallèlement à Vect(v),
3. pour tout e ∈ Rn tel que e∗ e = 1, on a

H(v + kvke)v = −kvke, et H(v − kvke)v = kvke.

On obtient ainsi un nouvel algorithme d’élimination. Soit A ∈ Mn,p (R). On note a(1) ∈ Rn la première colonne
de A et e1 ∈ Rn le premier vecteur de la base canonique de Rn . On pose H (1) := H(a(1) − ka(1) ke1 ). Alors
 (1) 
ka k
H (1) e(1) = ka(1) ke1 , et H (1) A =  ... ∗ ,
 

où cette dernière matrice est du même rang que A. Supposons qu’à l’étape k ∈ {1, . . . , n}, on a
 
T ∗
|
H (k−1) . . . H (1) A = 
 
,
(0) a(k) ∗ 
|

où T est triangulaire supérieure et a(k) ∈ Rn−k+1 est non nul si rg A > k. On pose alors
 
I 0
H (k) = k−1 ,
0 H(a(k) − ka(k) k(1, 0, . . . , 0)
de sorte que  
T ∗

 ka(k+1) k 


 0 .

 .. 
(0) . ∗
0
Remarque 3.4.8
• Si n = p et A est inversible, on obtient A = QR avec R ∈ Tnu (R) et Q := H (1) . . . H (n) ∈ On (R). On a également
R ∈ GLn (R).
• Si n > p et rg A = p, on obtient A = QR, avec Q ∈ On (R) et R ∈ Mn,p (R). La matrice R peut être mise sous
la forme  
R1
R= ,
0n−p


avec R1 ∈ Tpu (R) ∩ GLp (R), et la matrice Q peut se décomposer sous la forme Q = Q1 Q2 . De cette manière,
on a A = Q1 R1 . De plus, pour tout x ∈ Rp , on a

kAx − bk2 = kQRx − bk2 = kRx − Q∗ bk2 = kR1 x − Q∗1 bk2 + kQ∗2 bk2 .

On trouve donc la solution au sens des moindres carrés qui est x = R1−1 Q∗1 b, et le résidu est kAx − bk = kQ∗2 bk.
• Si n > p et r := rg A < p, on obtient Q = QR, où la matrice R se met sous la forme
 
R1 R2
R= ,
0 0


avec R1 ∈ Tru (R) et R2 ∈ Mn,p−r (R). De plus, la matrice Q se met sous la forme Q = Q1 Q2 , avec
Q1 ∈ Mn,r (R). Alors pour tout x = (x1 , x2 ) ∈ Rr × Rp−r , on a

kRx − Q∗ bk2 = kR1 x1 + R2 x2 − Q∗1 bk2 = kQ∗2 bk2 .

Cette quantité est minimale pour x1 = R1−1 Q∗1 b − R1−1 R2 x2 , i. e. pour

 −1 ∗  
−R1−1 R2 x2

R1 Q1 b
x= + , avec x2 ∈ Rp−r .
0 x2
| {z }
∈Ker A

Un défaut d’unicité appraît par la présence du vecteur x2 ∈ Rp−r en paramètre, qui est traité en choissant la
solution de norme minimale.

23
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

4.5 Décomposition en valeurs singulières

Définition
√ 3.4.9 (Valeurs singulières) Soit A ∈ Mn,p (C). Les valeurs singulières de A sont les réels positifs
σi := λi , où les λi sont les valeurs propres réelles non nulles de A∗ A.
Remarque 3.4.10 Comme les matrices A et A∗ A ont même noyau, elles ont le même rang, donc la matrice A∗ A admet
exactement r := rg A valeurs singulières, comptées avec leurs multiplicités.
Théorème 3.4.11 Soit A ∈ Mn,p (C) de rang r ≤ min(n, p). Alors A ademt exactement r valeurs singulières
comptées avec leurs multiplicités, et il existe U ∈ Un (C) et V ∈ Up (C) telles que A = U ΣV ∗ , avec

Σ := diag(σ1 , . . . , σr , 0) ∈ Mn,p (R),

où les réels σi sont les valeurs singulières de A.

Définition 3.4.12 (Pseudo-inverse) On appelle pseudo-inverse d’une matrice A ∈ Mn,p (C) la matrice

A+ := V Σ+ U ∗ , avec Σ+ := diag(σ1−1 , . . . , σr−1 , 0) ∈ Mp,n (R),

en reprenant les notations du théorème précédent.

Remarque 3.4.13 La pseudo-inverse ne dépend pas du choix de la décomposition en valeurs singulières.
Appliquons ceci à la méthode des moindres carrés. On écrit A = U ΣV ∗ la décomposition en valeurs singulières de
A. Soit x ∈ Cp . Alors
kAx − bk = kU ∗ Ax − U ∗ bk = kΣV ∗ x − U ∗ bk.
On pose c := U ∗ b ∈ Cn , et y := V ∗ x ∈ Cp . Alors
r
X n
X
kΣy − ck2 = |σi yi − ci |2 + |ci |2 ,
i=1 i=r+1

et ce terme est minimal si yi = σi−1 ci pour tout i ∈ {1, . . . , r}, i. e. y = Σ+ c+y0 , avec y0 := (0, . . . , 0, yr+1 , . . . , yp ) ∈ Cp .
D’où x = A+ b + V y0 . Parmis ces solutions, le vecteur x est de norme minimale si et seulement si le vecteur y l’est,
avec
X r n
X
2 −1 2
kyk <= |σi ci | + |ci |2 ,
i=1 i=1

i. e. si et seulement si x = A+ b.

5 Méthodes variationelles
5.1 Principe
5.2 Algorithme du gradient à pas fixe
5.3 Algorithme du gradient à pas optimal
On considère une suite (xj )j de Rn définie par

∀j ∈ N, xj+1 = xj − αj ∇f (xj ), avec αj := argmin f (xj − α∇f (xj )).

α∈R

Petit calcul : Soit A ∈ Sn++ (R). On va poser A1/2 ∈ Sn++ (R) vérifiant (A1/2 )2 = A. Ainsi,

∀y ∈ Rn , kyk2A = hAy, yi = hA1/2 y, A1/2 yi = kA1/2 yk22 ,

donc
1 1 1
∀y ∈ Rn , f (y) = (hAy, yi − 2hAx, yi) = (hA(y − x), y − xi − hAx, xi) = (ky − xk2A − kxk2A ),
2 2 2
où Ax = b.

24
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

Proposition 3.5.1 Pour tout espace affine E de Rn ,

argmin f = argmin ky − xk2A
E y∈E

est le projeté A-orthogonal de x sur E.

Preuve. On a
k−1
∀j ∈ N, kxj+1 − xkA ≤ kxj − xkA ,
k+1
avec K := cond2 A. Notons pour j ∈ N, ej := xh − x. Alors
kej+1 kA = min k xj − α(Axj − b) − x kA ≤ kxj − αR (Axj − b) − xkA ,
α∈R | {z }
:=e(R)
j+1

où αR est la constante optimale dans la méthode itérative de Richardson, αR := 2/(λmin + λmax ). Et

(R) (R)
kej+1 kA = kA1/2 ej+1 k = A1/2 RA−1/2 A1/2 ej k2 ≤ k |A1/2 RA
{z
−1/2 1/2 1/2
} k2 kA ej k2 = |ρ(A {z RA−1/2 ) kej kA ,
}
∈§n (R)
=ρ(R)= K−1
K+1

où R := (I − αR A). 

5.4 Espaces de Krylov

Définition 3.5.2 (Espace de Krylov) Soient r ∈ Rn , j ∈ N et A ∈ Mn (R). On appelle espace de Krylov de A
associé à r et d’ordre j l’ensemble
KJ := Vect(r, Ar, . . . , Aj−1 r),
noté parfois Kj (A, r). On pose par convention K0 := {0}.
Proposition 3.5.3 Soient r ∈ Rn , A ∈ Mn (R) et la suite (Kj )j vérifiant
K0 ⊂ K1 ⊂ · · · ⊂ Kl = Kl+1 = . . . ,
pour un certain l ≤ n. On a donc 
j pour j ≤ l
dim Kj =
l pour j > l.

Exercice 3.5.4 Considérer l’ensemble {j ∈ N | (r, Ar, . . . , Aj−1 r} est une famille libre}, et l son maximum.
Remarque 3.5.5
• ∀j ∈ N, Kj (A, r) = {P (A)r | P ∈ Rj−1 [X]}.
• Les méthodes GPF et GPO vérifient la propriété : en posant pour j ∈ N, rj := b − Axj , on a

rj ∈ Kj+1 (A, r0 )
xj ∈ x0 + Kj (A, r0 )
Ceci se démontre par récurrence, le cas j = 0 étant vérifié. Pour l’hérédité, on utilise le fait que pour j ∈ N,
xj+1 = xj +αj rj ,
|{z} |{z}
∈x0 +Kj ∈Kj+1

donc
rj+1 = rj −αj Arj ∈ Kj+2 (A, r0 ).
|{z} |{z}
∈Kj+1 ∈Kj+2

Pourrait-on modifier les algorithmes pécédents de sorte que xj = x = A−1 b pour j assez grand (typiquement, n),
avec toujours xj ∈ x0 + Kj (A, r0 ), i. e. construire P ∈ R[X] tel que x − x0 ∈ P (A)r0 ?
Proposition 3.5.6 Soient A ∈ GLn (R), x0 ∈ Rn et Kl (A, r0 ) l’espace de Krylov maximal associé à r0 := b − Ax0 .
Alors x = A−1 b ∈ x0 + Kl (A, r0 ), un sous-espace vectoriel affine de dimension l ≤ n.
Preuve. La famille liée (r0 , Ar0 , . . . , Al r0 ) par définition de l’entier l. Alors il existe α0 , . . . , αl−1 ∈ R tels que
Al r0 = α0 r0 + · · · + αl−1 Al−1 r0 .
Pour montrer que α0 6= 0, on procède par l’absurde. Si α0 = 0, alors
A(Al−1 r0 ) = A(α1 r0 + · · · + αl−1 al−2 r0 ),
ce qui est impossible car (r0 , . . . , Al−1 r0 ) est libre. Ainsi,
1
x= (Al−1 r0 − α1 r0 − · · · − Al−2 r0 ) + x0 ∈ x0 + Kl (A, r0 ).
α0


25
ANU CHAPITRE 3. MÉTHODES ITÉRATIVES

5.5 Algorithme du gradient conjugué

Soit A ∈ Sn++ (R). On veut construire récursivement les projetés A-orthogonaux de x = A−1 b sur les sous-espaces
affines croissants x0 + Kj (A, r0 ), pour j ∈ N à partir d’un choix x0 ∈ Rn , et ce sans connaître x. Par la proposition
précédente, le nème itéré (au plus tard) coïncide avec x. l’algorithme converge exactement en l étapes. Mais comment
construire un telle suite (xj )j , de sorte que
xj = argmin f.
x0 +Kj (A,r0 )

Pour obtenir xj+1 à partir de xj , on va utiliser les directions A-orthogonales des espaces de Krylov.
Soient x0 ∈ Rn , r0 := b − Ax0 , p0 := r0 , la direction de descente. La famille (p0 ) est une base A-orthogonale de
K1 (A, r0 ). L’objectif est de construire (p0 , . . . , pl−1 ) une base A-orthogonale de Kl (A, r0 ) en drapeau. Connaissant xj ,
on cherche xj+1 = xj + αj pj , pour un αj bien choisi. On sait que xj = pEj (x), avec Ej := x + Kj (A, r0 ). Alors

∀i ≤ j − 1, hxj+1 − x, pi iA = hxj+1 − xj , pi iA + hxj − x, pi iA = 0.

| {z } | {z }
=0 car hpj,πiA =0 =0 car xj =pEj (x)

Pour i = j,
hxj+1 − x, pj iA = αj kpj k2A + hxj − x, pj iA = αj kpj k2A + hAxj − b, pj i2 .
Ainsi, on pose
hrj , pj i
αj := ,
kpj k2A
avec rj := b − Axj .
Comment trouver (pj )j ? On utilise l’agorithme de Gram-Schmidt pour déduire pj de (p0 , . . . , pj−1 ) base de Kj et
de rj qui complète cette base en une base de Kj+1 . Cette orthogonalisation est « courte » :
j−1
X hrj , pi i
pj = rj − βi p i , avec βi = .
i=0
kpi k2A

Pour i ≤ j − 2, hrj , pi iA = hrj , Api i2 , avec Api ∈ Kj−1 . or rj = b − Axj ∈ Kj⊥A , donc hrj , piiA = 0. Donc
pj = rj − βj−1 pj−1 , avec
hrj , pj−1 iA krj k2
βj−1 = = − .
kpj−1 k1A krj−1 k2

26
Chapitre 4

Approximation spectrale

1 Motivation
Soit A ∈ Mn (K), K = R ou C. Comment trouver les éléments propres de A ? Ce problème pose plusieurs difficultés.
• Pour n grand, on n’a pas d’espoir de résoudre algébriquement.
• L’équation est non linéaire.
• Si on ne connaît A qu’approximativement, le résultat est-il suffisant ?

2 Méthodes numériques
2.1 Méthode de la puissance
Soit A ∈ Mn (R), pour trouver la valeur propre λ ∈ C tel que |λ| = ρ(A), on s’appuie sur l’observation vague Ak x '
λ x pour k grand (c’est faux). Plus précisément, Ak est obtenue à partir des sous-espaces propres ou caractéristiques
k

de A et des λki , λi ∈ σ(A).

On considère l’algorithme suivant. Soit x0 ∈ Cn , on pose pour tout k ∈ N,
Axk
xk+1 := , vk := hxk , Axk i.
kAxk k2

Théorème 4.2.1 Soit A ∈ Mn (C), de valeurs propres λ1 , . . . , λd , avec

|λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λd |.

On dit que λ1 est dominante. Soit x0 ∈ / ⊕di=2 Ker((A − λi In )mi ). On peut interpréter ceci en disant que x0 possède une
composante normale dans Ker((A − λ1 In )m1 ) dans Cn = ⊕di=1 Ker((A − λi In )mi ).
Alors les suites (xk )k et (vk )k sont bien définies et on a

lim vk = λ1 ,
k→+∞

et
lim qk = e ∈ Ker(A − λ1 In ),
k→+∞

où  ¯ k
λ1
qk := xk .
|λ1 |
De plus, (xk )k « converge »vers Ker(A − λ1 In ).

27