Méthodes d’Analyse Numérique pour l’Année

Préparatoire II

–Cours de Mathématiques Appliquées–

Analyse Numérique, Analyse Numérique matricielle

Mohamed ADDAM
Professeur de Mathématiques

École Nationale des Sciences Appliquées d’Al Hoceima

–ENSAH–
Année Universitaire 2021/2022
[Link]@[Link]
[Link]@[Link]

c Mohamed ADDAM.

21 février 2022
2
Table des matières

1 Analyse numérique matricielle 5

1.1 Spectre et rayon spectral d’une matrice, Matrice positive . . . . . . . . . . . . . . . . . . 5
1.1.1 Matrice positive et matrice définie positive . . . . . . . . . . . . . . . . . . . . . 6
1.1.2 Valeurs singulières d’une matrice . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.3 Décomposition en valeurs singulières . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.4 Pseudo-inverse de Moore-Penrose . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Normes matricielles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.1 Représentation des entiers et des réels sur ordinateur . . . . . . . . . . . . . . . . 10
1.3.2 Effet de la représentation des réels et les erreurs d’arrondi sur la résolution de Ax = b 13
1.3.3 Propriétés du conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

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

2.1 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 résolution du système triangulaire . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.2 Principe des méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.3 Élimination de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.4 Factorisation LU d’une matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.5 Factorisation de Cholesky où bien factoristion BB T . . . . . . . . . . . . . . . . 23
2.2 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.2 Comparaison des méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.3 Principales méthodes itératives classiques . . . . . . . . . . . . . . . . . . . . . . 28
2.2.4 Etude de la convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3 Interpolation et approximation polynômiale 35

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Interpolation polynomiale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.1 Interpolation polynomiale de Lagrange . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Détermination du polynôme d’interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3.1 Cas où n = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3.2 Cas général . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4 Interpolation par les différences divisées . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.5 Erreur d’interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3
4 TABLE DES MATIÈRES

4 Intégration et dérivation numérique 45

4.1 Introduction et outils de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Formule de quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3 Quadratures interpolatoires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4 La méthode des rectangles à gauche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4.1 Formule du rectangle ou du point milieu . . . . . . . . . . . . . . . . . . . . . . . 47
4.4.2 Formule du trapèze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.4.3 Formule de Cavalieri-Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5 Méthode des moindres carrés et optimisation quadratique 53

5.1 Maxima et minima de fonctions de deux variables . . . . . . . . . . . . . . . . . . . . . . 53
5.1.1 Gradient d’une application et Matrice hessienne d’une F.P.V . . . . . . . . . . . . 53
5.1.2 Approximations linéaire et quadratique : Formule de Taylor . . . . . . . . . . . . 54
5.1.3 Points critiques d’une application . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.1.4 Maxima et minima des fonctions de n variables . . . . . . . . . . . . . . . . . . . 56
5.2 Fonctions quadratiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.1 Forme linéaires et bilinéaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.2 Équivalence entre la résolution d’un système linéaire et la minimisation quadratique 58
5.3 Application aux moindres carrés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.1 Approximation par la droite des moindre carrés . . . . . . . . . . . . . . . . . . . 60
5.3.2 Interprétation géométrique : projection sur un sous-espace . . . . . . . . . . . . . 60
Chapitre 1

Analyse numérique matricielle

1.1 Spectre et rayon spectral d’une matrice, Matrice positive

Soit A = (ai,j )1≤i,j≤n une matrice carrée de taille n × n.
X n
1. La trace de A est tr(A) = ai,i .
i=1
2. Les valeurs propres de A sont les n racines réelles ou complexes (λi )1≤i≤n du polynôme caractéris-
tique P de A. Le spectre de A, noté Sp(A) est l’ensemble de tous les valeurs propres de A :
Sp(A) = {λi : 1 ≤ i ≤ n}

3. La matrice A est diagonale si ai,j = 0 pour i 6= j, on la note

A = diag(aii ) = diag(a11 , a22 , . . . , ann ).
On rappelle les propriétés suivantes :
Xn Yn
1. tr(A) = λi , dét(A) = λi .
i=1 i=1
2. tr(AB) = tr(BA), tr(A + B) = tr(A) + tr(B).
3. det(AB) = det(BA) = det(A)det(B).

Définition 1.1.1 On appelle le rayon spectral de la matrice A, noté %(A), le nombre réel positif
%(A) = max{|λi | : 1 ≤ i ≤ n}

Définition 1.1.2 Une matrice A est

1. Symétrique si A est réelle et A = AT ;
2. hermitienne si A = A∗ ;
3. Orthogonale si A est réelle et AAT = AT A = I ;
4. Unitaire si AA∗ = A∗ A = I ;
5. Normale si AA∗ = A∗ A.
une matrice A est dite singulière si elle n’est pas inversible.

Propriété 1.1.1 Si A et B sont deux matrices inversibles, alors (AB)−1 = B −1 A−1 , (AT )−1 = (A−1 )T ,
(A∗ )−1 = (A−1 )∗ .
5
6 CHAPITRE 1. ANALYSE NUMÉRIQUE MATRICIELLE

1.1.1 Matrice positive et matrice définie positive

Définition 1.1.3 Soit A une matrice
1. La matrice A est définie positive si
(Ax, x) > 0, ∀x ∈ E − {0E }
et (Ax, x) = 0, ⇔ x = 0E .
2. La matrice A est positive où semi-définie positive si
(Ax, x) ≥ 0, ∀x ∈ E − {0}
Théorème 1.1.1 Une matrice hermitienne A est définie positive (resp. positive), si et seulement si toutes
ses valeurs propres sont > 0 (resp. ≥ 0).
Démonstration. soit A une matrice hermitienne et x 6= 0 un vecteur dans E.
(Ax, x) = λ(x, x) = λkxkE


1.1.2 Valeurs singulières d’une matrice

Définition 1.1.4 Soit A une matrice. On appelle valeurs singulières d’une matrice A, les racines carrées
positives des valeurs propres de la matrice hermitienne A∗ A (ou AT A si la matrice A est réelle).
Remarque 1.1.1 Les valeurs propres de la matrice hermitienne A∗ A sont toujours ≥ 0 puisque
A∗ Ax = λx, x 6= 0 ⇒ (A∗ Ax, x) = λkxkE ,
les valeurs singulières sont toutes > 0 si et seulement si la matrice A est inversible

1.1.3 Décomposition en valeurs singulières

Toute matrice peut être réduite sous forme diagonale en la multipliant à droite et à gauche par des
matrices unitaires bien choisies. Plus précisément on a le résultat suivant :
Propriété 1.1.2 Soit A ∈ Cm×n . Il existe deux matrices unitaires U ∈ Cm×m et V ∈∈ Cn×n telles que
U ∗ AV = diag(σ1 , σ2 , . . . , σp ) ∈ Cm×n , p = min(m, n) (0.1)
et σ1 ≥ σ2 ≥ . . . ≥ σp ≥ 0.
La relation (0.1) est appelée décomposition en valeurs singulières (SVD) de A et les scalaire (σi ) sont
appelés valeurs singulières de A.
♣ Si A est une matrice réelle, alors U et V sont aussi des matrices réelles et on peut remplacer U ∗ par U T .
Les valeurs singulières sont caractérisées par
p
σi = λi , où λi ∈ Sp(A∗ A), i = 1, . . . , p. (0.2)
– Si A ∈ Cn×n est une matrice hermitienne de valeurs propres λ1 , λ2 , . . . , λn , alors d’après (0.2) les
A coïncident avec les modules des valeurs propres de A. En effet, puisque
valeurs singulières de p
∗ 2
AA = A , on a σi = λ2i = |λi | pour i = 1, . . . , n.
– Si
σ1 ≥ σ2 ≥ . . . ≥ σr ≥ σr+1 = . . . = σp = 0,
alors le rang de A est r (rang(A) = r), le noyau de A est le sous-espace vectoriel engendré par
les vecteurs colonnes de V (Ker(A) = {vr+1 , . . . , vn }), et l’image de A est le sous-espace vectoriel
engendré par les vecteurs colonnes de U (Im(A) = {u1, . . . , ur }).
1.2. NORMES MATRICIELLES 7

1.1.4 Pseudo-inverse de Moore-Penrose

Définition 1.1.5 Supposons A ∈ Cm×n soit de rang r et qu’elle admette une décomposition en valeurs
singulières du type U ∗ AV = Σ. La matrice A† = V Σ† U ∗ est appelée matrice pseudo-inverse de Moore-
Penrose, où
 
1 1 1
†
Σ = diag , , . . . , , 0, . . . , 0 ∈ Cm×n , p = min(m, n) (0.3)
σ1 σ2 σr

la matrice A† est aussi appelée matrice inverse généralisée de A, on a


† (AT A)−1 AT , si rang(A) = n < m,
A =
A−1 , si rang(A) = n = m.

Exemple 1.1.1 On peut se demander si l’inverse de Moore-Penrose peut être utilisée dans des cas pratique.
On cherche à estimer une valeur en x0 à partir de deux valeurs situées au même point x1 . On considère le
système simple     
a a w1 b
=
a a w2 b
où a et b sont des données réelles.
La matrice A étant singulière, on va recourir à l’inverse généralisée de Moore-Penrose. Soit
 2 
T T 2a 2a2
AA = A A = =C
2a2 2a2

On a
dét(C) = 0 ⇒ λ2 = 0;
tr(C) = 4a2 ⇒ λ1 = 4a2 ,
et une matrice de vecteurs propres normés
!
√1 − √12
Q= 2
√1 √1
2 2

La solution du système au sens de l’inverse généralisée est

b
† † T 2a
w = A b = QΣ Q b = b
2a

1.2 Normes matricielles

Soit K = R ou C.

Définition 1.2.1 Une norme matricielle est une application k . k : K(m×n) −→ [0, +∞[ telle que :
i) kAk ≥ 0, ∀A ∈ K(m×n) et kAk = 0 si et seulement si A = 0 ;
ii) kαAk = |α|kAk, ∀A ∈ K(m×n) , ∀α ∈ K (Propriété d’homogénéité) ;
iii) kA + Bk ≤ kAk + kBk, ∀A, B ∈ K(m×n) (Inégalité triangulaire).
8 CHAPITRE 1. ANALYSE NUMÉRIQUE MATRICIELLE

Définition 1.2.2 On dit que la norme matricielle k . k est compatible ou consistante avec la norme vecto-
rielle k · k si
kAxk ≤ kAkkxk, ∀A ∈ K(m×n) ∀x ∈ Kn .
Plus généralement, on dit que trois normes, toutes notées k . k et respectivement définies sur Km , Kn , et
K(m×n) , sont consistantes si ∀x ∈ Kn , ∀y ∈ Km et A ∈ K(m×n) tels que Ax = y, on a kyk ≤ kAkkxk. 
Pour qu’une norme matricielle soit intéressante dans la pratique, on demande généralement qu’elle
possède la propriété suivante :

Définition 1.2.3 On dit qu’une norme matricielle k . k est sous-multiplicative si ∀A ∈ K(n×m) , ∀B ∈

K (m×q)

kABk ≤ kAkkBk. (0.4)



Exemple 1.2.1 Soit k . kN l’application définie par

kAkN = max |aij |.
i,j

On considère deux matrices A et B données par

 
1 1
A= = B.
1 1
On peut vérifier facilement que k . kN est une norme et qu’elle ne satisfait pas l’inégalité (0.4) puisque
2 = kABkN > kAkN kBkN = 1.
D’où la norme k . kN n’est pas une norme sous-multiplicative.

Norme de Frobenius
La norme !1/2
m X
X n
p
kAkF = |aij |2 = tr(A∗ A)
j=1 i=1

est une norme matricielle appelée norme de Frobenius (ou norme euclidienne dans C(n×n) et elle est
compatible avec la norme vectorielle euclidienne k . k2. En effet,
m X n 2 m n n
!
X X X 2
X 2
kAxk22 = aij xj ≤ |aij | |xj | = kAk2F kxk22 .
i=1 j=1 i=1 j=1 j=1
√
On peut remarquer que pour cette norme, on a kIn kF = n.

Théorème 1.2.1 Soit k.k une norme vectorielle. La fonction

kAxk
kAk = sup (0.5)
x6=0 kxk
est une norme matricielle. On l’appelle norme matricielle subordonnée ou associée à la norme vectorielle
k.k. On l’appelle aussi parfois norme matricielle naturelle, ou encore norme matricielle induite par la
norme vectorielle k.k.
1.2. NORMES MATRICIELLES 9

Démonstration. Commençons par remarquer que (0.5) est équivalente à

kAk = sup kAxk. (0.6)

kxk=1

x
Pour tout x 6= 0, on peut définir un vecteur unitaire u = kxk
, de sorte que (0.5) s’écrive

kAk = sup kAuk = kAwk, kwk = 1.

kuk=1

cela étant, vérifions que (0.5) est effectivement une norme, en utilisant les conditions d’une norme matri-
cielle de la définition (1.2.1). 

Exemples de normes remarquables

D’autres exemples de normes remarquables, il s’agit de normes matricielles subordonnées fournies par
les p-normes :

kAxkp
kAkp = sup (0.7)
x6=0 kxkp

La 1-norme et la norme infinie se calculent facilement :

m
X
kAk1 = max |aij |, (0.8)
1≤j≤n
i=1

n
X
kAk∞ = max |aij |, (0.9)
1≤i≤m
j=1

On a les propriétés suivantes :

1. kAk1 = kAT k∞
2. Si A est matrice symétrique réelle, alors kAk1 = kAk∞ . 
La 2-norme ou norme spectrale mérite une discussion particulière vu son intérêt pour des applications
pratiques.

Théorème 1.2.2 Soit σ1 la plus grande valeur singulière de A. Alors

p p
kAk2 = %(A∗ A) = %(AA∗ ) = σ1 .

En particulier, si A est hermitienne (ou symmétrique réelle), alors

kAk2 = %(A),

tandis que si A est unitaire alors kAk2 = %(A) = 1.

Démonstration. Puisque A∗ A est hermitienne, alors il existe une matrice unitaire U telle que

U ∗ A∗ AU = diag(µ1 , µ2, . . . , µn )
10 CHAPITRE 1. ANALYSE NUMÉRIQUE MATRICIELLE

où µi sont les valeurs propres positive de A∗ A. Soit y = U ∗ x, alors

s s
(A∗ Ax, x) (U ∗ A∗ AUy, y)
kAk2 = sup = sup
x6=0 (x, x) x6=0 (y, y)
v
u n n
uX X q
= sup t 2
µi |yi | / |yi |2 = max |µi | = σ1 .
x6=0 1≤i≤n
i=1 i=1

Si A est hermitienne, les mêmes considérations s’appliquent directement à A. En fin si A est unitaire

kAxk22 = (Ax, Ax) = (x, A∗ Ax) = kxk22

et donc kAk2 = 1. Enfin, si A est unitaire 

Exercice 1.2.1 1. Soit B une matrice carrée. Montrer que les conditions suivantes sont équivalentes :
(a) lim B k = 0 ;
k→+∞

(b) lim B k x = 0 pour tout vecteur x ;

k→+∞
(c) %(B) < 1 ;
(d) kBk < 1 pour au moins une norme matricielle subordonnée k · k.
2. Soit B une matrice carrée, et k · k une norme matricielle quelconque. Montrer que

lim kB k k1/k = %(B).

k→+∞

1.3 Conditionnement
Les sources des erreurs numériques sont :
1. Erreur d’arrondi dues à la représentation des réels sur la machine,
2. Erreur sur les données (données qui proviennent d’autres calculs)
3. Erreur de troncature (faites dans les méthodes itératives : on remplace la valeur exacte par une valeur
approchée)

1.3.1 Représentation des entiers et des réels sur ordinateur

Représentation des entiers

Soit a ∈ N,
n−1
X
a= ai 2i , avec ai ∈ {0, 1}
i=0
la représentation binaire (représentation en base 2) de a est

a = an−1 an−2 . . . a2 a1 a0 .

Exemple 1.3.1 1. a = 13 = 23 + 22 + 20 = 1 × 23 + 1 × 22 + 0 × 21 + 1 × 20
a est représenté sur la machine par (1101)

13 = (1101)
1.3. CONDITIONNEMENT 11

2. 226 = 27 + 26 + 25 + 0 × 24 + 0 × 23 + 0 × 22 + 21 + 0 × 20
donc 226 = (11100010)

Le chiffre 0 où 1 représente un bit.

bit=binary disit=chiffre binaire.

Exemple 1.3.2 Sur l’ordinateur, les entiers sont représentés par un ensemble de bits.

1. a = 13 est représenté par 4 bits.

2. a = 226 est représenté par 8 bits

Exemple 1.3.3 Avec une représentation sur 16 bits.

a = a15 a14 . . . a2 a1 a0 , ai ∈ {0; 1}


1 entier négatif,
le binaire a15 désigne le signe
0 entier positif

Exemple 1.3.4 a = 226 est représente sur 16 bits par

a = (0000000011100010)

a = −226 est représenté sur 16 bits par

a = (1000000011100010)

le plus grand entier qu’on peut représenter sur la machine avec 16 bits est :

g = 214 = 213 + . . . + 22 + 21 + 1 = 215 − 1 = 32768 − 1 = 32767

Remarque 1.3.1 La somme de deux entiers sur la machine peut donner un nombre négatif

Exemple 1.3.5 Mathématiquement 32767 + 1 = 32768.

Sur la machine :

32767 = 0111111111111111
+
1 = 0000000000000001
32768 = 1000000000000000

ce nombre est négatif car le bit de a15 est 1.

12 CHAPITRE 1. ANALYSE NUMÉRIQUE MATRICIELLE

Représentation des réels

Représentation en vergule flottantes : une opération élémentaire sur des réels est appelée flop où floaling
point.
soit x ∈ R, x est représenté par
x = m × bp
où m est la montisse, b est la base et p est l’exposant, avec la condition
1
≤ |m| < 1
b
donc m = 0.d1 d2 . . . dt avec 0 ≤ di < b et d1 6= 0.
t désigne la précision (le nombre des chiffres après la vergule)
t
X
m= di b−i .
i=1

Exemple 1.3.6 En base de 10 c’est-à-dire b = 10.

Le réel x = 455.321 est représenté sur la machine avec x
e où

x
e = 0.455321 |{z}
e3
103

ici ”e” signifie l’exposant, et t = 6 est la précision de 6 chiffres.

Exercice 1.3.1 Supposons que ` ≤ p ≤ u. Calculer le plus grand et le plus petit réel.

Erreurs de représentation

La représentation des réels sur la machine n’est pas exacte.

Exemple 1.3.7 supposons que b = 10 et t = 3

1
1. y = 3
= 0.333333 . . . 3 est représenté par ye = 0.333
3254
2. y = 100
= 32.54 est représenté par ye = 0.325e2
on ne représente que les valeurs approchées.

Erreurs d’arrondies sur les opérations élémentaires (+, −, ∗, /)

les opérations flottants ne sont pas associatives

(e
x + ye) + ze 6= x
e + (e
y + ze).

Exemple 1.3.8 b = 10, t = 3, x = 10−3 , y = 1 et z = −1.

Mathématiquement, on a x + y + z = 10−3 .
Sur la machine, on a x]
+ y = 1.001 → 0.1001e1 = 0.1e1 car la précision t = 3.
1. (x]
+ y) + ze = 0
2. x y + ze) = 10−3 = 0.1e − 2.
e + (e
 1.3. CONDITIONNEMENT 13

1.3.2 Effet de la représentation des réels et les erreurs d’arrondi sur la résolution
de Ax = b

3x − 7.0001y = 0.9998,
Exemple 1.3.9 soit le système suivant : admet la solution unique
3x − 7y = 1

x = 31 ,
y = 1−0.9998
0.0001
   
3 −7.0001 0.9998
A= , b=
3 −7 1
supposons qu’on travaille sur une machine avec t = 4, alors 7.0001 sera représentée par 0.70001e1 →
0.7e1.
sur la machine, on a à résoudre le système suivant qui est singulier

0.3e1x − 0.7e1y = 0.9998,
0.3e1x − 0.7e1y = 1

ce système n’a pas de solution.

   
e= 3 −7 eb = 0.9998
A ,
3 −7 1

ceci montre qu’une perturbation sur la matrice A induit une matrice A e où le système n’a pas de solution.
Perturbation sur b : Une perturbation sur b peut conduire à des résultats qui ne sont pas justes :
   
3 −7.0001 1
A= , eb =
3 −7 1
 
3x − 7.0001y = 1, x = 31 ,
le système suivant admet la solution unique mais elle n’est pas une
3x − 7y = 1 y = 0
solution juste.
Les mauvais résultats obtenus sont dûs au fait que la matrice A est mal conditionnée.

Définition 1.3.1 une matrice est mal conditionnée si une petite perturbation (modification des données)
conduit à des résultats différents.

1er cas : Soit à résoudre le système Ax = b. Soit 4b la perturbation sur b, on résout donc le système

A(x + 4x) = b + 4b,
4x : perturbation sur x

x∗ = x + 4x est la solution obtenue après modification de b.

A(x + 4x) = b + 4b
Ax + A4x = b + 4b ⇒ A4x = 4b
⇒ 4x = A−1 4b

soit k · k une norme matricielle induite :

k4xk ≤ kA−1 kk4bk

 14 CHAPITRE 1. ANALYSE NUMÉRIQUE MATRICIELLE

d’autre part, on a
1 kAk
kbk = kAxk ≤ kAkkxk ⇒ ≤
kxk kbk
on déduit que
k4xk k4bk
≤ kAkkA−1 k
kxk | {z } kbk
| {z } κ(A)=cond(A) | {z }
erreur relative sur x erreur relative sur b

La quantité κ(A) = cond(A) = kAkkA−1 k s’appelle le conditionnement de A ≥ 1. Si κ(A) est plus

proche de 1 alors, une petite perturbation sur b ( k4bk
kbk
très petit) entraine des petites perturba-
tions sur x ( k4xk
kxk
très petit)
2eme cas : Supposons qu’on a une perturbation sur A

(A + 4A)(x + 4x) = b,
4A : perturbation sur A

Ax + A4x + 4A(x + 4x) = b

Ax + A4x + 4A(x + 4x) = b ⇒ 4A(x + 4x) = −A4x
⇒ 4x = −A−1 4A(x + 4x)
soit k · k une norme matricielle induite :
k4xk ≤ kA−1 kk4Akkx + 4xk
on déduit que
k4xk k4Ak
≤ kAkkA−1 k
kx + 4xk | {z } kAk
| {z } κ(A)=cond(A) | {z }
erreur relative sur x erreur relative sur A

Si κ(A) est plus proche de 1, alors une petite perturbation sur A entraine des petites perturbations sur
le résultat x.

Définition 1.3.2 Soit A une matrice.

1. La matrice A est bien conditionnée si κ(A) est plus proche de 1.
2. La matrice A est mal conditionnée si κ(A) est très grand.

1.3.3 Propriétés du conditionnement

Soit A une matrice matrice carrée inversible.
1. κ(A) ≥ 1, κ(A) = κ(A−1 ) et κ(αA) = κ(A) pour tout α ∈ K.
2. κ2 (A) = kAk2 kA−1 k2 = σσn1 où σ1 est la plus grande valeur singulière de A et σn est la plus petite
valeur singulière de A.
3. Le matrices orthogonales sont bien conditionnées, κ2 (A) = 1.

Exemple 1.3.10    
3 −7.0001 −1 1 −7 7.0001
A= , A =
3 −7 3.10−4 −3 3
10.0001×14.0001
On a κ∞ (A) = kAk∞ kA−1 k∞ = 3
× 104 qui est très grand.
Chapitre 2

Résolution de systèmes : Méthodes directes et

méthodes itératives

Soit A une matrice carrée d’ordre n, inversibles à coefficients dans R et soit b un vecteur à n compo-
santes. l’objectif de ce chapitre est de résoudre le système linéaire Ax = b à n équations.

2.1 Méthodes directes

On appelle méthode directe de résolution du système linéaire Ax = b, toute méthode permettant d’ob-
tenir la solution en un nombre fini d’opérations arithmitiques élémentaires sur des nombres réels (additions,
soustractions, multiplications, divisions ) et éventuellement l’extraction des racines carrées.

2.1.1 résolution du système triangulaire

Supposons que A est une matrice triangulaire suppérieure

 
a1,1 a1,2 . . . . . . a1,n
 0 a2,2 a2,3 a2,n 
 .. .. 
 .. .. .. 
A= . . . . . 
 .. .. .. 
 . 0 . . an−1,n 
0 ... ... 0 an,n

La résolution du système Ax = b s’effectue par (back substitution où backward substitution).

Elle consiste à calculer xi en fonction de xn , xn−1 , . . . , xi+1 .
Yn
A étant une matrice inversible, alors det(A) = ai,i 6= 0 donc ai,i 6= 0 pour tout 1 ≤ i ≤ n.
i=1
bn
– Calculons d’abord xn =
an,n
– On reporte la valeur dans l’équation précédente pour calculer

bn−1 − an−1,n xn bn−1 an,n − an−1,n bn

xn−1 = =
an−1,n−1 an,n an−1,n−1
15
16 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

– Ainsi de suite, plus généralement, on obtient

n
X
bi − ai,j xj
j=i+1
xi = , pour i = n, n − 1, . . . , 1
ai,i
le coût de la résolution est mesuré en nombre d’opérations élémentaires appelé aussi coût de calcul. le
calcul de xi (1 ≤ i ≤ n) nécessite 2(n − i) + 1 opérations. On déduit que le coût de calcul des xi est
n
X n
X
2
(2(n − i) + 1) = 2n − 2 i + 2n = 2n2 − n(n + 1) + n = n2
i=1 i=1

opérations.

2.1.2 Principe des méthodes directes

Les résolutions directes consistent à transformer le système Ax = b en le système Rx = c où R est une

matrice triangulaire supérieure qu’on sait résoudre facilement.
On multiplie A par des matrices bien choisies, soit M le produit des ces matrices, on transforme alors le
système Ax = b en un nouveau système

MAx = Mb = c,

On détermine M de telle sorte que MA soit triangulaire où diagonale.

2.1.3 Élimination de Gauss

Pour résoudre le système Ax = b, le principe de la méthode de Gauss consiste à :

1. Phase délimination : prémultiplier A et le second membre b par des matrices bien choisies pour
transformer le système Ax = b en un système triangulaire supérieur (Rx = c) donc facile à résoudre
(ohase de triangularisation)
2. Remontée par back substitution : Résoudre par la méthode de la remontée back substitution du
système Rx = c obtenu.
Description de la méthode :


 a1,1 x1 + a1,2 x2 + . . . + a1,n xn = b1 , (1)

 a2,1 x1 + a2,2 x2 + . . . + a2,n xn = b2 , (2)
Ax = b ⇔ (S (1) ) : .. .


 . = .. (i)
 a x + a x + ...+ a x = bn , (n)
n,1 1 n,2 2 n,n n

avec A = (ai,j ) 1≤i≤n , x = (x1 , x2 , . . . , xn )T et b = (b1 , b2 , . . . , bn )T .

1≤j≤n

– 1er -étape : élimination de x1 des équations (2) jusqu’à (n). On suppose que a1,1 6= 0 (sinon on
permute l’équation (1) avec une équation (i) du système tel que ai,1 6= 0).
Puisque a est inversible, alors il existe i tel que ai,1 6= 0 et pour éliminer x1 de l’équation (i) (2 ≤ i ≤
n), on effectue la combinaison linéaire suivante entre l’équation (1) et l’équation (i) : on remplace
l’équation (i) par l’équation
2.1. MÉTHODES DIRECTES 17

équation(1)
équation(i) − . ai,1 ,
a11
Donc l’équation (i) devient sous la forme suivante
     
a1,1 a1,2 a1,j
ai,1 − ai,1 x1 + ai,2 − ai,1 x2 + . . . + ai,j − ai,1 xj
a1,1 a1,1 a1,1
 
a1,n b1
+ ai,n − ai,1 xn = bi − ai,1
a1,1 a1,1
Posons
ai,1 (2) (2)
`i,1 = , ai,j = ai,j − `i,1 a1,j et bi = bi − `i,1 b1
a1,1
alors le système devient :


 a1,1 x1 + a1,2 x2 + . . . + a1,n xn = b1 ,



 (2) (2) (2)
(2)
(S ) : 0 x1 + a2,2 x2 + . . . + a2,n xn = b2 ,

 .. .


 . = ..
 (2) (2) (2)
0 x1 + an,2 x2 + . . . + an,n xn = bn ,

Formulation matricielle : Posons

 
1 0 ... ... 0
 −`2,1 1 0 ... 0 
 . .. ai,1
 .. .. 
L1 =  .. 0 . . . avec `i,1 = , 2≤i≤n
 . .. .. ..  a1,1
 .. . . . 0 
−`n,1 0 ... 0 1

A(1) = A, b(1) = b
A(2) = L1 A,
le système (S (2) ) équivalent à

A(2) x = L1 A(1) x = L1 b(1) = b(2)

D’où on obtient le système

Ax = b ⇔ A(2) x = b(2)
où    
a1,1 a1,2 . . . ... a1,n b1
 (2) (2)
0 a2,2 a2,3
(2)
a2,n   (2) 
   b2 
 .. .. .. .. ..   .. 
A(2) =
 . . . . . 
 et b(2) = . 
 .. .. ..   (2) 
 . . . a(2)   bn−1 
n−1,n
(2) (2) (2)
0 an,2 . . . . . . an,n bn
(2)
– 2eme -étape : élimination de x2 des équations (3) jusqu’à (n). On suppose que a2,2 6= 0 (sinon on
(2)
permute l’équation (2) avec une équation (i) du système tel que ai,2 6= 0).
Pour éliminer x2 de l’équation (i) (3 ≤ i ≤ n) du système A(2) x = b(2) , on effectue la combinaison
linéaire suivante entre l’équation (2) et l’équation (i) : on remplace l’équation (i) par l’équation
18 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

équation(2) (2)
équation(i) − (2)
. ai,2
a2,2
Posons
(2)
ai,2 (3) (2) (2) (3) (2) (2)
`i,2 = (2)
, ai,j = ai,j − `i,2 a2,j et bi = bi − `i,2 b2
a2,2
alors le système devient :

 a1,1 x1 + a1,2 x2 + . . . + . . . + a1,n xn = b1 ,





 (2) (2) (2)

 0 x1 + a2,2 x2 + . . . + . . . + a2,n xn = b2 ,

(S (3) ) :

 (3) (3)
0 x1 + 0 x2 + a3,3 x3 + . . . + a3,n xn = b3 ,
(3)



 .. .


 . = ..
 (3) (3) (3)
0 x1 + 0 x2 + an,3 x3 + . . . + an,n xn = bn ,

Formulation matricielle : Posons

 
1 0 ... ... ... 0
 0 1 0 ... ... 0 
 
 0 −`3,2 1 0 ... 0  (2)
ai,2
 ..
L2 =  .. ..  avec `i,2 = , 3≤i≤n
 0 −`4,2 0 . . . (2)
a2,2
 .. .. .. .. 
 0 . . . . 0 
0 −`n,2 0 ... 0 1

A(1) = A, b(1) = b
A(2) = L1 A, et A(3) = L2 A(2) = L2 L1 A,
le système (S (3) ) équivalent à

A(3) x = L2 L1 A(1) x = L2 L1 b(1) = b(3)

D’où on obtient le système

Ax = b ⇔ A(3) x = b(3)
où  
  b1
a1,1 a1,2 a1,3 . . . a1,n  (2) 
 0 a(2) b2
a2,3 . . . a2,n   
(2) (2)
 2,2   (3) 
 . (3)   b3 
A(3) = . 0
(3)
a3,3 . . . a3,n  et b(3) = 
 .   ..
. 
 .. .. .. .. ..   
 . . . . .   (3) 
(3) (3)  bn−1 
0 0 an,3 . . . an,n (3)
bn
(k)
– k eme -étape : élimination de xk des équations (k + 1) jusqu’à (n). On suppose que ak,k 6= 0 (sinon
(k)
on permute l’équation (k) avec une équation (i) (i > k) du système tel que ai,k 6= 0).
Pour éliminer xk de l’équation (i) (k ≤ i ≤ n) du système A(k) x = b(k) , on effectue la combinaison
linéaire suivante entre l’équation (k) et l’équation (i) : on remplace l’équation (i) par l’équation
2.1. MÉTHODES DIRECTES 19

équation(k) (k)
équation(i) − (k)
. ai,k
ak,k
Posons
(k)
ai,k (k+1) (k) (k) (k+1) (k) (k)
`i,k = (k)
, ai,j = ai,j − `i,k ak,j et bi = bi − `i,k bk
ak,k
alors le système devient :


 a1,1 x1 + a1,2 x2 + . . . + . . . + . . . + a1,n xn = b1 ,





 (2) (2) (2)

 0 x1 + a2,2 x2 + . . . + . . . + . . . + a2,n xn = b2 ,





 (3) (3) (3)
(S (k+1)
): 0 x1 + 0 x2 + a3,3 x3 + . . . + . . . + a3,n xn = b3 ,
 .. .


 . = ..

 (k) (k) (k)

 0 x1 + 0 x2 + . . . + ak,k xk + . . . + ak,n xn = bk ,


 .. .


 . = ..

 0 x + 0 x + . . . + a(k) x + . . . + a(k) x (k)
1 2 n,k k n,n n = bn ,

Formulation matricielle : Posons

 
1 0 ... ... ... ... ... 0
 .. 
 0 1 0 . 
 . . . . . 
 .. . . . . .. .. 
 
 . .. ..  (k)
 .. 0 1 . .  ai,k

Lk =  .  avec `i,k = , k+1≤i≤n
. . .  (k)
 .. .. −`k+1,k 1 .. ..  ak,k
 . .. .. 
 . . . . . .. 
 . . . 0 . . . 
 . .. .. .. . . . . 
 .. . . . . . 0 
0 ... 0 −`n,k 0 ... 0 1

A(1) = A, b(1) = b
A(2) = L1 A, ... , A(k+1) = Lk A(k) = Lk . . . L2 L1 A,
le système (S (k+1) ) équivalent à

A(k+1) x = Lk . . . L2 L1 A(1) x = Lk . . . L2 L1 b(1) = b(k+1)

D’où on obtient le système

Ax = b ⇔ A(k+1) x = b(k+1)
Au bout de (n − 1)-étape, on aboutit au système suivant

A(1) = A, b(1) = b

A(2) = L1 A, ... , A(n) = Ln−1 A(n−1) = Ln−1 . . . L2 L1 A,

le système (S (n) ) équivalent à

A(n) x = Ln−1 . . . L2 L1 A(1) x = Ln−1 . . . L2 L1 b(1) = b(n)

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

D’où on obtient le système

Ax = b ⇔ A(n) x = b(n)
avec A(n) est une matrice triangulaire supérieure.
Pour résoudre le système Ax = b, on résout le système équivalent Rx = c avec
R = A(n) = Ln−1 Ln−2 . . . L2 L1 A = M A où M = Ln−1 Ln−2 . . . L2 L1

et
c = b(n) = Ln−1 Ln−2 . . . L2 L1 b = M b où M = Ln−1 Ln−2 . . . L2 L1 .
Dans la pratique : On ne calcule pas le produit Ln−1 Ln−2 . . . L2 L1 mais plutot R = Ln−1 Ln−2 . . . L2 L1 .
Le calcul de R se fait par étape (où bien selon plusieurs algorithmes).
On ne stocke que la matrice A (n2 éléments de type réel) à la fin de la triangularisation on n’aura plus besoin
de la matrice A, mais, plutot on travaillera sur la matrice R. Donc A sera stochée dans la partie mémoire
resrvée pour le stockage de A.
L’algorithme (étapes de calcul) est donné dans un language naturel. Pour pouvoir le mettre en œuvre sur un
ordinateur il faut le traduire en un language de programmation tel que : Matlab, C, C++, Fortron, Java, . . ..
1. Algorithme 1 :
Pour chaque étape k : k = 1, . . . , n − 1

Lk A(k) : A(k+1) ← Lk A(k) ,
on calcule
Lk b(k) : b(k+1) ← Lk b(k) ,
fin pour.
Pour pouvoir traduire l’algorithme 1 en language de programmation, on doit l’affiner, c’est-à-dire
l’écrire à l’aide des opérations élémentaires, soit l’algorithme 2 obtenu après avoir affiné l’algorithme 1.
2. Algorithme 2 : On considère les deux phrases mathématiques suivantes :
(k)
(k+1) (k) ai,k (k)
(∗) ai,j ← ai,j − (k)
.ak,j
ak,k
(k)
(k+1) (k) ai,k (k)
(∗∗) bi ← bi − (k)
.bk
ak,k

Po u r k = 1 . . . n −1 , f a i r e
Po u r c h a q u e l i g n e i ( k< i <= n ) , f a i r e
p o u r c h a q u e i n d i c e j de l a l i g n e i , f a i r e
é c r i r e i c i l a phrase mathématique ( ∗ )
Fin pour j
é c r i r e i c i l a phrase mathématique (∗∗)
Fin pour i
Fin pour k
en réalité, les instructions
(k) (k)
(k+1) (k) ai,k (k) (k+1) (k) ai,k (k)
ai,j ← ai,j − .a
(k) k,j
et bi ← bi − (k)
.bk
ak,k ak,k
seront remplacées par les instructions (affectation)
ai,k ai,k
ai,j ← ai,j − .ak,j et bi ← bi − .bk
ak,k ak,k
2.1. MÉTHODES DIRECTES 21

(k)
(k+1) (k) ai,k (k)
Le signe ← où bien := veut dire que x = ai,j sera remplacé par y = ai,j − (k) .ak,j Autrement
ak,k
dit : dans la zone mémoire qui représente x on met la valeur de y.
Po u r k = 1 . . . n −1 , f a i r e
Po u r c h a q u e l i g n e i ( k< i <= n ) , f a i r e
p o u r c h a q u e i n d i c e j de l a l i g n e i , f a i r e
a [ i , j ] : = a [ i , j ]− a [ i , k ] / a [ k , k ] ∗ a [ k , j ]
Fin pour j
b [ i ] : = b [ i ]− a [ i , k ] / a [ k , k ] ∗ b [ k ]
Fin pour i
Fin pour k

3. Algorithme 3 :
Fo r k = 1 . . n −1 ,
Fo r i =k + 1 . . n ,
Fo r j =k + 1 . . n ,
a [ i , j ] : = a [ i , j ]− a [ i , k ] / a [ k , k ] ∗ a [ k , j ]
End j
b [ i ] : = b [ i ]− a [ i , k ] / a [ k , k ] ∗ b [ k ]
End i
End k

2.1.4 Factorisation LU d’une matrice

(k)
Supposons qu’à chaque étape k de la méthode de Gauss ak,k 6= 0. Dans ce cas, on a montré qu’il existe
une matrice M telle que :
MA = R

où R est une matrice triangulaire supérieur et M = Ln−1 Ln−2 . . . L2 L1 avec

 
1
0 ... ... ... ... ... 0
 .. .. 
 0 . 0 . 
 . . .. .. 
 .. . . 1 . . 
 
 . .. .. 
 .. 0 1 . . 

Lk =  . 
.. .. .. 
 .. . −`k+1,k 1 . . 
 . .. .. 
 . . . . . .. 
 . . . 0 . . . 
 . . . .. .. .. 
 .. .. .. . . . 0 
0 ... 0 −`n,k 0 ... 0 1

Propriété 2.1.1 (Propriété de la matrice M)

Qn−1
1. M est inversible puisque det(M) = k=1 det(Lk ) = 1 car det(Lk ) = 1 6= 0.
22 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

2. M −1 est une matrice triangulaire inférieure et

 
1 0 ... ... ... ... ... 0
 .. ..

 `2,1 . 0 .
 . ..
 ` .. ... ..
. .
 3,1 
 . . .. ..
 .. .. 1 .  .
−1
M = .  
.. 
..
 .. `k+1,k 1 . 
.
 . ..  ..
 . .. .. 
 . . `k+2,k+1 . .  .
 . . .. .. .. 
 .. .. . . . 0 
`n,1 . . . . . . `n,k `n,k+1 . . . `n,n−1 1

En effet, M −1 = L−1 −1 −1 −1
1 L2 . . . Ln−2 Ln−1 avec
 
1 0 ... ... ... ... ... 0
 . .. 
 0 .. 0 . 
 . . .. 
 .. . . . . . ..
. . 
 
 . .. .. 
 .. 0 1 . . 
−1
Lk =  . 
.. .. .. 
 .. . `k+1,k 1 . . 
 . .. .. 
 .. . . . . .. 
 . . 0 . . . 
 . . . .. .. .. 
 .. .. .. . . . 0 
0 ... 0 `n,k 0 ... 0 1

de la relation M A = R on déduit que A = M −1 R.

Posons L = M −1 Lower : triangulaire inférieure
et R = U Upper : triangulaire supérieure.

Finalement on a le théorème suivant :

(k)
Théorème 2.1.1 Si à chaque étape k de la méthode de Gauss on a le pivot ak,k 6= 0, alors il existe une
matrice triangulaire inférieur L et une matrice triangulaire supérieure U telles que A = L U qui s’appelle
la factorisation L U de la matrice A avec `i,i = 1.

D’une manière générale, on a le théorème suivant :

Théorème 2.1.2 Si A est inversible, alors il existe une matrice de permutation P telle que P A = L U où L
est une matrice triangulaire inférieur L et U une matrice triangulaire supérieure U avec L est une matrice
à diagonale unitée `i,i = 1.

Théorème 2.1.3 Soit A une matrice inversible.

A possède la factorisation L U, avec L est une matrice triangulaire inférieur L et U une matrice triangu-
laire supérieure U avec L est une matrice à diagonale unitée `i,i = 1, si et seulement si toutes les matrices
principales de A sont inversible.

Théorème 2.1.4 Si une A est inversible et possède la factorisation A = L U, où L est une matrice triangu-
laire inférieure et U une matrice triangulaire supérieure avec L est une matrice à diagonale unitée `i,i = 1,
alors cette factorisation est unique.
2.1. MÉTHODES DIRECTES 23

Interprétation : La méthode de Gauss utilisée pour résoudre un système régulier Ax = b (A inversible)

consiste à factoriser la matrice P A (P est une matrice de permutation bien choisie) en un produit L U avec
L est une matrice triangulaire inférieure et U une matrice triangulaire supérieure avec L est une matrice à
diagonale unitée `i,i = 1, en suite résoudre les systèmes suivants :

Ly = P b,
Ux = y

2.1.5 Factorisation de Cholesky où bien factoristion BB T

Définition 2.1.1 Soit A une matrice inversible

Une factorisation régulière de Cholesky de A est une factorisation A = BB T où B est une matrice trian-
gulaire inférieure régulière.
*Si les coefficients diagonaux de L sont positifs, on dit que l’on a une factorisation positive de Cholesky.

Lemme 2.1.1 Si A est une matrice symétrique définie positive, alors A possède la factorisation L U.

Démonstration. Montrons que toutes les matrices principales d’ordre 1 ≤ k ≤ n sont inversible
 
A11 A12
A=
A21 A22

décomposition en bloc A11 est une matrice principale d’ordre k.

Montrons
 que y T A11 y > 0 pour tout y 6= 0 avec y = (y1 , . . . , yk )T . Posons x = (x1 , . . . , xn )T avec
xi = yi , 1 ≤ i ≤ k
xi = 0, i > k
Alors xT Ax = y T A11 y et comme A est définie positive, on déduit que y T A11 y > 0 ce qui montre que A11
est définie positive et A11 est symétrique puisque A l’est.
D’où on a la factorisation L U de A. 

Théorème 2.1.5 Soit A une matrice régulière (inversible).

A possède une factorisation régulière de Cholesky A = B B T si et seulement si A est symétrique définie
positive.
*Dans ce cas, elle possède une factorisation positive unique.

Démonstration. ⇒c Supposons que A = B B T . Soit x ∈ Rn , x 6= 0 :

xT Ax = xT B B T x = (B T x)T B T x = (B T x, B T x) = kB T xk ≥ 0

comme x 6= 0 et B est régulière alors LT x 6= 0, par conséquent on a

(B T x)T B T x = (B T x, B T x) = kB T xk > 0

donc A est symétrique et définie-positive.

⇐c Supposons que A est symétrique et définie-positive, alors les valeurs propres µii > 0 pour 1 ≤ i ≤ n.
Soit D = diag(µ11 , . . . , µnn ) et R = D −1 U triangulaire supérieure à diagonale unité, donc on déduit que

A = LDR
24 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

est une décomposition unique.

Comme A est symétrique, alors on a

L D R = A = AT = RT D T LT = RT D LT

puisque la factorisation L U est unique, on déduit que L = RT , par conséquent on a

A = L D LT
√ √
Posons Λ = diag( µ11 , . . . , µnn ), alors

A = L Λ Λ LT = BB T avec B = LΛ et bii > 0

et la factorisation B B T de A est unique. 

Méthode pratique pour calculer B

Le calcul de la factorisation de Cholesky A = BB T peut se faire par identification. Comme A est

symétrique, on refait l’identification de coefficients de la partie triangulaire inférieure de A, on a A =
(aij )1≤i,j≤n et B = (bij )1≤i,j≤n :
n
X
aij = bik bjk , pour i ≥ j avec (bij = 0 pour i < j)
k=1

Calcul de B par colonne :

– Pour j = 1 :

√
a11 = b211 ⇒ b11 = a11
a21
a21 = b21 b11 ⇒ b21 =√
a11
..
.
ai1
ai1 = bi1 b11 ⇒ bi1 = √
a11
..
.
an1
an1 = bn1 b11 ⇒ bn1 = √
a11

– Pour j = 2
q
a22 = b21 b21 + b22 b22 ⇒ b22 = a22 − b221
a32 − b31 b21
a32 = b31 b21 + b32 b22 ⇒ b32 =
b22
d’une manière générale, on trouve la relation

ai2 − bi1 b21

ai2 = bi1 b21 + bi2 b22 ⇒ bi2 =
b22
2.2. MÉTHODES ITÉRATIVES 25

– Pour j > 1, la j eme colonne est calculée à partir des colonnes 1, . . . , (j − 1) de la façon suivante :
v
u j−1
u X
bjj = t ajj − b2jk
k=1
j−1
X
aij − bjk bik
k=1
bij = pour i ≥ j + 1.
bjj

– Le calcul de bjj nécessite 2(j − 1) opérations élémentaires et une racine carrée.

– Le calcul de bij nécessite 2(j − 1) opérations élémentaires (+, −) et une division.
Donc le calcul de la colonne j nécessite :

2(n − j + 1)(j − 1) + n − j (+, −, /) et une racine carrée

Le coût total ou la complexité de la méthode est

n3 n2 5n
+ − (+, −, /) et n racines carrées.
3 2 6

2.2 Méthodes itératives

2.2.1 Généralités
Soit A ∈ K(n×n) et b ∈ Kn . On veut résoudre le système linéaire Ax = b où A est une matrice inversible.
Résoudre le système linéaire Ax = b par une méthode itérative consiste à construire une suite de vecteurs
(xk )k∈N où xk ∈ Kn tel que
lim xk = x = A−1 b.
k→+∞

Les méthodes itératives classiques sont de la forme :

 o
x , valeur initiale donnée
xk+1 = B xk + c

où B ∈ K(n×n) et c ∈ Kn sont donnés en fonction de A et b.

Définition 2.2.1 Une méthode itérative est dite convergente si pour toute valeur initiale xO , on a

lim xk = x = A−1 b.
k→+∞

Pratiquement : on s’arrête lorsque kxk − xk est suffisamment petit.

Définition 2.2.2 Une méthode itérative est dite consistante si la limite de la suite xk lorsqu’elle existe, est
solution du système linéaire Ax = b.

Définition 2.2.3 Une méthode itérative classique de la forme xk+1 = B xk + c est consistante si et seule-
ment si 
c = (I − B)A−1 b,
I − B inversible.
26 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Remarque 2.2.1 Une méthode itérative consistante n’est pas toujours convergente.
Exemple 2.2.1 Soit à résoudre le système Ax = b avec a = 21 I par la méthode itérative xk+1 = B xk + c
avec 
c = −b,
B = I + A.
– cette méthode est consistante car B = I + A ⇒ I − B = −A qui est inversible et c =
(I − B)A−1 b = −b.
– cette méthode n’est pas convergente car :
3
xk+1 − x = B(xk − x) = (xk − x)
2
Par itération du schéma on obtient
 k+1
k+1 3
x −x= (xO − x) qui diverge
2
Théorème 2.2.1 On considère une méthode itérative consistante xk+1 = B xk + c. Le système linéaire
Ax = b admet une solution unique si les assertions suivantes sont équivalentes :
1. %(B) < 1,
2. lim B k x = 0 pour tout vecteur x,
k→+∞
3. La méthode itérative est convergente.
Démonstration. 1) ⇒ 2) : %(B) < 1 ⇒ il existe au moins une norme matricielle telle que kBk < 1, car
%(B) = inf{kBk / k · k est une norme matricielle induite}
⇒ limk→+∞ kB k k = 0 car kB k k ≤ kBkk ,
⇒ limk→+∞ B k = 0.
2) ⇒ 3) :
ek = xk − x = B xk−1 − x + c
= B xk−1 − x + (I − B)A−1 b
= B xk−1 − x + (I − B)x
= B(xk−1 − x)
= Bek−1
⇒ ek = B k e0 (e0 donnée)
⇒ lim ek = lim (B k e0 ) = 0 (car lim B k = 0
k→+∞ k→+inf ty k→+∞

d’où limk→+∞ xk = x,
en suite la méthode itérative est convergente.
3) ⇒ 1) : Pour tout x0 donnée, on a limk→+∞ xk − x = 0,
d’où pour tout x0 , limk→+∞ ek = 0 avec ek = xk − x.
Soit λ ∈ Sp(B) associée au vecteur propre v.
alors Bv = λv ⇒ B k v = λk v
pour x0 = x − v, on obtient
B k (x − x0 ) = λk (x − x0 ) ⇒ B k e0 = λk e0
or B k e0 = ek , on en déduit que limk→+∞ B k = 0 car limk→+∞ ek = 0,
comme B k e0 = λk e0 , on déduit que limk→+∞ λk = 0
ce qui prouve que |λ| < 1, ∀λ ∈ Sp(B), d’où %(B) < 1. 
2.2. MÉTHODES ITÉRATIVES 27

2.2.2 Comparaison des méthodes itératives

On considère la méthode itérative convergente définie par

 O
x , valeur initiale :donnée
xk+1 = B xk + c

posons ek = xk − x l’erreur à la k eme itérative.

On s’arrête lorsque kek k est suffisamment petit avec e0 = x0 − x est l’erreur initiale.
On obtient donc ek = B ek−1 et par conséquence ek = B ek−1 .

Méthode convergente ⇔ lim B k = 0 ⇔ lim ek = 0.

k→+∞ k→+∞

Soit k · k une norme matricielle : kek k ≤ kB k kke0 k

et donc
kek k k kB k yk
≤ kB k = sup .
ke0 k y6=0 kyk

Il est important de pouvoir estimer le nombre d’itérations (vitesse de convergence) nécessaires à l’obtention
d’une approximation acceptable de la solution.
kk
Problème : Trouver k tel que ke ke0 k
≤ ε : erreur permise. Il suffit que : kB k k ≤ ε, c’est-à-dire :

ln(kB k k) ≤ ln(ε)
k ln(kB k k1/k ) ≤ ln(ε)
ln(ε)
k ≥ , car kBk < 1
ln(kB k k1/k )

Définition 2.2.4 On considère la méthode itérative convergente suivante


xo , valeur initiale donnée
xk+1 = B xk + c

1. On appelle le taux moyen de convergence pour la k ième itérative, d’une méthode itérative conver-
gente, le nombre
Rk (B) = − ln(kB k k1/k ).

2. On appelle la vitesse de convergence, le nombre

v(B) = lim Rk (B).

k→+∞

Théorème 2.2.2 Soit B une matrice carrée. Pour toute norme matricielle, on a :

%(B) = lim kB k k1/k

k→+∞

et, par suite

v(B) = − ln(%(B)).

Démonstration.
28 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

– Montrons que ∀ε > 0, ∃kε tel que : ∀k ≥ kε , on a |kB k k1/k − %(B)| < ε ?
Soit λ ∈ Sp(B) une valeur propre associée à un vecteur propre x 6= 0, alors an a

Bx = λx implique B k x = λk x,

par passage à la norme, on déduit que

kλk xk ≤ kB k k kxk implique |λk | kxk ≤ kB k k kxk

d’où |λ| ≤ kB k k1/k , pour tout λ ∈ Sp(B), finalement %(B) ≤ kB k k1/k .

1
– Soit ε > 0, posons Bε = %(B)+ε B, alors on a

%(B)
%(Bε ) ≤ <1 ⇔ lim Bεk = 0
%(B) + ε k→+∞

par conséquent : limk→+∞ kBεk k = 0.

Ceci veut dire que :
∀ε0 > 0, ∃k0 , ∀k ≥ k0 , on a : kBεk k < ε0 .
Pour ε0 = 1, pour ε > 0, on a ∃k0 , ∀k ≥ k0 , on a : kBεk k < 1

kB k k
kBεk k = ⇒ kBεk k1/k < %(B) + ε
(%(B) + ε)k
d’où ∀ε > 0, ∃k0 , ∀k ≥ k0 , on a

%(B) ≤ kB k k1/k ≤ %(B) + ε.

Soit M1 : xk+1 = B1 xk + c1 , et M2 : xk+1 = B2 xk + c2 ,
deux méthodes itératives convergente proposées pour résoudre le système Ax = b.
Si %(B1 ) < %(B2 ), alors on a v(B1 ) > v(B2 )
donc on dit que la méthode M1 converge plus rapide que la méthode M2 .


2.2.3 Principales méthodes itératives classiques

Les méthodes itératives classiques sont basées sur la décomposition de A sous la forme A = M − N où
A est inversible, avec M ∈ K(n×n) est une matrice facille à inverser.

Ax = b ⇔ Mx = Nx + b
⇔ x = M −1 Nx + M −1 b
⇔ x = Bx + c avec B = M −1 N et c = M −1 b

On peut associer à la résolution du système A x = b, la méthode itérative suivante


xo , valeur initiale donnée
xk+1 = B xk + c

la méthode itérative ainsi construite est consistante car I − B = I − M −1 N = M −1 A est inversible, et

c = M −1 b = M −1 AA−1 b = (I − B)A−1 b.
Puisque cette méthode itérative est consistante alors pour qu’elle soit convergente, il suffit que

%(B) < 1.
2.2. MÉTHODES ITÉRATIVES 29

Dans ce cas, lim xk = x avec x est la solution du système linéaire Ax = b.

k→+∞

%(B) < 1 ⇔ %(M −1 N) < 1.

Méthode de Jacobi : Elle consiste à décomposer A sous la forme

A= D−E −F =M −N

avec M = D et N = E + F Où
     
a11 0 . . . 0 0 a12 . . . a1n 0 0 ... 0
 .. .. ..   .. .. ..   .. .. .. 
 0 . . .   0 . . .   a . . . 
D= . . .  , −F =  . .. .. , −E =  .21 .. .. 
 .. .. .. 0   .. . . a(n−1)n   .. . . 0 
0 . . . 0 ann 0 ... 0 0 an1 . . . an(n−1) 0

On suppose que aii 6= 0, pour tout 1 ≤ i ≤ n ; dans ce cas la matrice M serait inversible.
La méthode itérative s’écrit

xk+1 = B xk + c = D −1 (E + F )xk + D −1 b.

La matrice B = J = D −1 (E + F ) est appelée la matrice de Jacobi associée à A.

Proposition 2.2.1 La méthode de Jacobi converge si et seulement si %(J) < 1.

Comme A = D − E − F alors E + F = D − A.
La méthode itérative peut aussi s’écrire

xk+1 = D −1 (D − A)xk + D −1 b = (I − D −1 A)xk + D −1 b.

À partir de x ∈ Kn , on construit la suite (xk ) par ses composantes : les composantes de (xk ) sont données
par :  
n
X
1  
xk+1
i = bi − aij xkj 
aii j=1
j6=i

qu’on peut aussi écrire sous la forme :

n
!
1 X rik
xk+1
i − xki = bi − aij xkj =
aii j=1
aii

où rik est la ieme composante du vecteur r k défini par : r k = b − Axk appelé vecteur résidu à la k eme
itération.
Méthode de Gauss-Seidel : Elle consiste à prendre M = D − E (M est triangulaire inférieure) et N = F :

A= M −N =D−E−F
La matrice M est inversible si aii 6= 0, ∀1 ≤ i ≤ n.
La méthode itérative s’écrit

xk+1 = B xk + c = (D − E)−1 F xk + (D − E)−1 b

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

avec B = M −1 N = (D − E)−1 F et c = M −1 b = (D − E)−1 b.

d’où on obtient :
(D − E)xk+1 = F xk + b.
Supposons les composantes xk+1 k+1
1 , . . . , xi−1 sont calculées, alors la i
eme
composante xk+1
i est obtenue par
i−1 n
!
1 X X
xk+1
i = bi − aij xk+1
j − aij xk+1
j
aii j=1 j=i+1

On pose L1 = (D − E)−1 F , cette matrice s’appelle la matrice de Gauss-Seidel.

Proposition 2.2.2 La méthode de Gauss-Seidel converge si et seulement si %(L1 ) < 1.

2.2.4 Etude de la convergence

Définition 2.2.5 Soit A = (aij )1≤i,j≤n une matrice.
n
X
1. Une matrice A est dite à diagonale strictement dominante si on a |aii | > |aij |.
j=1
j6=i
n
X
2. Une matrice A est dite à diagonale dominante si on a |aii | ≥ |aij |.
j=1
j6=i

Théorème 2.2.3 Si A est à diagonale strictement dominante, alors la méthode de Jacobi converge.

Démonstration.A est à diagonale strictement dominante alors A est inversible.

Soit J = D −1 (E + F ) = I − D −1 A la matrice de Jacobi, alors
 a
ij

 − , si i 6= j
aii
Jij =


0, si i = j
n
X Xn n
aij 1 X
|Jij | = = |aij | < 1
j=1 j=1
aii |aii | j=1
j6=i j6=i

On déduit que !
n
X
max |Jij | = kJk∞ < 1
1≤i≤n
j=1

d’où %(J) < 1. D’où la méthode de Jacobi est convergente. 

Théorème 2.2.4 Soit A une matrice hermitienne (resp. symétrique) et définie positive.
Considérons la décomposition A = M − N avec M inversible.
Si la matrice M ∗ + N (resp. la matrice M T + N) est définie positive, alors %(M −1 N) < 1.
De plus, le schéma itératif classique xk+1 = B xk + c avecB = M −1 N et c = M −1 b est convergente.

Démonstration.
1. D’abord on a M ∗ + N = M ∗ + M − A.
2.2. MÉTHODES ITÉRATIVES 31

2. On a (M ∗ + N)∗ = (M ∗ + M − A)∗ = (M ∗ )∗ + M ∗ − A∗ = M + M ∗ − A = M ∗ + M − A, donc

M ∗ + N est hermitienne.
Pour montrer que %(M −1 N) < 1, on construit une norme vectorielle, notée k · kA , telle que la norme
matricielle induite vérifie kM −1 Nk < 1.
Comme A est définie positive, on pose kxkA = (x∗ Ax)1/2 , alors k · kA est bien une norme vectorielle.
Considérons la norme matricielle induite

kM −1 Nk = sup kM −1 NxkA = sup kx − M −1 AxkA ,

kxkA =1 kxkA =1

montrons que pour tout x ∈ Kn tel que kxkA = 1, on a kx − M −1 AxkA < 1 ?

Soit x ∈ Kn tel que kxkA = 1, posons M −1 Ax = y ce qui est équivalent à A x = M y.

kx − M −1 Axk2A = kx − yk2A
= (x − y)∗ A(x − y)
= x∗ Ax − x∗ Ay − y ∗ Ax + y ∗ Ay
= kxk2A − y ∗M ∗ y − y ∗My + y ∗ Ay
= 1 − y ∗(M ∗ + M − A)y
= 1 − y ∗(M ∗ + N)y

Pour x 6= 0, on a y 6= 0 (car A et M sont inversibles), alors

y ∗ (M ∗ + N)y > 0 car M ∗ + N est définie positive

comme kx − M −1 Axk2A > 0, alors 1 − y ∗ (M ∗ + N)y > 0, donc on déduit que

kx − M −1 Axk2A < 1, ∀x ∈ Kn , kxkA = 1

d’où on obtient kM −1 Nk < 1, ce qui prouve que %(M −1 N) < 1.



Corollaire 2.2.1 Soit A une matrice hermitienne (resp. A une matrice symétrique).
Si la matrice 2D − A est définie positive, alors la méthode de Jacobi converge si et seulement si la matrice
A est définie positive.

Cas des matrices tridiagonales : Soit A une matrice tridiagonale.

 
a1 b1 0 . . . 0
 . . 
 c1 a2 b2 . . .. 
 . . 
A=  0 c2 . . . . 0


 . . 
 .. . . . . . . . . bn 
0 . . . 0 c2 an

et les matrices E et F définient par

   
0 0 0 ... 0 0 b1 0 . . . 0
 . .. ..   . . . 
 c1 . . 0 . .   0 . . b2 . . .. 
 . ..   .. .. 
E=  0 c2 . . . 0  et F =
 0 0 . . 0 

 . .   . . 
 .. .. ... ..
. 0 
. .
 .. . . . . . . bn 
0 ... 0 cn 0 0 ... 0 0 0
32 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES

Lemme 2.2.1 Soit µ ∈ C∗ et Aµ = D − µE − µ1 F une matrice.

Alors, la matrice Aµ est inversible et on a det(Aµ ) = det(A).
Démonstration. On considère la matrice Qµ donnée par
 
µ 0 0 ... 0
 .. . 
 0 µ2 0 . .. 
 .. .. 
Qµ =  0 0 . . 0 ,
 µ 6= 0.
 . . . . 
 .. . . . . . . 0 
0 . . . 0 0 µn
On peut vérifier que Aµ = Qµ AQ−1
µ , d’où on obtient

det(Aµ ) = det(Qµ )det(A)det((Qµ )−1 ) = det(A).


Théorème 2.2.5 Soit A une matrice tridiagonale dont les éléments diagonaux sont non nul, alors
%(L1 ) = (%(J))2 .
La méthode de Gauss-Seidel et la méthode de Jacobi convergent ou divergent simultaniment.
Démonstration. Soit J = D −1 (E + F ) et L1 = (D − E)−1 F .
On calcule le polynôme caractéristique de J et de L1 :
PJ (λ) = det(D −1 (E + F ) − λIn )
= det(D −1 (E + F − λD))
= det(D −1 )det(E + F − λD)
= (−1)n det(D −1 )det(λD − E − F )
= (−1)n det(D −1 )det(λD + E + F ) (d’après le lemme : µ = −1)
= (−1)n det(λIn + D −1 (E + F ))
= (−1)n det(−(−λ)In + D −1 (E + F ))
= (−1)n PJ (−λ)
donc λ est une valeur propre de J si et selement si (−λ) est une valeur propre de J.
PL1 (λ) = det((D − E)−1 F − λIn )
= det((D − E)−1 )det(F − λD + λE)
= (−1)n det((D − E)−1 )det(λD − λE − F )
 
λ
n −1
= (−1) det((D − E) )det λD − E − µF , ∀µ ∈ C∗
µ
√
ceci d’après le lemme, on choisit µ = λ (où λ 6= 0), λ1/2 est un nombre complexe vérifiant (λ1/2 )2 = λ.
Pour λ 6= 0, on a


PL1 (λ) = (−1)n det((D − E)−1 )det λD − λ1/2 E − λ1/2 F


= (−1)n λn/2 det((D − E)−1 )det λ1/2 D − (E + F )


= λn/2 det((D − E)−1 )det(D)det D −1 (E + F ) − λ1/2 In


= λn/2 det J − λ1/2 In


= λn/2 PJ λ1/2
2.2. MÉTHODES ITÉRATIVES 33

car “det((D − E)−1 )det(D) = 1” puisque (D − E) est triangulaire inférieure et D est diagonale, ces deux
matrices ont la même diagonale, dans ce cas on a : det(D − E) = det(D).
ceci montre que λ est une valeur propre de L1 , avec λ 6= 0 alors {λ1/2 , −λ1/2 } ⊂ Sp(J).
Réciproquement : Si β 6= 0 est une valeur propre de J, alors β 2 est une valeur propre de L1 . D’où, on en
déduit que %(L1 ) = (%(J))2 . 

Corollaire 2.2.2 On a
%(J) < 1 ⇔ %(L1 ) < 1.
alors les méthodes de Jacobi et de Gauss-Seidel convergent où divergent simultanément.
Lorsque les deux méthodes convergent, alors la méthode de Gauss-Seidel converge deux fois plus vite que
la méthode de Jacobi : 2 ln(%(J)) = ln(%(L1 )) ce qui est équivalent : 2 v(J) = v(L1 ).
34 CHAPITRE 2. RÉSOLUTION DE SYSTÈMES LINÉAIRES
Chapitre 3

Interpolation et approximation polynômiale

3.1 Introduction

Ce chapitre traite de l’approximation d’une fonction dont on ne connaît les valeurs qu’en certains points.
Plus précisément, étant donné n + 1 couples (xi , yi ), le problème consiste à trouver une fonction Φ =
Φ(x) telle que Φ(xi ) = yi pour i = 0, . . . , m, où les yi sont des valeurs données.

Définition 3.1.1 On dit alors que Φ interpole {yi } aux nœuds {xi }.
On parle d’interpolation polynomiale quand Φ est un polynôme, d’approximation trigonométrique
quand Φ est un polynôme trigonométrique et d’interpolation polynomiale par morceaux (ou d’interpolation
par fonctions splines ou d’interpolation par fonctions à base radiales) si Φ est polynomiale par mor-
ceaux.

Les quantités yi peuvent, par exemple, représenter les valeurs aux nœuds xi d’une fonction f connue ana-
lytiquement ou des données expérimentales. Dans le premier cas, l’approximation a pour but de remplacer
f par une fonction plus simple en vue d’un calcul numérique d’intégrale ou de dérivée. Dans l’autre cas,
le but est d’avoir une représentation synthétique de données expérimentales dont le nombre peut être très
élevé. Nous étudions dans ce chapitre l’interpolation polynomiale, polynomiale par morceaux et les splines
paramétriques. Nous aborderons aussi les interpolations trigonométriques et interpolations basées sur les
polynômes orthogonaux.

Vision sur une fonction

Soient [a, b] un intervalle de R, S = (xi )0≤i≤n une subdivision de [a, b] et f : [a, b] → R une fonction
connue aux (n + 1) points xi (i = 0, . . . , n) de la subdivision S, c’est-à- dire qu’on connaît les valeurs

yi = f (xi ), pour i = 0, . . . , n.

Définition 3.1.2 On dit qu’un polynôme P de degrée inférieur ou égal à n (i.e., deg(P) ≤ n) interpole f
(ou encore interpole les valeurs y0 , . . . , yn aux (n + 1) points x0 , . . . , xn s’il vérifie les conditions d’inter-
polation suivantes :
P(xi ) = f (xi ), (où encore yi = P(xi ) i = 0, . . . , n
35
36 CHAPITRE 3. INTERPOLATION ET APPROXIMATION POLYNÔMIALE

3.2 Interpolation polynomiale

Considérons n + 1 couples (xi , yi ). Le problème d’interpolation consiste de trouver un polynôme Pm ∈
Pm ou Pm est l’espace des polynômes de degré m. Le polynôme Pm est appelé polynôme d’interpolation
ou polynôme interpolant, tel que
Pm (xi ) = a0 + a1 xi + . . . + am xm
i = yi , i = 0, . . . , n.
Les points xi sont appelés nœuds d’interpolation.

3.2.1 Interpolation polynomiale de Lagrange

Soit [a, b] un intervalle borné de R. Soit a = x1 < x2 < . . . < xn < xn+1 = b une subdivision de
l’intervalle [a, b]. On se donne n + 1 réels yi .
Pour tout i = 1, . . . , n + 1, on appelle polynôme de Lagrange 1 d’indice i, le polynôme `i défini par
Yn
x − xj
`i (x) = .
j=0
xi − xj
j6=i

Pour tout 0 ≤ i ≤ n, on a `i ∈ Pn et que le système {`0 , `1 , . . . , `n } forme une base de l’espace Pn . On

l’appelle base de Lagrange de Pn .
Théorème 3.2.1 Etant donné n + 1 points distincts x0 , . . . , xn et n + 1 valeurs correspondantes y0 , . . . , yn ,
alors il existe un unique polynôme Pn ∈ Pn tel que Pn (xi ) = yi pour i = 0, . . . , n.
Dé[Link] prouver l’existence, on va construire explicitement Pn . Posons
Yn
x − xj
`i ∈ Pn : `i (x) = , i = 0, . . . , n
j=0
xi − xj
j6=i

{`0 , `1 , . . . , `n } est une base de l’espace Pn , alors on a

n
X
Pn (x) = ci `i (x),
i=0

d’où n
X
Pn (xi ) = cj `j (xi ) = yi , i = 0, . . . , n
j=0

Il est facile de voir que 

1, si i = j,
`j (xi ) = δij =
0, si i 6= j.
on en déduit immédiatement que ci = yi pour tout i = 0, . . . , n.
Par conséquent, le polynôme d’interpolation existe et s’écrit sous la forme suivante
n
X
Pn (x) = yi `i (x). (0.1)
i=0

Pour montrer l’unicité, supposons qu’il existe un autre polynôme Ψm de degré m ≤ n tel que Ψm (xi ) = yi
pour i = 0, . . . , n. La différence Pn − Ψm s’annule alors en n + 1 points distincts xi , elle est donc nulle.
Ainsi, Pn = Ψm . 
1. [Link] est un mathématicien Franco-Italien(1736 − 1813)
3.3. DÉTERMINATION DU POLYNÔME D’INTERPOLATION 37

Corollaire 3.2.1 On a n
X ωn+1 (x)
Pn (x) = 0
yi
i=0
(x − xi )ωn+1 (xi )
où ωn+1 est le polynôme nodal de degré n + 1 défini par
n
Y
ωn+1 (x) = (x − xi ).
i=0

la relation (0.1) est appelée formule d’interpolation de Lagrange, et les polynômes `i (x) sont les poly-
nômes caractéristiques (de Lagrange).

Exemple 3.2.1 On considère l’intervalle [−1, 1] et une subdivision

x0 = −1 < x1 < . . . < xn = 1.

On choisit n = 4, on a x0 = −1, x1 = −0.5, x2 = 0, x3 = 0.5 et x4 = 1

Y4
x − xj x − x1 x − x2 x − x3 x − x4
`0 (x) = = ,
j=0
x0 − xj x0 − x1 x0 − x2 x0 − x3 x0 − x4
j6=1

Y4
x − xj x − x0 x − x2 x − x3 x − x4
`1 (x) = = ,
j=0
x1 − xj x1 − x0 x1 − x2 x1 − x3 x1 − x4
j6=2

Yn
x − xj x − x0 x − x1 x − x3 x − x4
`2 (x) = = ,
j=0
x2 − xj x2 − x0 x2 − x1 x2 − x3 x2 − x4
j6=3

Y4
x − xj x − x0 x − x1 x − x2 x − x4
`3 (x) = = ,
j=0
x3 − xj x3 − x0 x3 − x1 x3 − x2 x3 − x4
j6=4

Y4
x − xj x − x0 x − x1 x − x2 x − x3
`4 (x) = = .
j=0
x4 − xj x4 − x0 x4 − x1 x4 − x2 x4 − x3
j6=4

3.3 Détermination du polynôme d’interpolation

3.3.1 Cas où n = 2
Il s’agit de déterminer une polynôme d’interpolation P de degré n = 2, c’est-à-dire P ∈ P2 . On écrit

P(x) = a0 + a1 x + a2 x2 .

L’interpolation se fait en trois nœuds x0 , x1 et x2 , alors on écrit le système suivant :


 P(x0 ) = y0 ,
(E) : P(x1 ) = y1 ,

P(x2 ) = y2 .
38 CHAPITRE 3. INTERPOLATION ET APPROXIMATION POLYNÔMIALE

Représentation des polynômes de Lagrange

1.2

0.8

0.6

0.4

li(x)
0.2

−0.2

−0.4

−0.6
−1 −0.5 0 0.5 1
x

F IGURE 3.1 – Exemples de polynômes de Lagrange `2 , `3 et `4 .

y0 , y1 et y2 sont des données réelles et x0 , x1 et x2 sont des nœuds fixés.

Ainsi,le système (E) est équivalent au système suivant dons les inconnus sont a0 , a1 et a2 tles que

 a0 + a1 x0 + a2 x20 = y0 ,
(E) : a0 + a1 x1 + a2 x21 = y1 ,

a0 + a1 x2 + a2 x22 = y2 .
Ce système analytique est équivalent à un système matricielle qu’on peut écrire sous la forme suivante :
    
1 x0 x20 a0 y0
 1 x1 x21   a1  =  y1  .
1 x2 x22 a2 y2
| {z } | {z } | {z }
A x b

Ce système admet une unique solution si dét(A) 6= 0 (c’est-à-dire le discriminant ∆ 6= 0). On obtient un
système matricielle Ax = b à résoudre dont l’inconnu est le vecteur x, on a
   −1  
a0 1 x0 x20 y0
 a1  =  1 x1 x21   y1  .
a2 1 x2 x22 y2
Le discriminant ∆ = dét(A) qu’on calcule selon
1 x0 x20
dét(A) = 1 x1 x21 = (x2 − x1 )(x2 − x0 )(x1 − x0 ).
1 x2 x22
Une fois les coefficients a0 , a1 et a2 sont calculés, on a :
P(x) = a0 + a1 x + a2 x2
 
a0
= (1 x x2 )  a1 
a2
 −1  
1 x0 x20 y0
= (1 x x2 )  1 x1 x21   y1 
1 x2 x22 y2
 
y0
= (ξ0 ξ1 ξ2 )  y1 
y2
3.3. DÉTERMINATION DU POLYNÔME D’INTERPOLATION 39
 −1
1 x0 x20
avec (ξ0 ξ1 ξ2 ) = (1 x x2 )  1 x1 x21  . D’après le technique de transposition, on obtient le
1 x2 x22
système matricielle suivant :
   −1  
ξ0 1 1 1 1
 ξ1  =  x0 x1 x2   x 
ξ2 x20 x21 x22 x2

En résolvant ce système pour ontenir

(x − x1 )(x − x2 )
ξ0 =
(x0 − x1 )(x0 − x2 )

(x − x0 )(x − x2 )
ξ1 =
(x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 )
ξ2 = .
(x2 − x0 )(x2 − x1 )
On pose `0 (x) = ξ0 , `1 (x) = ξ1 et `2 (x) = ξ2 qui sont respectivement les polynômes de Lagrange associés
à x0 , x1 et x2 . Dans ce cas, on a :

P(x) = a0 + a1 x + a2 x2 = y0 `0 (x) + y1 `1 (x) + y2 `2 (x).

Exemple 3.3.1 Pour des valeurs données, de la fonction x 7→ f (x), aux points x = x0 , x2 et x2 , on exprime
le polynôme d’interpolation par

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

P2 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Pour f (1) = 0, f (−1) = −3 et f (2) = 4 avec x0 = 1, x1 = −1 et x2 = 2 on obtient

(x + 1)(x − 2) (x − 1)(x − 2) (x − 1)(x + 1)

P2 (x) = ×0+ × (−3) + ×4
(1 + 1)(1 − 2) (−1 − 1)(−1 − 2) (2 − 1)(2 + 1)

Après simplification, on a
1

P2 (x) = 5x2 + 9x − 14 .
6
On a bien P2 (1) = 0, P2 (−1) = −3 et P2 (2) = 4.

3.3.2 Cas général

On cherche P ∈ Pn vérifiant le problème d’interpolation P(xi ) = yi pour i = 0, . . . , n, c’est à dire, on

cherche a0 , . . . , an tels que P(x) = a0 + a1 x + . . . + an xn et vérifiant les (n + 1) équations ci-dessous :


 a0 + a1 x0 + . . . + an xn0 = y0 ,

 a0 + a1 x1 + . . . + an xn = y1 ,
1
(E) : .. .


 . = ..
 a + a x + . . . + a xn = y .
0 1 n n n n
40 CHAPITRE 3. INTERPOLATION ET APPROXIMATION POLYNÔMIALE

La forme matricielle de ce système est :

    
1 x0 x20 . . . xn0 a0 y0
 1 x1 x2 . . . xn1   a1   y1 
 1    
 .. .. .. .   .. = .. 
 . . . . . . ..   .   . 
1 xn x2n . . . xnn an yn
et on peut facillement montrer que
1 x0 x20 . . . xn0
1 x1 x21 . . . xn1 Y
∆ = dét(A) = .. .. .. . = (xi − xj ).
. . . . . . .. 0≤i<j≤n
1 xn x2n . . . xnn
√
Exemple 3.3.2 Cherchons P le polynôme interpolant la fonction f (x) = x2 + 1 aux points 0, 1, 4 et 9.
En notant x0 = 0, x1 = 1, x2 = 4 et x3 = 9 et en utilisant la base de Lagrange, on peut affirmer que :
P(x) = f (x0 )`0 (x) + f (x1 )`1 (x) + f (x2 )`2 (x) + f (x3 )`3 (x).
√ √ √
On a f (x0 ) = 1, f (x1 ) = 2, f (x2 ) = 17, f (x3 ) = 82. On a bien
x − x1 x − x2 x − x3 1
`0 (x) = =− (x − 1)(x − 4)(x − 9),
x0 − x1 x0 − x2 x0 − x3 36
x − x0 x − x2 x − x3 1
`1 (x) = = x(x − 4)(x − 9),
x1 − x0 x1 − x2 x1 − x3 24
x − x0 x − x1 x − x3 1
`2 (x) = = − x(x − 1)(x − 9),
x2 − x0 x2 − x1 x2 − x3 60
x − x0 x − x1 x − x2 1
`3 (x) = = x(x − 1)(x − 4).
x3 − x0 x3 − x1 x3 − x2 360
Donc on obtient √ √ √
P(x) = `0 (x) + 2`1 (x) + 17`2 (x) + 82`3 (x)

3.4 Interpolation par les différences divisées

Dans cette section, nous discuterons une autre méthode d’interpolation polynômiale dans laquelle nous
utiliserons les différences divisées. Cette méthode est due à Isaac Newton.
On note par x0 , x1 , . . . , xn les (n + 1) points distincts et on désigne par Pn le polynôme d’interpolation de
la fonction f , qui est donné par l’expression suivante :
Pn (x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + . . . + an (x − x0 )(x − x1 ) . . . (x − xn−1 ).
Par la substitution x = x0 , x1 , . . . , xn dans l’expression explicite de Pn , on obtient le système suivant :
f (x0 ) = a0 ,
f (x1 ) = a0 + a1 (x1 − x0 ),
f (x2 ) = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 ),
..
.
f (xn ) = a0 + a1 (xn − x0 ) + a2 (xn − x0 )(xn − x1 ) + . . .
+an (xn − x0 )(xn − x1 ) . . . (xn − xn−1 ).
3.4. INTERPOLATION PAR LES DIFFÉRENCES DIVISÉES 41

Ce système d’équations nous permet de déterminer les coefficients a0 , a1 , . . . , an de façon unique. La

première équation détermine a0 = f (x0 ), la deuxième détermine a1 = f (x 1 )−a0
x1 −x0
si x1 − x0 6= 0, si
f (x2 )−[a0 +a1 (x2 −x0 )]
(x2 − x0 )(x2 − x1 ) 6= 0 alors a2 = (x2 −x0 )(x2 −x1 ) . Finallement, si on connait a0 , a1 , . . . , an−1 alors la
dernière équation du système nous donne an puisque (xn − x0 )(xn − x1 ) . . . (xn − xn−1 ) 6= 0.
Le polynôme Pn s’écrit donc de façon unique et la formule reste permanente. Mais, si on ajoute un autre
point xn+1 à la suite x0 , x1 , . . . , xn alors on construit un polynôme d’interpolation Pn+1 lié à la nouvelle
suite x0 , x1 , . . . , xn , xn+1 de la forme suivante :
Pn+1 (x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) + . . . + bn (x − x0 ) . . . (x − xn−1 )
+bn+1 (x − x0 ) . . . (x − xn−1 )(x − xn ).
En utilisant le même principe qu’avant, on trouve b0 = a0 , b1 = a1 , . . . , bn = an car par substitution des
points x = x0 , x1 , . . . , xn , xn+1 dans l’expression Pn+1 , on a :
f (x0 ) = b0 ,
f (x1 ) = b0 + b1 (x1 − x0 ),
f (x2 ) = b0 + b1 (x2 − x0 ) + b2 (x2 − x0 )(x2 − x1 ),
..
.
f (xn ) = b0 + b1 (xn − x0 ) + b2 (xn − x0 )(xn − x1 ) + . . .
+bn (xn − x0 )(xn − x1 ) . . . (xn − xn−1 ),
f (xn+1 ) = b0 + b1 (xn − x0 ) + b2 (xn − x0 )(xn − x1 ) + . . .
+bn+1 (xn+1 − x0 )(xn+1 − x1 ) . . . (xn+1 − xn ).
Par la résolution de ce système d’équations, on aboutit à ce que
bi = ai , ∀1 ≤ i ≤ n
En général, on a l’expression
aj = f [x0 , x1 , . . . , xj ]
avec la relation suivante : n
X f (xi )
f [x0 , x1 , . . . , xn ] = n .
Y
i=0
(xi − xj )
j=0
j6=i

Cette formule consiste à écrire les valeurs f (xi ) sous la forme d’une combinaison linéaire
Exemple 3.4.1 Pour n = 2, on a
f (x0 ) f (x1 ) f (x2 )
f [x0 , x1 , x2 ] = + +
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
Est-ce qu’on peut écrire l’expression f [x0 , x1 , x2 , x3 ] sous une forme plus ssimple ? Trouver α et β tels que
f [x0 , x1 , x2 , x3 ] = αf [x0 , x1 , x2 ] + βf [x1 , x2 , x3 ].
Par comparaison des deux termes, on obtient
1 1
α= et β=
x0 − x3 x3 − x0
f [x1 , x2 , x3 ] − f [x0 , x1 , x2 ]
f [x0 , x1 , x2 , x3 ] = .
x3 − x0
42 CHAPITRE 3. INTERPOLATION ET APPROXIMATION POLYNÔMIALE

La formule générale des différences divisées est donnée par l’expression récurrente :
f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ]
f [x0 , x1 , . . . , xn ] =
xn − x0
Nous avons écris f [xi ] au lieu de f (xi ). D’où, par exemple, on calcule f [x2 , x3 ] par

x0 f [x0 ]
f [x0 , x1 ]
x1 f [x1 ] f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f [x2 ] f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f [x3 ]
TABLE 3.1 – Schéma pour le calcul des différences divisées

f [x3 ] − f [x2 ]
f [x2 , x3 ] =
x3 − x2
et pour f [x1 , x2 , x3 ] par
f [x2 , x3 ] − f [x1 , x2 ]
f [x1 , x2 , x3 ] = .
x3 − x1
D’après l’expression générale du polynôme d’interpolation de Newton, on obtient
Pn (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) + . . .
+f [x0 , . . . , xn ](x − x0 )(x − x1 ) . . . (x − xn−1 ).
Exemple 3.4.2 En utilisant le principe des différences divisées pour obtenir le polynôme d’interpolation
tel que
f (1) = 0, f (−1) = −3, f (2) = 4
Comme l’ulistre le tableau, on remarque
x f (x) Différences divisées
1 0
−3 − 0 3
=
−1 − 1 2
7
3
− 23 5
−1 −3 =
2−1 6
4 − (−3) 7
=
2 − (−1) 3
2 4
TABLE 3.2 – Schéma pour le calcul des différences divisées

3 7 f [x1 , x2 ] − f [x0 , x1 ] 5
f [x0 , x1 ] = , f [x1 , x2 ] = , f [x0 , x1 , x2 ] = = .
2 3 x2 − x0 6
Ainsi, on obtient
P2 (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )
3 5
= 0 + (x − 1) + (x − 1)(x + 1)
2 6
1

= 5x2 + 9x − 14 .
6
3.5. ERREUR D’INTERPOLATION 43

3.5 Erreur d’interpolation

Dans cette section, nous évaluons l’erreur d’interpolation faite quand on remplace une fonction f donnée
par un polynôme Pn qui l’interpole aux nœuds x0 , x1 , . . . , xn .

Théorème 3.5.1 Soient x0 , x1 , . . . , xn , (n+1) nœuds distincts et soit x un point appartenant au domaine de
définition de f . On suppose que f esst de classe C n+1 sur le plus petit intervalle Ix contenant x0 , x1 , . . . , xn
et x. L’erreur d’interpolation au point x est donnée par
n
f (n+1) (ξ) Y
En (x) = f (x) − Pn (x) = (x − xi ),
(n + 1)! i=0

où ξ ∈ Ix .

Démonstration. Comme f (xi ) = Pn (xi ) pour i = 0, 1, . . . , n, alors il existe une fonction réelle φ telle que
n
Y
f (x) − Pn (x) = φ(x) (x − xi ),
i=0
Q
Pour x fixé et x 6= xi (i = 0, 1, . . . , n), posons g(t) = f (t) − Pn (t) − c ni=0 (x − xi ) où c est un paramètre
réel choisi tel que g(x) = 0. Dans ce cas, la condition g(x) = 0 entraîne

f (x) − Pn (x)
c= n .
Y
(x − xi )
i=0

Maintenant, remarquons que la fonction g est de classe C n+1 sur [a, b] et de plus

g(x) = 0,
g(xi ) = 0, pour i = 0, . . . , n,

c’est à dire que g admet n + 2 racines distinctes deux à deux et est de de classe C n+1 sur [a, b]. Cela entraîne
que g 00 admet n + 1 racines distinctes deux à deux et est de de classe C n sur [a, b].
Et en réitérant, on a g 00 admet n racines distinctes deux à deux et est de de classe C n−1 sur [a, b] et finalement,
on arrive à vérifier que g (n+1) admet une racine sur ]a, b[. Ainsi,

∃ξx ∈]a, b[, tel que g (n+1) (ξx ) = 0.

(n+1) (n+1)
Or, g (n+1) (t) = f (n+1) (t)−Pn (t)−c(n+1)!, et puisque Pn (t) = 0 (car deg(Pn ) ≤ n) et g (n+1) (ξx ) =
0, alors on vérifie que
n
f (n+1) (ξx ) Y
c= (x − xi ),
(n + 1)! i=0
c’est-à- dire que
n
f (n+1) (ξ) Y
f (x) − Pn (x) = (x − xi ).
(n + 1)! i=0

44 CHAPITRE 3. INTERPOLATION ET APPROXIMATION POLYNÔMIALE

Corollaire 3.5.1 Dans les conditions du théorème ci-dessus, on

Mn+1 (b − a)n+1
sup |f (x) − Pn (x)| ≤ ,
x∈[a,b] (n + 1)!

où Mn+1 = supx∈[a,b] |f (n+1) (x)|.

Dé[Link] suffit d’utiliser la technique de majoration. 

Remarque 3.5.1 La majoration de l’erreur d’interpolation donnée ci-dessus peut laisser croire que si le
nombre d’abscisses d’in- terpolation n est grand, alors le polynôme d’interpolation Pn de f aux points
d’interpolation x0 , x1 , . . . , xn tend vers f . En fait, on n’a pas nécessairement lim f (x) − Pn (x) = 0 pour
n→+∞
tout x ∈ [a, b] car cette limite dépend aussi de la façon dont la quantité Mn se comporte lorsque la valeur
de n devient large. L’exemple ci-dessous permet d’illustrer cette remarque.

Exemple 3.5.1 (Phénomène de Runge)

Supposons qu’on approche la fonction suivante
1
f (x) = , a ≤ x ≤ b.
1 + 9x2
Chapitre 4

Intégration et dérivation numérique

4.1 Introduction et outils de base

Nous présentonsdans ce chapitre les méthodes les plus couramment utilisées pour l’intégration numé-
rique. Bien que nous nous limitons essentiellement aux intégrales sur des intervalles bornés. Eventuelle-
ment, nous oborderons des extensions aux intervalles non bornés.

Théorème 4.1.1 (Première formule de la moyenne)

Soient u et v deux fonctions continues sur [a, b] telles que u est de signe constant dans [a, b]. Alors
Z b Z b
∃η ∈]a, b[ tel que u(x)v(x)dx = v(η) u(x)dx.
a a

4.2 Formule de quadrature

Soit f une fonction réelle intégrable sur un intervalle [a, b]. Le calcul explicite de l’intégrale définie
Rb
I(f ) = a f (x)dx peut être difficile, voire impossible.

Définition 4.2.1 On appelle formule de quadrature ou formule d’intégration numérique toute formule
permettant de calculer une approximation de I(f ).

Définition 4.2.2 Soit I =]a, b[ un intervalle avec −∞ ≤ a < b ≤ +∞. On appelle formule (ou méthode)
de quadrature à n+1 points sur C(I) la donnée d’une forme linéaire ϕn définie sur C(I) par :
n
X
∀f ∈ C(I), ϕn = λn,k f (xn , k)
k=0

où n est un entier naturel non nul, (xn,k )0≤k≤n une suite de points deux à deux distincts dans l’intervalle I
et (yn,k )0≤k≤n une suite de réels non tous nuls.
Rb
Si π ∈ C(I, R+∗ ) et f ∈ C(I) telle que a f (x)π(x)dx soit convergente, on cherche des formules de
quadrature, avec des coefficients xn,k et λn,k indépendants de f, permettant d’approximer une telle intégrale.
On écrira : Z b n
X
f (x)π(x)dx ' ϕn (f ) = λn,k f (xn , k)
a k=0
45
46 CHAPITRE 4. INTÉGRATION ET DÉRIVATION NUMÉRIQUE

Définition 4.2.3 Une formule de qudrature à n+1 points sur C(I) est dite d’ordre p si elle est exacte sur
Rp [x] et iinexacte pour au moins un polynôme de Rp+1[x].

Définition 4.2.4 On dit qu’une méthode de quadrature définie par une suite de formes linéaires (ϕn )n est
convergente si pour toute fonction f ∈ C(I) la suite réelles (ϕn (f ))n converge vers ϕ(f ).

Une possibilité consiste à remplacer f par une approximation fn , où n est un entier positif, et calculer
I(fn ) au lieu de I(f ). En posant In (f ) = I(fn ), on a
Z b
In (f ) = fn (x)dx, n ≥ 1. (0.1)
a

La dépendance par rapport aux extrémités a et b sera toujours sous-entendue.

On écrira donc In (f ) au lieu de In (f ; a, b).

Si f ∈ C 0 ([a, b]), l’erreur de quatrature En (f ) = I(f ) − In (f ) satisfait

Z b
|En (f )| ≤ |f (x) − fn (x)|dx ≤ (b − a) sup |f (x) − fn (x)| = (b − a)kf − fn k∞ .
a x∈[a,b]

Donc, si pour un certain n, kf − fn k∞ < ε, alors |En (f )| ≤ ε(b − a).

L’approximation fn doit être facilement intégrable, ce qui est le cas si, par exemple, fn ∈ Pn . Une
méthode naturelle consiste à prendre fn = Πn f , le polynôme d’interpolation de Lagrange de f sur un
ensemble de n + 1 nœuds distincts {xi , i = 0, . . . , n}. Ainsi, on définit de (0.1) que
n
X Z b
In (f ) = f (xi ) `i (x)dx, n ≥ 1. (0.2)
i=0 a

où `i est le polynôme caractéristique de Lagrange de degré n associé au nœuds xi . On note que (0.2) est un
cas particulier de la formule de quadrature suivante :
n
X
In (f ) = αi f (xi ), n ≥ 1. (0.3)
i=0
Rb
où les coefficients αi de la combinaison linéaire sont donnés par αi = a
`i (x)dx.

Le système {(f (x1 ), α1 ), . . . , (f (xn ), αn )} est un système pondéré de somme In (f ). On dit aussi que
la formule (0.3) est une somme pondérée des valeurs de f aux points xi , pour i = 0, . . . , n. On dit que
ces points sont les nœuds de la formule de quadrature et que les nombres αi ∈ R sont ses coefficients ou
encore ses poids. Les poids est les nœuds dépendent en général de n.

Définition 4.2.5 On appelle le degré d’exactitude d’une formule de quadrature, le plus grand entier r ≥ 1
pour lequel
In (f ) = I(f ), ∀fn ∈ Pr

Exercice 1
Montrer que l’ordre maximum d’une formule de quadrature à n+1 points est 2n+1.
4.3. QUADRATURES INTERPOLATOIRES 47

4.3 Quadratures interpolatoires

4.4 La méthode des rectangles à gauche

Pour I = [a, b], xn,k = a + k b−a
n
et λn,k = b−a
n
, on a la formule des rectangles à gauche :
n−1
b−aX
Rn (f ) = f (xn,k )
n k=0

La méthode des rectangles à gauche est d’ordre 0 et convergente.

Theorem 4.4.1 Soit f une fonction de classe C 1 sur [a,b] et n un vecteur naturel non nul. Avec les notations
qui précèdent, il existe un réel cn ∈ [a, b] tel que :
Z b
(b − a)2 0
f (x)dx − Rn (f ) = f (cn )
a 2n
et on a : Z b
kf 0 k∞ (b − a)2
f (x)dx − Rn (f ) ≤
a 2n

4.4.1 Formule du rectangle ou du point milieu

Cette formule est obtenue en remplaçant f par une constante égale à la valeur de f au milieu de [a, b],
ce qui donne
 
a+b
I0 (f ) = (b − a)f . (1.1)
2

Le poids est donc α0 = b − a et le nœud x0 = (a + b)/2.

Si f ∈ C 2 ([a, b]), l’erreur de quadrature est

h3 00 b−a
E0 (f ) = f (ξ), h= , (1.2)
3 2
où ξ ∈]a, b[. En effet, le développement de Taylor au second ordre de f en c = (a + b)/2 s’écrit

f (x) = f (c) + f 0 (c)(x − c) + f 00 (η(x))(x − c)2 /2

en intégrant sur [a, b] et on obtient

Z b
h3 00
f (x)dx = (b − a)f (c) + f (ξ).
a 3

Exemple 4.4.1 On considère la fonction f : x 7→ f (x) = sin(x), on veut intégrer cette fonction sur
l’intervalle [1, 1.2].
La formule du point milieu implique
Z 1.2
sin(x)dx = (1.2 − 1) sin((1 + 1.2)/2) ≈ 0.2 × 0.8912 = 0.1782.
1
48 CHAPITRE 4. INTÉGRATION ET DÉRIVATION NUMÉRIQUE

Par contre la valeur théorique est

Z 1.2
sin(x)dx = − cos(1.2) + cos(1) = 0.1779.
1

L’erreur relative est

|0.1779 − 0.1782|
.100% = 0.17%
0.1779
Théorème 4.4.1 (Formule de la moyenne discrète)
Soit u ∈ C 0 ([a, b]), soient xj les (s + 1) points de [a, b] et δj les (s + 1) constantes, toutes de même signe.
Alors, il existe η ∈]a, b[ tel que
Xs X s
δj u(xj ) = u(η) δj .
j=0 j=0

Dé[Link] um = minx∈[a,b] u(x) = u(xm ) et uM = maxx∈[a,b] u(x) = u(xM ), où xm et xM sont

deux points de [a, b]. Alors
s
X s
X s
X
um δj ≤ δj u(xj ) ≤ uM δj . (1.3)
j=0 j=0 j=0
Ps P
On pose σs = j=0 δj u(xj ) et on considère la fonction continue U(x) = u(x) sj=0 δj . D’après (1.3), on
a U(xm ) ≤ σs ≤ U(xM ). Le théorème de la moyenne implique l’existence d’un point η ∈]a, b[ tel que
U(η) = σs . d’où le résultat 

*Formule composite du rectangle : Supposons maintenant qu’on approche l’intégrale I(f ) en remplaçant
f par son interpolation polynomiale composite de degré 0 sur [a, b], construite sur m sous-intervalles de
largeur h = (b − a)/m, avec m ≥ 1. En introduisant les nœuds de quadrature xk = a + (2k + 1)h/2, pour
k = 0, . . . , m − 1, on obtient la formule composite du point milieu :
m−1
X
Im = h f (xk ), m ≥ 1. (1.4)
k=0

Si f ∈ C 2 ([a, b]), l’erreur de quadrature Em (f ) = I(f ) − Im (f ) est donnée par

b − a 2 00
Em (f ) = h f (ξ), (1.5)
24
où ξ ∈]a, b[. On déduit de (1.5) que (1.4) a un degré d’exactitude égal à 1 ; on peut montrer (1.5) en utilisant
(1.2) et la linéarité de l’intégration.
En effet, pour k = 0, . . . , m − 1 et ξk ∈]a + kh, a + (k + 1)h[,
m−1
X m−1
X
00 3 h2 b − a
Em = f (ξk )(h/2) /3 = f 00 (ξk ) ,
24 m
k=0 k=0

d’après la formule de la moyenne discrète, pour u = f 00 et δj = 1 pour j = 0, . . . , m − 1, on obtient

b − a 2 00
Em = h f (ξ).
24
D’où
Z b m−1
X b − a 2 00
f (x)dx = h f (a + (2k + 1)h/2) + h f (ξ), où ξ ∈]a, b[.
a 24
k=0
4.4. LA MÉTHODE DES RECTANGLES À GAUCHE 49

4.4.2 Formule du trapèze

Cette formule est obtenue en remplaçant f par Π1 f , son polynôme d’interpolation de Lagrange de
degrée 1 aux noeuds x0 = a et x1 = b. Les noeuds de la formule de quadrature sont alors x0 = a, x1 = b et
ses poids α0 = α1 = (b − a)/2 :
Z b
b−a
I1 (f ) = f (x)dx = [f (a) + f (b)].
a 2

En effet, soit A(a, f (a)) et B(b, f (b)) les deux points dont les abscisses respectifs a et b. Les points
M(x, y) du segment [AB] ont une abscisse x ∈ [a, b] et une ordonnée

f (b) − f (a)
y(x) = f (a) + (x − a)
b−a
L’approximation par la méthode des trapèze implique
Z b Z b
f (x)dx ≈ y(x)dx
a a

Analytiquement, on obtient
Z b Z b 
f (b) − f (a) b−a
y(x)dx = f (a) + (x − a) dx = (f (a) + f (b)) .
a a b−a 2

Cette formule est exacte pour les polynômes de degré 1 mais pas pour x2 , elle est donc d’ordre 1.
Si f ∈ C 2 ([a, b]), l’erreur de quadrature est donnée par

h3 00
E1 (f ) = − f (ξ), h = b − a,
12
où ξ est un point de l’intervalle d’intégration. C’est-à-dire que

I(f ) = I1 (f ) + E1 (f ).

En effet, l’expression de l’erreur d’interpolation implique

Z b Z b
1
E1 (f ) = (f (x) − Π1 f (x))dx = − f 00 (ξx )(x − a)(b − x)dx.
a 2 a

Comme ω2 (x) = (x − a)(x − b) < 0 sur ]a, b[, alors d’après le théorème de la moyenne on a
Z b
1 (b − a)3
E1 (f ) = f 00 (ξ) ω2 (x)dx = f 00 (ξ) ,
2 a 12

pour un ξ ∈]a, b[, d’où le résultat.♦

Exemple 4.4.2 On considère la fonction f : x 7→ f (x) = ex , on veut intégrer cette fonction sur l’intervalle
[1, 2]. La valeur théorique de cette intégrale est
Z 2
ex dx = e2 − e1 = e(e − 1) ≈ 4.6708.
1
50 CHAPITRE 4. INTÉGRATION ET DÉRIVATION NUMÉRIQUE

Maintenant, on calcule la valeur approximative de l’intégrale où bien l’intégration numérique en utilisant

la méthode du trapèze. Soit A(1, e) et B(2, e2 ) les deux points du plan ayant pour abscisses respectifs 1 et
2, les points du segmant [AB] ont une abscisse x ∈ [1, 2] et une ordonnées y(x) = e[1 + (e − 1)(x − 1)].
On a bien y(1) = e et y(2) = e2 .
Z 2 Z 2
y(x)dx = e[1 + (e − 1)(x − 1)]dx = (2 − 1)/2(e1 + e2 ) ≈ 5.0537
1 1

D’où l’erreur relative est

|4.6708 − 5.0537|
.100% = 7.4%
4.6708

Exemple 4.4.3 On considère la fonction f : x 7→ f (x) = sin(x), on veut intégrer cette fonction sur
l’intervalle [1, 1.2].
La formule des trapèze implique
Z 1.2
sin(x)dx = (1.2 − 1)/2(sin(1) + sin(1.2)) ≈ 0.1 × 1.7735 = 0.1774.
1

Par contre la valeur théorique est

Z 1.2
sin(x)dx = − cos(1.2) + cos(1) = 0.1779.
1

L’erreur relative est

|0.1779 − 0.1774|
.100% = 0.28%
0.1779

On peut observer que par la relation de Chasles pour les intégrales, on a [a, b] = [a, c]∪[c, b] où c = (a+b)/2,
et on a Z Z Z b c b
f (x)dx = f (x)dx + f (x)dx.
a a c

D’après la formule de trapèze pour deux points, on obtient

Z b 
  
c−a b−c
f (x)dx ≈ (f (a) + f (b)) + (f (b) + f (c))
a 2 2

or, c − a = b − c = (b − a)/2, alors

Z b
b−a
f (x)dx ≈ [f (a) + 2f (c) + f (b)] .
a 4

*Formule composite d’intégration par la méthode de trapèze : Soit [a, b] un intervalle de R et n ≥ 1 un

entier. Soit h = b−a
n
et on définit xi = a + ih pour i = 0, . . . , n. On considère la subdivision suivante

a = x0 < x1 < . . . < xn−1 < xn = b

Alors on a
Z b Z x1 Z x2 Z xn n−1 Z
X xi+1
f (x)dx = f (x)dx + f (x)dx + . . . + f (x)dx = f (x)dx
a x0 x1 xn−1 i=0 xi
4.4. LA MÉTHODE DES RECTANGLES À GAUCHE 51

D’après la formule des trapèzes pour deux points xi et xi+1 on a

Z xi+1
xi+1 − xi h
f (x)dx = [f (xi+1 ) + f (xi )] = [f (xi+1 ) + f (xi )]
xi 2 2
Alors
Z n−1 Z
X n−1
b xi+1
hX
f (x)dx = f (x)dx = [f (xi+1 ) + f (xi )]
a i=0 xi 2 i=0
finalement, " #
Z b n−1
X
h
f (x)dx ≈ f (x0 ) + 2 f (xi ) + f (xn )
a 2 i=1

où bien " #
Z b n−1
X
h b − a 2 00
f (x)dx = f (a) + 2 f (xi ) + f (b) − h f (ξ),
a 2 i=1
12
où ξ ∈]a, b[. Le degré d’exactitude est à nouveau égal à 1.
La méthode composée est la méthode des trapèzes définie par :
n−1
!
b−a f (a) + f (b) X
Tn (f ) = + f (xi )
n 2
k=1

Theorem 4.4.2 Soit f une fonction de classe C 2 sur [a, b] et n un entier naturel non nul. Avec les notations
qui précèdent, il existe un réel x i ∈ [a, b] tel que :
Z b
(b − a)3 00
f (x)dx − Tn (f ) = − f (ξ)
a 12n2
et on a : Z b
(b − a)3 00
f (x)dx − Tn (f ) ≤ kf k∞
a 12n2

4.4.3 Formule de Cavalieri-Simpson

Soit [a, b] un intervalle de R. On définit c = (b + a)/2 et soit p1 (x) un polynôme de degré ≤ 2 qui
interpole f aux points a, c et b. Alors,

p2 (x) = A + Bx + Cx2 ,

où A, B et C sont des constantes à déterminer tels que

p2 (a) = f (a), p2 (c) = f (c) p2 (b) = f (b)

Après avoir déterminé A, B et C, on obtient

Z b
b−a
p2 (x)dx = [f (a) + 4f (c) + f (b)]
a 6
et le rôle de Simpson est
Z b Z b
b−a
f (x)dx ≈ p2 (x)dx = [f (a) + 4f (c) + f (b)] .
a a 6
52 CHAPITRE 4. INTÉGRATION ET DÉRIVATION NUMÉRIQUE

On peut montrer que, si f ∈ C 4 ([a, b]), l’erreur de quadrature est

h5 (4) (b − a)5 (4) b−a
E2 (f ) = − f (ξ) = − f (ξ), h = ,
90 2880 2
où ξ ∈]a, b[. On en déduit que le degré d’exactitude est 3.
Z b
b−a (b − a)5 (4)
f (x)dx = [f (a) + 4f (c) + f (b)] − f (ξ), ξ ∈]a, b[.
a 6 2880
*Formule composite d’intégration par la méthode de Simpson :
Pour n ≥ 0, soit h = (b − a)/n le pas d’une subdivision S de l’intervalle [a, b]. On définit xi = a + ih pour
i = 0, . . . , n où x0 = a et xn = b.
Z b Xn Z xi
f (x)dx = f (x)dx
a i=1 xi−1
n 
X 
h
≈ f (xi−1 ) + 4f (xi− 1 ) + f (xi )
6 i=1
2

où xi− 1 = 12 (xi + xi−1 ).

2
Si f ∈ C 4 ([a, b]) alors l’erreur de quadrature associée à cette intégration numérique est donnée par
b − a 4 (4)
E(f ) = − h f (ξ)
180 × 24
où ξ ∈]a, b[. Et, dans ce cas on écrit
Z b
h X  b−a
n
f (x)dx = f (xi−1 ) + 4f (xi− 1 ) + f (xi ) − h4 f (4) (ξ).
a 6 i=1 2 2880
*Autre formule composite de Simpson :
Si on introduit les nœuds de quadrature xk = a + k2 h pour k = 0, . . . , 2n, avec h = (b − a)/n le pas d’une
subdivision S de l’intervalle [a, b]. on a alors
Z b X n Z xi
f (x)dx = f (x)dx
a i=1 xi−1
n−1 n−1
!
h X X
≈ f (a) + 2 f (x2i ) + 4 f (x2i+1 ) + f (b)
6 i=1 i=0
où x0 = a et x2n = b.
La méthode de Simpson est définie par :
h X 
n
Sn (f ) = f (xi−1 ) + 4f (xi− 1 ) + f (xi )
6 i=1 2

n−1 n−1
!
h X X
= f (a) + 2 f (x2i ) + 4 f (x2i+1 ) + f (b) .
6 i=1 i=0
4
Theorem 4.4.3 Soit f une fonction de classe C sur [a, b] et n un entier naturel non nul. Avec les notations
qui précèdent, il existe un réel ξ ∈ [a, b] tel que :
Z b
(b − a)5 (4)
f (x)dx − Sn (f ) = − f (ξ)
a 2880n4
et on a : Z b
(b − a)5 (4)
f (x)dx − Sn (f ) ≤ f ∞
a 2880n4
Chapitre 5

Méthode des moindres carrés et optimisation

quadratique

Nous nous intéresserons dans cette partie au problème de maximisation ou minimisation des fonctions
de plusieurs variables. Commençons par les fonctions de deux variables. Nous rappelons ici un théorème
fondamental, qui va servir dans ce chapitre

Théorème 5.0.2 Une matrice symétrique A est définie positive si et seulement si det(Ak ) > 0 pour toutes
les sous-matrices principales Ak de A.

5.1 Maxima et minima de fonctions de deux variables

5.1.1 Gradient d’une application et Matrice hessienne d’une F.P.V
Soit f : Rn −→ Rp une application définie par
f (X) = (f1 (X), f2 (X), . . . , fp (X))
où X = (x1 , . . . , xn ) ∈ Rn . Supposons que l’application est différentiable sur Rn−
. −→
– Pour p = 1. On appelle Gradient de f au point X0 , noté ∇f (X0 ) = gradf (X0 ), est le vecteur
pointant dans la direction de croissance maximale de la fonction f où bien il est le vecteur définie par
 ∂f 
∂x1
(X0 )
 
 ∂f 
 
∇f (X0 ) =  ∂x2 (X0 ) 
 .. 
 . 
∂f
∂xn
(X0 )
– Pour p > 1. On appelle Gradient de f au point X0 , noté ∇f (X0 ), la matrice jacobienne de taille
(n × p) définie par
 ∂f1 ∂f1 ∂f1 
∂x1
(X0 ) ∂x 2
(X0 ) . . . ∂x n
(X0 )
 
 ∂f 
 (X0 )
2 ∂f 2
(X0 ) . . . ∂xn (X0 )
∂f 2
 ∂x1 ∂x2 
Jf (X0 ) =  . . . . 
 .. .. .. .. 
 
 
∂fp ∂fp ∂fp
∂x1
(X0 ) ∂x2 (X0 ) . . . ∂xn (X0 )
53
54 CHAPITRE 5. MÉTHODE DES MOINDRES CARRÉS ET OPTIMISATION QUADRATIQUE

Définition 5.1.1 On appelle matrice hessienne de f au point X0 , noté Hf (X0 ), la matrice de taille (n×n)
définie par  ∂2f 
∂2f ∂2f
2
∂x1
(X 0 ) ∂x1 ∂x2
(X 0 ) . . . ∂x1 ∂xn
(X 0 )
 
 
 ∂2f 2
∂ f ∂f2 
 ∂x2 ∂x1 (X0 ) 2 (X0 ) . . . ∂x2 ∂xn (X0 )
Hf (X0 ) = 
 ..
∂x 2
.. .. ..


 . . . . 
 
 
∂2f ∂2f ∂2f
∂xn ∂x1
(X0 ) ∂xn ∂x2 (X0 ) . . . ∂x2
(X0 )
n

5.1.2 Approximations linéaire et quadratique : Formule de Taylor

On considère f : Rn −→ R une fonction de plusieurs variables définie par

f (X) = f (x1 , . . . , xn )

Supposons que f est de classe C 2 au voisinage d’un point X0 . Soit h ∈ Rn .

La formule de Taylor de f en le point x0 (ou la formule d’approximation) est donnée par l’expression
suivante :
1
f (X0 + h) = f (X0 ) + (∇f (X0 ), h) + (Hf (X0 )h, h) + o(khk2 ).
2
Que se passe-t-il si ∇f (X0 ) = 0Rn ? si (Hf (X0 )h, h) < 0 ? si (Hf (X0 )h, h) > 0 ? si (Hf (X0 )h, h) = 0 ?

5.1.3 Points critiques d’une application

Considérons une fonction de deux variables définie par f (x, y). Nous supposerons que f est de classe
C 3 (i.e. f et ses dérivées partielles d’ordre≤ 3 sont continues).

Définition 5.1.2 On dit que (x0 , y0 ) est un point critique de f si fx (x0 , y0) = 0 et fy (x0 , y0 ) = 0 où fx et
fy dénotent les dérivées partielles de f par rapport à x et y respectivement.

Soit X0 = (x0 , y0 ) un point critique de f . En utilisant le formule de Taylor on obtient :

1 2 
f (x0 + h, y0 + k) = f (X0 ) + hfx (X0 ) + kfy (X0 ) + h fxx (X0 ) + 2hkfxy (X0 ) + k 2 fyy (X0 ) + R
2
où R = 61 [h3 fxxx (X) + 3h2 kfxxy (X) + 3hk 2 fxyy (X) + k 3 fyyy (X)] avec X = (x0 + θh, y0 + θk) (0 <
θ < 1).
Or fx (x0 , y0) = 0 et fy (x0 , y0 ) = 0 puisque X0 est un point critique. Il en résulte que
1 2 
f (x0 + h, y0 + k) = f (X0 ) + ah + 2bhk + ck 2 + R
2
où on a posé a = fxx (X0 ), b = fxy (X0 ) et c = fyy (X0 ). Si h et k sont suffisamment petits, on montre alors
que
1
|R| ≤ |ah2 + 2bhk + ck 2 |
2
Par conséquent, f (x0 + h, y0 + k) − f (x0 , y0 ) a le même signe que ah2 + 2bhk + ck 2 pour h et k petits, ce
qui entraîne que (x0 , y0 ) est minimum local si

ah2 + 2bhk + ck 2 > 0, pour tout(h, k),

5.1. MAXIMA ET MINIMA DE FONCTIONS DE DEUX VARIABLES 55

est un maximum local si

ah2 + 2bhk + ck 2 < 0, pour tout(h, k),
ou est un point selle si ah2 + 2bhk + ck 2 prend des valeurs positives et des valeurs négatives.
En d’autres mots, on a minimum local si la forme quadratique

q(h, k) = ah2 + 2bhk + ck 2

est définie positive,un maximum local si elle est définie négative ou un point selle si elle est indéfinie. La
forme quadratique s’écrit   
fxx (X0 ) fxy (X0 ) h
q(h, k) = (h, k)
fyx (X0 ) fyy (X0 ) k
La matrice précédente  
fxx (X0 ) fxy (X0 )
H(X0 ) =
fyx (X0 ) fyy (X0 )
définissant la forme quadratique est appelée la matrice hessienne de f en X0 . Soient λ1 et λ2 les valeurs
propres de H(X0). On déduit du théorème 5.0.2 :
1. (x0 , y0) est un minimum local si λ1 > 0 et λ2 > 0.
2. (x0 , y0) est un maximum local si λ1 < 0 et λ2 < 0.
3. (x0 , y0) est un point selle si λ1 λ2 < 0.
Comme on l’a vu au théorème 5.0.2, la matrice hessienne H(X0 ) est définie positive si ses déterminants
principaux fxx (X0 ) et ∆ = det(H(X0 )) = fxx (X0 )fyy (X0 ) − (fxy (X0 )))2 sont > 0.
La matrice hessienne est définie négative si −H(X0 ) est définie positive, ce qui sera le cas si ses détermi-
nants principaux −fxx (X0 ) et det(−H(X0 )) = det(H(X0 )) = ∆ sont > 0, c’est-à-dire fxx (X0 ) < 0 et
∆ > 0. Il résulte de ce qui précède que H(X0 ) est indéfinie si ∆ < 0. Ainsi, on obtient :
1. Si ∆ > 0 et fxx (X0 ) > 0, alors X0 = (x0 , y0) est un minimum local.
2. Si ∆ > 0 et fxx (X0 ) < 0, alors X0 = (x0 , y0) est un maximum local.
3. Si ∆ < 0, alors X0 = (x0 , y0 ) est un point selle.

Remarque 5.1.1 Si ∆ = 0, on peut rien conclure quant à la nature du point critique.

Exemple 5.1.1 Considérons la fonction donnée par

1
f (x, y) = x3 + xy 2 − 4xy + 1.
3
Déterminons d’abord les points critiques.
Ce sont les solutions du système suivant :
 ∂F
fx (x, y) = ∂x
(x, y) =0
∂F
fy (x, y) = ∂y
(x, y) =0

On a  ∂F
fx (x, y) = ∂x
(x, y) = x2 + y 2 − 4y
∂F
fy (x, y) = ∂y
(x, y) = 2xy − 4x
On obtient les quatre points critiques : X1 = (0, 0), X2 = (0, 4), X3 = (2, 2) et X4 = (−2, 2). La matrice
hessienne de f est  
∂2F 2x 2y − 4
H(x, y) = (x, y) = .
∂x∂y 2y − 4 2x
56 CHAPITRE 5. MÉTHODE DES MOINDRES CARRÉS ET OPTIMISATION QUADRATIQUE

En évaluant la matrice hessienne aux points X1 , . . . , X4 on a

   
0 −4 0 4
H(0, 0) = , H(0, 4) = ,
−4 0 4 0
   
−4 0 4 0
H(2, 2) = , H(−2, 2) = .
0 −4 0 4
Les valeurs propres de ces matrices sont :
1. Pour H(0, 0) sont λ1 = 4 et λ2 = −4.
2. Pour H(0, 4) sont λ1 = 4 et λ2 = −4.
3. Pour H(2, 2) sont λ1 = 4 et λ2 = 4.
4. Pour H(−2, 2) sont λ1 = −4 et λ2 = −4.
Par conséquent, X1 et X2 sont des points selles, X3 est minimum local et X4 est un maximum local.

5.1.4 Maxima et minima des fonctions de n variables

Considérons une fonctions de n variables définie par F (x) = F (x1 , . . . , xn ). Un point a = (a1 , . . . , an )
∂F
est un point critique si Fxi (a) = 0 pour i = 1, . . . , n, où Fxi = ∂x i
désigne la dérivée partielle par rapport à
xi .
C’est-à-dire que Les points critiques de l’application F : (x1 , . . . , xn ) 7−→ F (x1 , . . . , xn ) sont les solutions
de l’équation vectorielle suivante :
∇F (x1 , . . . , xn ) = 0
ce qui équivalent au système d’équations suivant
 ∂F

 ∂x1
(x1 , . . . , xn ) = 0





 ∂F
∂x2
(x1 , . . . , xn ) = 0



 ..

 .

 ∂F
∂xn
(x1 , . . . , xn ) = 0

Utilisant, comme pour les fonctions de deux variables, la formule de Taylor on montre qu’un point critique
X
1. est un minimum local si la matrice hessienne H(X) est définie positive,
2. est un maximum local si H(X) est définie négative
3. est un point selle si H(X) est indéfinie.
Ici, la matrice hessienne est la matrice symétrique n × n définie par
 
∂2F
H(X) = (Fxi xj (X)) = (x1 , . . . , xn ) ,
∂xi ∂xj
où i = 1, . . . , n et j = 1, . . . , n.
5.2. FONCTIONS QUADRATIQUES 57

Exemple 5.1.2 On considère la fonction définie par

F (x, y, z) = x2 + xz − 3 cos(y) + z 2 .

Déterminons les points critiques. Les dérivées partielles sont

 ∂F
 ∂x (x, y, z) = 2x + z
∂F
∂y
(x, y, z) = 3 sin(y)
 ∂F
∂z
(x, y, z) = x + 2z

Ainsi, les points critiques sont les solutions de


 2x + z = 0
3 sin(y) = 0

x + 2z = 0

Il en résulte que les points critiques sont les points Xk = (0, kπ, 0), k ∈ Z. La matrice hessienne au point
Xk est  
2 0 1
H(Xk ) = 0 3(−1)k 0 .
1 0 2
– Si k est pair, les valeurs propres de H(Xk ) sont 3, 3 et 1. Ainsi, H(Xk ) est définie positive si k est
pair.D’où il s’ensuit que Xk est un minimum local si k pair.
– si k est impair, les valeurs propres sont 3, −3 et 1. Ainsi, H(Xk ) est indéfinie si k est impair. D’où
Xk est un point selle si k est impair.

5.2 Fonctions quadratiques

5.2.1 Forme linéaires et bilinéaires

Définition 5.2.1 Une forme linéaire ` sur un espace vectoriel réel E est une application linéaire de E dans
R.
C’est-à-dire pour tout x, y ∈ E et α, β ∈ R on a

f (αx + βy) = αf (x) + βf (y).

Si E est dimension finie n, (e1 , . . . , en ) sa base canonique et L = (`(e1 ), . . . , `(en )) le vecteur ligne des
images de ` dans la base (e1 , . . . , en ). Soit x = x1 e1 + . . . + xn en un vecteur dans E et X = (x1 , . . . , xn ),
alors
`(x) = (L, X)
En dimension finie, toute forme linéaire est représentée par un produit scalaire.

Définition 5.2.2 Une forme bilinéaire a sur un espace vectoriel réel E est une application de E × E dans
R, linéaire par rapport à chacune de ses deux arguments.
Soit a la forme bilinéaire, on a donc :

a(u, αv1 + βv2 ) = αa(u, v1 ) + βa(u, v2) (1.1)

a(αu1 + βu2 , v) = αa(u1 , v) + βa(u2 , v) (1.2)
58 CHAPITRE 5. MÉTHODE DES MOINDRES CARRÉS ET OPTIMISATION QUADRATIQUE

Représentation matricielle

Toute forme bilinéaire sur un espace E de dimension finie n se représente, dans une basee (e1 , . . . , en )
par une matrice carrée d’ordre n. Les coefficients Aij de la matrice A représentant l’application a sont
donnés par
Aij = a(ei , ej )
Soit u = u1 e1 + . . . + un en et v = v1 e1 + . . . + vn en on a

a(u, v) = (AU, V ) = (U, AT V )

où U = (u1 , . . . , un )T et V = (v1 , . . . , vn )T .
Si (·, ·) représente le produit scaalaire usuel de Rn et AT est la matrice transposée de A définie par

ATij = Aji , ∀i, j = 1, . . . , n

Définition 5.2.3 1. Une forme bilinéaire a sur un espace vectoriel réel E est symétrique si :

∀u, v ∈ E, a(u, v) = a(v, u).

2. Une forme bilinéaire a sur un espace vectoriel réel E est définie positive si :

∀u ∈ E, a(u, u) ≥ 0 et a(u, u) = 0 ⇔ u = 0

Les formes bilinéaires symétriques sont représentées par des matrices symétriques Aij = Aji .
Les formes bilinéaires symétriques définies positives sont représentées par des matrices symétriques définies
positives, qui vérifient donc

(AU, U) ≥ 0, ∀U ∈ Rn et (AU, U) = 0 ⇒ U = 0

Théorème 5.2.1 (Résultat important)

1. Les matrices symétriques réelles ont des valeurs propres réelles, sont diagonalisables et admettent
une base de vecteurs propres orthonormés.
2. Les matrices symétriques réelles définies positives ont des valeurs propres strictement positives, et
donc sont inversibles.

5.2.2 Équivalence entre la résolution d’un système linéaire et la minimisation qua-

dratique

On considère a : Rn × Rn → R une forme bilinéaire et A la matrice associée à a relativement à la base

(e1 , . . . , en ) :
Aij = a(ei , ej ), i, j = 1, . . . , n
On définit la fonctionnelle J : Rn → R donnée pour tout v ∈ Rn par
1
J(v) = (Av, v) − (b, v)
2
où b ∈ Rn est un vecteur donné.
5.3. APPLICATION AUX MOINDRES CARRÉS 59

Théorème 5.2.2 Si A est une matrice symétrique définie positive, il y a équivalence entre les trois pro-
blèmes suivants : 
trouver X ∈ Rn tel que
(1)
AX = b

trouver X ∈ Rn tel que
(2)
(AX, Y ) = (b, Y ), ∀Y ∈ Rn

trouver X ∈ Rn tel que
(3)
J(X) = 21 (AX, X) − (b, X) soit minimale.

Démonstration. (1) ⇒ (2) est évident.

(2) ⇒ (1) en prenant pour Y les vecteurs de base ei de Rn .
(2) ⇒ (3) On calcule J(X + λY ) pour tout λ réel et tout Y ∈ Rn , on obtient

1
J(X + λY ) = J(X) + λ ((AX, Y ) − (b, Y )) + λ2 (AY, Y )
2
en utilisant la symétrie de la matrice A.
On en déduit, si (AX, Y ) − (b, Y ) = 0 que J(X + λY ) = J(X) + 21 λ2 (AY, Y )
d’où en utilisant le fait que A est définie positive :

J(X + λY ) > J(X)

si λ et Y sont non nuls. Donc on a montré que si X vérifie (2), X minimise J.

(3) ⇒ (2), car X minimise J, on a

1
λ ((AX, Y ) − (b, Y )) + λ2 (AY, Y ) ≥ 0, ∀λ, ∀Y
2
le trinôme en λ ci-dessus doit être toujours positif. Ceci entraîne que son discriminant soit toujours négatif
ou nul. Or ce discriminant est
4 = ((AX, Y ) − (b, Y ))2
ceci implique (2). 

5.3 Application aux moindres carrés

Considérons le problème générale d’un système linéaire sur-déterminé, c’est-à-dire dans lequel il y a
plus d’équations que d’inconnues. C’est en particulier le cas dans le calcul de la droite des moindre carré
ou plus généralement de polynômes d’approximation au sens des moindres carrés. On ne peut pas obtenir
exactement l’égalité
AX = b
où A est une matrice rectangulaire de n lignes et m colonnes avec n > m. On définit

J(X) = kAX − bk2 = (AX − b, AX − b).

On essaie alors de minimiser l’écart entre les vecteurs AX et B de

min J(X)
X∈Rm
60 CHAPITRE 5. MÉTHODE DES MOINDRES CARRÉS ET OPTIMISATION QUADRATIQUE

en minimisant la norme euclidienne de leur différence, ou ce qui revient au même le carré de cette norme.
On utilise les propriétés classiques du produit scalaire (AU, V ) = (U, AT V ) pour obtenir :

J(X) = (AT AX, X) − 2(AT b, X) + (b, b)

la matrice AT A est une matrice carrée (m × m) symétrique définie positivedès lors que la matrice rectan-
gulaire A est de rang m. Le théorème (5.2.2) nous donne l’équivalence de ce problème de moindre carrés
avec la résolution du système linééaire
AT AX = AT b
On retrouve ainsi le système carré de m équations à m inconnues, dit “système des équations normales”.

5.3.1 Approximation par la droite des moindre carrés

Par exemple, dans le cas de la droite des moindre carrés, il s’agit de trouver la fonction affinie

y = α + βx

qui représente “au mieux” une collection de n valeurs yi associées aux n abscisses xi . Au sens des moindres
carrés , ceci revient à minimiser la somme
n
X
(yi − (α + βxi ))2
i=1

donc à minimiser la norme euclidienne de la différence

J(a) = kAa − bk2

où l’on a noté b le vecteur des n valeurs yi , A la matrice rectangulaire à n lignes et 2 colonnes

 
1 x1
1 x2 
 
 .. .
. 
A = . . 
 
1 xn−1 
1 xn
 
α
et a = . En appliquant les résultats précédents, on obtient la solution en résolvant le système (2 × 2)
β
suivant
AT Aa = AT b

5.3.2 Interprétation géométrique : projection sur un sous-espace

Reprenons le problème de l’approximation au sens des moindre carrés

J(x) = kAx − bk2

On peut interpréter ce problème comme celui de la recherche du vecteur de la forme Ax le plus proche au
sens de la norme euclidienne d’un vecteur b donné dans Rn .
Les vecteurs de la forme Ax sont une combinaison linéaire des m vecteurs colonnes de la matrice A. Ces
5.3. APPLICATION AUX MOINDRES CARRÉS 61

vecteurs colonnes sont indépendants car A est supposée de rang m.

Ils engendrent donc un sous-espace vectoriel F de Rn de dimension m. Et le problème s’interprète comme
la recherche du vecteur du sous-espace F le plus proche du vecteur b.
On obtient donc x en écrivant que Ax est la projection orthogonale de b dans F donc

(Ax − b, V ) = 0, ∀V ∈ F

ceci car Rn = F ⊕ F ⊥ et que Ax − b ∈ F ⊥ .

ce qui est équivalent, puisque les colonnes de A engendrent F , à

(Ax − b, Aj ) = 0, ∀j = 1, . . . , n

où Aj est le j ime vecteur colonne de A. On retrouve ainsi le résultat

AT Ax = AT b.