ANALYSE NUMERIQUE MASTER GEOMATIQUE

May 15, 2025

Contents

1 Rappels 3

1.1 Formule de Taylor . . . . . . . . . . . . . . . . . . . 3

1.1.1 Comment evaluer le polynôme a0 + a1x +
a2x2 + · · · + anxn . . . . . . . . . . . . . . . 5

2 Interpolation polynomiale 6

2.1 Interpolation polynomiale de Lagrange . . . . . . . . 6

2.1.1 Méthode de Newton . . . . . . . . . . . . . . 8
2.1.2 Méthode progressive de Newton . . . . . . . . 11
2.1.3 Méthode regressive de Newton . . . . . . . . 11
2.2 Interpolation inverse de Lagrange . . . . . . . . . . . 13
2.3 Interpolation de Hermite . . . . . . . . . . . . . . . 15

3 Résolution des systèmes linéaires 19

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Méthodes directes de résolution de systèmes linéaires
Au = b . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2.1 Systèmes triangulaires . . . . . . . . . . . . 20
3.2.2 Elimination naïve de Gauss . . . . . . . . . . 23
3.2.3 Triangularisation . . . . . . . . . . . . . . . . 23
3.2.4 Factorisation LU . . . . . . . . . . . . . . . 28
3.2.5 Elimination de Gauss avec stratégie de pivot 31
3.2.6 Système linéaires spéciaux : factorisations
de Cholesky, LDM t ... . . . . . . . . . . . . . 37

1
3.3 Méthodes itératives de résolution de systèmes linéaires
Au = b . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.1 Méthodes itératives standard : Jacobi, Gauss
Seidel, méthodes de relaxation . . . . . . . . 41

2
Chapter 1
Rappels

La formule de Taylor constitue un outil fondamental d'analyse tout

au long de ce cours. Sa formulation universelle est celle dite avec
reste intégral. Cette forme avec reste intégral se prête bien aux
généralisations aux cas des fonctions vectorielles et des fonctions à
plusieurs variables. Elle peut s'énoncer :

1.1 Formule de Taylor

Théorème 1.1.1 Soit I un intervalle ouvert non vide de RI . Soit

I ∗. Soit f ∈ C m(I) et x ∈ I , alors pour tout y ∈ I on a
m ∈N
m−1 y
(y − x)k (k) (y − t)m−1 (m)
X Z
f (y) = f (x) + f (t)dt
k! x (m − 1)!
k=0

Remarque :
En posant h = y − x, on a
m−1 x+h
hk (k) (x + h − t)m−1 (m)
X Z
f (x + h) = f (x) + f (t)dt
k! x (m − 1)!
k=0

m−1 1
hk (k) hm
X Z
f (x + h) = f (x) + (1 − λ)m−1f (m)(x + λh)dλ
k! (m − 1)! 0
k=0

3
Remarque :
f ∈ C m(I) s'entend une fonction m fois continûment dérivable à
valeur dans R
In
Exemple 1 f (t) = eit = cos t + i sin t t ∈ [0, 2π]. Si la formule
avec reste de Lagrange est vraie , ∃c ∈]0; 2π[ tel que :
f (2π) − f (0) = f 0(c)(2π − 0)
Or f 0(t) = ieit et |f 0(t)| = 1 ∀ t contradiction.
Exemple 2 f (x) = e ∀i, f (x) = e et on a
x i x

n
X f (i)(x0) f (n+1)(ξ)
i
f (x) = f (x0) + (x − x0) + (x − x0)(n+1)
i=1
i! (n + 1)!
n
x0
X ex0 eξ
i
=e + (x − x0) + (x − x0)(n+1)
i=1
i! (n + 1)!
n
x0
X (x − x0)i eξ
= e [1 + ]+ (x − x0)(n+1)
i=1
i! (n + 1)!
Posons h = x − x0 , on a:
n
x0
X hi eξ
f (x0 + h) = e [1 + ]+ h(n+1)
i=1
i! (n + 1)!
Pour x0 = 0 , On a
n
X hi eξ
f (h) = [1 + ]+ h(n+1)] 0 < ξ < h
i=1
i! (n + 1)!
eξ (n+1) eξ
|Rn| = | h |≤ |h|(n+1)
(n + 1)! (n + 1)!
Si on veut avoir |Rn| <  on a (n+1)! |h|(n+1) ≤  , 0 < h < 1 ,
e ξ

0 <  < 1 1 < eξ < e < 3.

Il sut que (n + 1)! ≥ 3 pour  = 10−6 on a (n + 1)! ≥ 310−6 on
prend n = 9
4
1.1.1 Comment evaluer le polynôme a0 + a1 x + a2 x2 + · · · + an xn

Méthode de Horner

P (x) = (((anx + an−1)x + an−2)x + · · · )x + a0

Algorithmique

P ←− an
P ←− P x + an−1
...
P ←− P x + a1
P ←− P x + a0
Exemple
n
h
X hi
e =1+ + Rn
i=1
i!
Pour n connu on a:
P ←− nh
h
P ←− (P + 1) n−1
...
P ←− (P + 1) i!h
P ←− (P + 1) h1
P ←− (P + 1)
Exercice : Donner une approximation de ln 2 en utilisant le développe-
ment de Taylor à l'ordre 5 de f (x) = ln(x + 1) au voisinage de 0
puis majoré l'erreur commise.

5
Chapter 2
Interpolation polynomiale

2.1 Interpolation polynomiale de Lagrange

Théorème

Soit (xi, yi)ni=0 , n + 1 couples de réels tels que xi 6= xj si i 6= j alors

il existe un polynôme Pn et un seul de dégre ≤ n tel que Pn(xi) = yi
et on a n
X
Pn(x) = yiLi(x)
i=0
avec n
Y x − xj
Li(x) = ( )
xi − xj
j=1,j6=i

(Li)ni=0 est la base de Lagrange et Pn est dit polynôme d'interpolation

de Lagrange relatif aux donées (xi, yi)ni=0.
Si f est la fonction qui x 7−→ y = f (x), on dit que Pn est le
polynôme d'interpolation de la fonction f relativement aux n + 1
n÷uds distints (x0, · · · , xn)
Remarque p ∈ Pn ⇐⇒ ∃(a0 , a1 , · · · , an ) ∈ R
n+1
tel que :
n
X
p(x) = aixi
i=0

Le problème d'interpolation de Lagrange se ramène à chercher

6
(a0, a1, · · · , an) ∈ Rn+1 tel que
n
X
aixik = yk ∀k ∈ {0, 1, · · · , n}
i=0

D'après le théorème, le système

n
X
aixik = yk ∀k ∈ {0, 1, · · · , n}
i=0

admet une solution unique trouvée à travers le polynôme d'interpolation

n
X
P = y i Li .
i=0

Comme la dimension de Pn est n + 1 alors (Li)ni=0 constitue un

système générateur minimal. C'est-à-dire une base de Pn
Exemple 1 On considère le polynôme f déni par la table suivante
:
x 1.4 1.25
y 3.4 3.9
Evaluer f (1.3) en utilisant le polynôme d'interpolation de La-
grange.
Solution : n + 1 = 2 =⇒ n = 1
x − 1.25 x − 1.25
L1(x) = =
1.4 − 1.25 0.15

x − 1.4 x − 1.4
L2(x) = =−
1.25 − 1.4 0.15

P (x) = 3.4L1(x) + 3.9L2(x)

7
d'où
P (x) = − 43 x + 16.7
3

et on a

f (1.3) = p(1.3) = − 43 (1.3) + 16.7

3

Exemple 2 On considère le polynôme f déni par la table suivante

:
x 1 2 3
y 2 -1 7
Détermine le polynôme d'interpolation P de f
Solution : n + 1 = 3 =⇒ n = 2
(x − 2)(x − 3)
L1(x) =
(1 − 2)(1 − 3)

(x − 1)(x − 3)
L2(x) =
(2 − 1)(2 − 3)

(x − 1)(x − 2)
L3(x) =
(3 − 1)(3 − 2)

Donc P = 2L1 − L2 + 7L3 par suite P (x) = 112 x2 − 392 x + 16

2.1.1 Méthode de Newton

Soit f une fonction réelle dénie sur [a, b] (a ≤ b) x0, x1, x2, · · · , xn
n+1 points distints de [a, b]. On cherche pk pour k ∈ N I ∗ (k ≤ n).
pk ∈ Pk polynôme d'interpolation de f aux points x0, x1, x2, · · · , xk .
On veut construire le pn par recurrence sur k .
Pk+1(xi) − Pk (xi)

8
C'est-à-dire
k
Y
Pk+1(x) = Pk + Ck+1 (x − xi)
i=0
avec
Ck+1 = f [x0, x1, x2, · · · , xk , xk+1]
diérence divisée d'ordre k + 1 aux points x0, x1, x2, · · · , xk , xk+1.
Pn interpolant f aux points x0, x1, x2, · · · , xn sera donné par :
Pn(x) = f [x0]+f [x0, x1](x−x0)+f [x0, x1, x2](x−x0)(x−x1)(x−x2)+
k−1
Y
+ · · · f [x0, x1, x2, · · · , xk ] (x − xi)
i=0
n−1
Y
+f [x0, x1, x2, · · · , xn] (x − xi).
i=0
Donc
n
X k−1
Y
f (x) = f [x0] + f [x0, x1, x2, · · · , xk ] (x − xi)
k=1 i=0

f [x0, x1, x2, · · · , xk ] est le coecient de xk dans Pk . f [x0] = f (x0)

plus généralement f [xi] = f (xi) et p1(x) = f [x0] + f [x0, x1](x − x0)
p∗1 (x) = f [x1] + f [x1, x0](x − x1) de fait de l'unicité du polynôme
d'interpolation on a f [x1, x0] = f [x0, x1] et
f [x1] − f [x0]
f [x0, x1] =
x1 − x0
Théorème

f [xi1, xi2, xi3, · · · , xik ] − f [xi0, xi1, xi2, · · · , xik−1]

f [xi0, xi1, xi2, · · · , xik ] =
xik − xi0
Démonstration 1 Soit xi0, xi1, xi2, · · · , xin distincts. Pk−1 ∈ Pk−1
polynôme de Lagrange de f aux points xi0, xi1, xi2, · · · , xik−1 Pk−1
∗
∈
9
Pk−1 polynôme de Lagrange de f aux points xi0, xi1, xi2, · · · , xik−1,xik .
D'après l'algorithme de Neville on a :
∗
(x − xi0)Pk−1 (x) − (x − xik Pk−1(x))
Pk (x) =
xik − xi0
f [xi0, xi1, xi2, · · · , xik−1] est le c÷cient de xk−1 dans Pk−1.
f [xi0, xi1, xi2, · · · , xik ] est le c÷cient de xk−1 dans Pk−1
∗
.
f [xi0, xi1, xi2, · · · , xik ] est le c÷cient de xk dans Pk .
Remarque : Table de diérence divisée

Ordre 0 Ordre 1 Ordre 2 Ordre 3

xi f [xi] f [xi, xi+1] f [xi, xi+1, xi+2] f [xi, xi+1, xi+2, xi+3]
x0 f (x0)

f [x0, x1]

x1 f (x1) f [x0, x1, x2]

f [x1, x2] f [.....]

x2 f (x2) f [x1, x2, x3]

f [x2, x3] f [.....]

x3 f (x3) f [x2, x3, x4]

f [x3, x4]

x4 f (x4)
x 1.25 1.4
Exemple 3 Soit f la fonction dénie par :
y 3.9 3.7

10
Donner une approximation de f (1.3) en utilisant la méthode de
Newton progressive.
1.3 − xi xi f [xi] f [xi, xi+1]
0.05 1.25 3.9

3.7−3.9
1.4−1.25 = − 34

−0.1 1.4 3.7

Ainsi d'après la formule de Neuton :
4
P (x) = 3.9 − (x − 1.25)
3
4 16.7
P (x) = − x +
3 3
Donc P (1.3) = 3.9 − 43 (0.05)

2.1.2 Méthode progressive de Newton

P (x) = f [x0] + f [x0, x1](x − x0) + f [x0, x1, x2](x − x0)(x − x1) + · · ·
n−1
Y
+f [x0, x1, x2, · · · , xn] (x − xi).
i=0

2.1.3 Méthode regressive de Newton

P (x) = f [xn]+f [xn−1, xn](x−xn)+f [xn, xn−1, xn−2](x−xn)(x−xn−1)+

n−1
Y
+f [x0, x1, x2, · · · , xn] (x − xi).
i=0

x 1 2 3
Exemple 4 On considère la fonction f donnée par
y 2 -1 7
Estimer f (1.5) et f (2.5) par la méthode de Newton.
11
 1.5 − xi 2.5 − xi xi f [xi] f [xi, xi+1] f [, , ]
0.5 1.5 1 2

−3

On a 11
2
−0.5 0.5 2 −1

−1.5 −0.5 3 7
Ainsi
P (1.5)=2 − 3 × 0.5 + 11
2 × (0.5)(−0.5)

= 2 − 1.5 − 112 × 0.25

P (1.5)= −7
8
P (2.5)=7 + 8 × (−0.5) + 11 2 × (0.5)(−0.5)

= 7 − 4 − 112 × 0.25

P (2.5)= 13
8

x 1 -1 2
Exercice 1 Soit la fonction f donnée par
y 0 -3 4
Donner une expression du polynôme de Lagrange de f
1. En construisant la base de Lagrange
2. Par l'algorithme de Neville

12
 3. Par la méthode progressive de Newton
4. Par la méthode regressive de Newton
Exercice 2 Soit la fonction f donnée par
x 1 1.3 1.6 1.9 2.2
y 0.7652 0.6201 0.4564 0.2818 0.1104

1. Evaluer f (1.5) en utilisant le polynôme de Lagrange à travers

(a) L'algorithme de Neville
(b) La méthode progressive de Newton
(c) La méthode regressive de Newton
(d) La base de Lagrange
2. Dans chaque cas précicer le nombre d'opérations ottant.
Remarque Soit n ∈ N
I ∗ f ∈ Pn k ∈ N I tel que k > n
x0 < x1 < x2 · · · xk ∈ [a, b]. Pk polynôme de Lagrange inter-
polant f aux points x0 < x1 < x2 · · · xk

2.2 Interpolation inverse de Lagrange

Soit f dénie sur [a; b]. Soit a = x0 < x1 < · · · < xn = b, (n + 1)

points distincts de [a; b]. On suppose que f est injective sur [a; b] et
continue. Les yi = f (xi), i = 0, . . . , n sont distincts.
On appelle polynôme d'interpolation inverse de lagrange de f aux
points x0, x1, · · · , xn, le polynôme d'interpolation de la fonction
f −1, réciproque de la fonction f aux points y0, . . . , yn.

y ∈ f ([a; b]) 7−→ x = f −1(y) ∈ [a; b].

13
 x 1 −1 2
Exemple 3 Soit la fonction f dénie par :
y 0 −3 4
Calculons f −1(−1) et f −1(2) On a :
x −1 1 2
y −3 0 4
Nous obtenons la table suivante :
y x f −1[ ; ] f −1[ ; ; ] −1 − y
−3 −1 2

2
3
1 2
−
0 1 4 3 = 5 −1
4 − (−3) 84

1
4

4 2 −5
2 5 19
f −1(−1) = −1 + (2) − (2)(−1) =
3 84 42
2 5 73
f −1(2) = −1 + (5) − (5)(2) =
3 84 42
Exercice 3

x3 + x2 + x − 1 = 0 dans [0; 1].

En utilisant l'interpolation inverse de Lagrange aux points 0; 12 et 1,
donner une approximation de l'unique solution de cette équation.
14
Résultats

α existe et est l'unique solution de l'équation f (x) = 0 donc f (α) =

0.
Si p interpole f −1, alors α est approché par f −1(0).
Nous obtenons la table suivante :

y i xi f [yi; yi+1] f [yi; yi+1; yi+2] 0 − yi

−1 0 1

4
7
4 4
1 1 −
− 17 7 = − 40 1
8 2 2+1 357 8
4
17

2 1 −2

4 40 1 195
p(0) = 0 + (1) − (1)( ) =
7 357 8 357
195
donc α =
357

2.3 Interpolation de Hermite

Soit x0, x1, . . . , xk , (k+1) points distincts de [a; b]. Soit α0, α1, . . . , αk ∈
I . Soit f une fonction dénie sur [a; b] telle que ∀i ∈ {0, . . . , k},
N
k
f est αi fois dérivable en xi. Soit n = k + αi. Alors il existe un
X

i=0

15
unique polynôme de dégré n vériant :
p(j)xi = f (j)xi, ∀j ∈ {0, 1, . . . , αi}, ∀i ∈ {0, 1, . . . , k}
Remarque

Le système d'équations linéaires p(j)xi = f (j)xi, ∀j ∈ {0, 1, . . . , αi}, ∀i ∈

k
{0, 1, . . . , k} contient (αi + 1) équations pour n + 1 inconnues.
X

i=0

f (j−1)(xi)
f [x
| 1, .{z
. . , x}1] =
(j − 1)!
j fois
Exemple 4 Déterminons le polynoôme d'interpolation de Hermite
x 0 1
f (x) −1 0
de la fonction f dénie par : 0
f (x) −2 10
f 00(x) × 40
Nous dressons la table de diérence divisée correspondante comme
suit :
x f (x) f [yi; yi+1] f [yi; yi+1; yi+2] f [yi; yi+1; yi+2] f [; ; ; ]
0 −1 1

0 1 −2
3
1 0 1 6
9 5
1 0 10 11
20
1 0 10
p(x) = −1−2(x−0)+3(x−0)2 +6(x−0)2(x−1)+5(x−0)2(x−1)2
p(x) = −1 − 2x + 2x2 − 4x3 + 5x4
16
Exercice 4 L'observation d'un mobile est consignée sur la table
suivante où x désigne sa position sur une droite
t 0 3 5 8 13
x(t) 0 225 385 623 993
ẋ(t) 75 77 × 74 72
3 7 3
ẍ × − ×
5 5 5
Donner une estimation de la position du mobile, ainsi que sa vitesse
à la date t = 10.
Exercice 5 Construire à l'aide d'une feuille excel, la table des dif-
férences divisées de la fonction de la table
x −2 −1 0 1 2 3
y 1 4 11 16 13 −4

17
Théorème 1 Soit x0, x1, . . . , xk , (k + 1) points distincts de [a; b].
k
Soit (α0, α1, . . . , αk ) ∈ N
X
(k+1)
I .n= αi + k.
i=0
Soit f ∈ C [a; b] et P ∈ Pn, le polynôme de Hermite tel que
n+1

p(j)(xi) = f (j)(xi), ∀j ∈ {0, 1, . . . , αi}, ∀i ∈ {0, 1, . . . , k}

alors ∀x ∈ [a; b], ∃ξ ∈ [xmin; xmax]
k
1 (n+1)
Y
f (x) − p(x) = f (ξ) (x − xi)
(n + 1)! i=0
où 
xmin = min(x, x0, . . . , xk )
xmax = max(x, x0, . . . , xk )

x 0 1
Exemple 5 En considérant la table suivante, on a : f (xi) f0 f1
f 0(xi) f00 f10

x f (x) f [yi; yi+1] f [yi; yi+1; yi+2] f [yi; yi+1; yi+2]

0 f0
f00
0 f0
f1 − f0 (f1 − f0) − f00
1 f1
f10 f10 − (f1 − f0) f10 + f00 − 2(f1 − f0)
1 f1
donc
p(x) = f0 + f00 x + [(f1 − f0) − f00 ]x2 + [f10 + f00 − 2(f1 − f0)]x2(x − 1)

18
Chapter 3
Résolution des systèmes linéaires

3.1 Introduction

La résolution d'un système linéaire de la forme Au = b est au centre

de la plupart des traitements de problèmes numériques, provenant
de divers domaines. Le problème de calcul de reseau électrique ci
dessous en est un exemple.
On considère un reseau électrique simple contenant un certain
nombre de résistances et une force electromotrice (une batterie)
disposées suivant le schéma ci dessous.
En utilisant les loi de Kirchho et d'Ohm on peut établir les
équations qui gouverne ce circuit. En désignant par x1, x2, x3 et x4
les courants de boucle comme le montre la gure on a le système



 15x1 − 2x2 − 6x3 = 300
−2x1 + 12x2 − 4x3 − x4 = 0


 −6x1 − 4x2 + 19x3 − 9x4 = 0
− x2 − 9x3 + 21x4

 = 0
La solution de ce système est
x1 = 26.5 x2 = 9.35 x3 = 13.3 x4 = 6.13
Le but de ce chapitre est de passer en revue quelques techniques
ayant fait leur preuve sur des systèmes du genre même si ceux ci
ont des centaines d'inconnues.
19
 Ces techniques sont classées en deux catégories : les méthodes
directes et les méthodes itératives.
Nous étudirons d'abord les méthodes directes c'est à dire celles
qui conduisent en un nombre ni d'opérations élémentaires à la
solution exacte lorsque l'arithmétique utilisée est à précision in-
nie. Nous aborderons ensuite les méthodes itératives consistant à
rechercher la solution du système comme limite d'une suite vecto-
rielle. Ici la solution est toujours entachée d'erreurs.

3.2 Méthodes directes de résolution de systèmes linéaires

Au = b

3.2.1 Systèmes triangulaires

Ils sont de deux types selon que le matrice A est une matrice tri-
angulaure supérieure ou une matrice triangulaure inférieure. Dans
ces cas on utilise de simples techniques de substitution.
Substitution par la remontée

Le système est de la forme



 a11x1 + a12x2 + . . . ... + a1nxn = b1

 a22x2 + . . . ... + a2nxn = b2
... ... ...






aiixi + ... + ainxn = bi


 ... ... ...





 an−1,n−1xn−1 + an−1,nxn = bn−1
 annxn = bn
Le système est inversible si pour tout i = 1, 2, . . . , n, aii 6= 0.
Ainsi la procédure de remontée démarre par
bn
xn =
ann
20
et par une reccurence regressive on a
 
n
1  X
xi = bi − aij xj  i = n − 1, n − 2, . . . , 1
aii j=i+1

D'où l'algorithme :
real array (aij )n×n , (bi)n , (xi)n
integer i, j, n

real sum

xn ← bn/ann
for i = n − 1 to 1 step -1 do

begin

sum ← bi
for j = i + 1 to n do

begin

sum ← sum − aij xj

end

xi ← sum/aii
end

Le nombre de couples (multiplication - addition) que necessite

2
cet algorithme est de l'ordre de n2

21
Substitution par la descente

Le système est de la forme



 a11x1 = b1

 a21x1 + a22x2 = b2
... ... ...






ai1x1 + . . . + aiixi = bi


 ... ... ...





 an−1,1x1 + ... + . . . + an−1,n−1xn−1 = bn−1
 a x + ... + . . . + an,n−1xn−1 + annxn = bn
n1 1

Le système est inversible si pour tout i = 1, 2, . . . , n, aii 6= 0.

Ainsi la procédure de descente démarre par
b1
x1 =
a11
et par une reccurence progressive on a
 
i−1
1  X
xi = bi − aij xj  i = 2, 3, . . . , n
aii j=1

D'où l'algorithme :
real array (aij )n×n , (bi)n , (xi)n
integer i, j, n

real sum

x1 ← b1/a11
for i = 2 to n do

begin

sum ← bi
for j = 1 to i − 1 do

begin

sum ← sum − aij xj

22
 end

xi ← sum/aii
end

Il necessite le même nombre d'opérations que l'algorithme précé-

dent.
Exemple

Résoudre le système
6x1 = 12
3x1 +6x2 = −12
4x1 −2x2 +7x3 = 14
5x1 −3x2 +9x3 +21x4 = −2

3.2.2 Elimination naïve de Gauss

Le processus d'élimination naïve de Gauss comporte deux étapes :

ˆ la triangularisation qui consiste à transformer le système initial
en un système triangulaire supérieur
ˆ la subtitution par la remontée pour résoudre le système trian-
gulaire supérieur obtenu

3.2.3 Triangularisation

Elle est obtenue en appliquant des transformations élémentaires

successives de Gauss à la matrice initiale. Commençons d'abord
par nous faire une idée précise d'une transformation de Gauss.

23
Transformations de Gauss

Elles ont pour fonction d'annuler les composantes d'un vecteur à

partir d'un certain rang. Soit n, k ∈ N
I ∗ avec k < n. On considère
un vecteur x ∈ KI n tel que
xT = [x1, x2, . . . , xk , xk+1 . . . , xn] , xk 6= 0
En posant
xi
αi = pour i = k + 1, . . . , n
xk

αT = [0, 0, . . . , 0, αk+1 . . . , αn]

et  
1 ... 0 0 0 ... 0
 ... . . . ... ... ... ... ... 
 
0 ... 1 0 0 ... 0 
 
T
Mk = I − αek =  0 ... 0 1 0 ... 0 
 
0 . . . 0 −αk+1 1 . . . 0 
 
 .. ... ... ... . . . ... 
. 
0 . . . 0 −αn 0 . . . 1
on a
Mk x = x − xk α ⇒ (Mk x)T = [x1, x2, . . . , xk , 0 . . . , 0]
Exercice 1 Soit α = [0, . . . , 0, αk+1, . . . , αn]T . Montrer que
si Mk = I − αeTk alors Mk−1 = I + αeTk
Exercice 2 Soit M = Mk . . . M1 où Mi = I − α(i)eTi avec
α(i) = [0, . . . , 0, −li+1,i, . . . , −lni]T
Montrer que eTj α(i) = 0 pour j < i.
Montrer que M est une matrice triangulaire inférieure de diago-
nale unité. Plus explicitement on a
24
  
1
l 1 
 21 
l l32 1 
 31
 .. ...

.

M =

 lk1 lk2 lk3 . . . lkk−1 1


 
 lk+1,1 lk+1,2 lk+13 . . . lk+1,k 1
. ...

 ..


ln1 ln2 . . . lnk 0 ... 1
Montrer que
k
X
M −1 = (I + α(1)eT1 ) . . . (I + α(k)eTk ) = I + α(i)eTi
i=1

Exercice 3 Déterminer la transformation de Gauss M ∈ R

I 3×3
telle que    
2 2
M 3  =  7 
4 8
Théorème 3.2.1 On suppose que A ∈ K I m×n et que pour k <
min(m, n) on a déterminé les transformations de Gauss M1, . . . , Mk−1 ∈
I m×m telles que
K
" #
(k−1) (k−1)
A11 A12 k−1
A(k−1) = Mk−1 . . . M1A = (k−1)
0 A22 m−k+1
k−1 n−k+1

où A(k−1)
11 est une matrice triangulaire supérieure. Si
 
(k−1) (k−1)
akk ... akn
 .. ... 
=.
(k−1)
A22 
(k−1) (k−1)
amk . . . amn
25
et a(k−1)
kk 6= 0 alors les multiplicateurs
(k−1)
aik
lik = (k−1)
i = k + 1, . . . , m
akk

sont dénis et en posant Mk = I−α(k)eTk avec α(k) = (0, . . . , 0, lk+1,k , . . . , lmk )

on a
" #
(k) (k)
A11 A12 k
A(k) = Mk A(k−1) = (k)
0 A22 m−k
k n−k
avec A(k)
11 matrice triangulaire supérieure.

Preuve

Théorème 3.2.2 Soit A ∈ K I n×n une matrice non singulière. On

suppose que l'on a déterminé les transformations de Gauss succec-
cives M1, . . . , Mn−1 telles que
U = A(n−1) = Mn−1A(n−2) = Mn−1 . . . M1A
soit une matrice triangulaire supérieure. Alors on a
Ax = b ⇔ U x = (Mn−1 . . . M1) b
On peut alors utiliser l'algorithme de substitution par la remontée
En résumé l'algorithme de triangularisation ou d'élimination de
Gauss peut s'écrire :
real array (aij )n×n , (bi)n
integer i, j, k, n

for k = to n − 1 do

26
begin

for i=k+1 to n do

begin

for j=k to n do

begin

aij ← aij − (aik /akk )akj

end

bi ← bi − (aik /akk )bk

end

end

Remarque :

Dans cet algorithme nous observons que :

1) le multiplicateur aakk
ik
ne depend pas de j donc peut sortir de
la boucle j et on obtient
real array (aij )n×n , (bi)n
integer i, j, k, n
real xmult
for k= to n−1 do

begin

for i=k+1 to n do

begin

xmult ← aik /akk

aik ← xmult
for j = k to n do

begin

aij ← aij − (xmult)akj

end

bi ← bi − (xmult)bk
27
 end

end

2) les termes aij pour i > j qui devaient être mis à 0 n'ont pas été
mis à jour, opération complètement superue. L'espace occupé par
ces élements pourra être utilisé dans l'algorithme de factorisation
LU
Exercice 4 Résoudre par élimination de Gauss les systèmes :

 3x1 +4x2 +3x3 = 10
(a) x +5x2 −x3 = 7
 1
6x1 +3x2 +7x3 = 15

 3x1 +2x2 −5x3 = 0
(b) 2x −3x2 x3 = 0
 1
x1 +4x2 −x3 = 4
    
1 −1 2 1 x1 1
3 2 1 4   x2   1 
    
(c) =
5 8 6 3   x3   1 
  
4 2 5 3 x4 −1

3.2.4 Factorisation LU

Théorème 3.2.3 Une matrice régulière A ∈ K I n×n possède une

factorisation A = LU où L est triangulaire inférieure à diagonale
unité et U est triangulaire supérieure, si et seulement si toutes les
sous matrices principales A(k) de A sont régulière.
Preuve

28
 Si A = LU , nous pouvons avoir une décomposition par bloc des
trois matrices
     
k A11 A12 L11 0 U11 U12 k
= .
(n − k) A21 A22 L21 L22 0 U22 (n − k)
k (n − k) k (n − k) k (n − k)
où A11 est la sous matrice principale d'ordre k de A c'est à dire
A(k). Le produit par bloc donne
A11 = L11U11 et det A11 = det L11 det U11
L11 est triangulaire inférieure à diagonale unité donc det L11 = 1.
U étant une matrice triangulaire supérieure régulière on a
det U11 = ki=1 U (i, i) 6= 0.
Q

Ainsi  
(k)
det A = det A11 = det U11 6= 0
Réciproquement si pour tout k det A(k) 6= 0 on va montrer par
reccurence que l'algorithme naïf de Gauss s'applique.
En eet pour k = 1 la sous matrice est le scalaire non nul A11.
Donc l'etape 1 de l'algorthme naïf de Gauss est réalisable.
Supposons que nous ayons réalisé l'étape k − 1 de l'algorithme
naï de Gauss. Alors on
A = A1 = L1L2 . . . Lk−1A(k) = RA(k)
où R = L1L2 . . . Lk−1 est une matrice triangulaire inférieure àdiag-
onale unité. Posons B = A(k) Par une decomposition par bloc de
dimension k et n − k on a
     
A11 A12 R11 0 B11 B12
= .
A21 A22 R21 I B21 B22
et A11 = R11B11 avec
det B11 = det A11 = det A(k) puisque det R11 = 1
29
 Etant donnée la structure de B = A(k) alors B11 est est une
matrice triangulaire supérieure avec tous les termes diagonaux non
nuls. En particulier A(k)
kk est non nuls et peut être choisi comme le
k eme
pivot.
Théorème 3.2.4 Si une matrice régulière A ∈ K I n×n possède une
factorisation A = LU où L est triangulaire inférieure à diagonale
unité et U est triangulaire supérieure, alors cette factorisation est
unique.
Exemple

Résoudre par l'élimination naïve de Gauss puis déduire la fac-

torisation de la matrice A du système
    
6 −2 2 4 x1 16
 12 −8 6 10   x2   26 
    
=
 3 −13 9 3   x3   −19 
  
−6 4 1 −18 x4 −34

Exercice 5 En utilisant l'élimination naïve de Gauss, calculer la

factorisation LU de A
 
1
  1 0 0 3
3 0 3
 0 1 3 −1 
 
a) A =  0 −1 3  b) A = 
 3 −3 0 6 

1 −3 0
0 2 4 −6

Exercice 6 On considère
 
2 2 1
A=1 1 1
3 2 1

30
 a) Montrer aue A n'admet pas une factorisation LU avec L une
matrice triangulaire inférieure à diagonale unité et U une matrice
triangulaire supérieure
b)Permuter les lignes pour que la factorisation soit possible.

3.2.5 Elimination de Gauss avec stratégie de pivot

Il s'agit ici de modier l'approche naïve an de prendre en compte

une stratégie de choix des pivots. Il faut par exemple commencer
par résoudre le cas d'inopérabilité où l'étape k est bloqué par A(k−1)
kk =
0. On peut avoir satisfaction, soit à travers une permutation de la
ligne k avec une ligne en dessous, c'est à dire une permutation
d'équations; soit à travers une permutation d'inconnues.
Nous allons regarder essentiellement une stratégie dite de pivot
partiel normalisé. Cette stratégie prend en compte des critères

de stabilité du processus de triangularisation.

étape k = 1
On normalise chaque équation du système, en recherchant tout
simplement le vecteur de normalisation s(1) = [s1, s2, . . . , sn] avec
si = max |aij |
1≤j≤n

On regarde ensuite les rapports

 
|ai1|
i = 1, 2, . . . , n
si
On choisit l'indice j correspondant à la première occurence de la
plus grande valeur de cet ensemble. On procede alors à la permuta-
tion des lignes (équations) indexées 1 et j . L'étape se termine par
l'application de la transformation de Gauss M1. En résumé on a
" #
(1) (1)
A11 A12
A(0) = A → A
g (0) = P A(0) → A(1) = M P A(0) =
1 1 1 (1)
0 A22
31
avec A(1) = A
(0)
I le pivot 1.
11 ∈ K
g
11
Pour l'étape 2 on applique l'étape 1 à la sous matrice
 
(1) (1)
a22 ... a2n
 .. ...
=.
(1)
A22 I (n−1)×(n−1)
 ∈K

(1) (1)
an2 . . . ann

. On construit le vecteur de normalisation s(2) = [s2, . . . , sn] avec

(1)
si = max aij
2≤j≤n

On regarde ensuite les rapports

 (1) 
 ai2 
i = 2, . . . , n
 si 

On choisit l'indice j correspondant à la première occurence de la

plus grande valeur de cet ensemble. On procède alors à la permuta-
tion des lignes (équations) indexées 2 et j . Soit P I (n−1)×(n−1)
f2 ∈ K
cette permutation. On obtient
 
(1) (1)
 . 22 . . .
a a2n 
f f
. ...  ∈ K
.
(1) (1) (n−1)×(n−1)
A22 = f
P2A22 =  I
g
(1) (1)
an2 . . . ann
f f

. En posant P2 = diag(1, f
P2) on a
" #
(2) (2)
A11 A12
A(2) = M2P2A(1) = (2)
0 A22

où A(2)
11 ∈ K
I 2×2
est triangulaire supérieure et A
(2)
I (n−2)×(n−2);
22 ∈ K
M2 étant une transformation de Gauss.

32
 En supposant ce processus réalisé jusqu'à l'étape k − 1, c'est à
dire existence de transformation de Gauss M1, . . . , Mk−1 ∈ K
I n×n
et de permutations P1, . . . , Pk−1 ∈ K
I n×n telles que
" #
(k−1)
(k−1)
A11
A12
A(k−1) = Mk−1Pk−1 . . . M1P1A = (k−1)
0 A22
avec A(k−1)
11 I (k−1)×(k−1) triangulaire
∈K supérieure et A(k−1)
22 I (n−k+1)×(n−k+
∈K
étape k
On construit le vecteur de normalisation s(k) = [sk , . . . , sn] avec
(k−1)
si = max aij
k≤j≤n
On regarde ensuite les rapports
 (k−1) 
 aik 
i = k, . . . , n
 si 

On choisit l'indice j correspondant à la première occurence de

la plus grande valeur de cet ensemble. On procède alors à la
permutation des lignes (équations) indexées k et j . Soit P fk ∈
I (n−k+1)×(n−k+1) cette permutation. On obtient
K
 
](k−1) ](k−1)
a
 . 22 ... a2n 
. ...
.
(k) fk A(k−1)  ∈K (n−k+1)×(n−k+1)
A22 = P = I
g
22 
](k−1) ](k−1)
an2 ... ann
. En posant Pk = diag(Ik−1, P
fk ) on a
" #
(k) (k)
A11 A12
A(k) = Mk Pk A(k−1) = (k)
0 A22
I (k−1)×(k−1) et Mk = In − α(k)eTk avec
Ik−1 ∈ K
h iT
(k)
α = 0, . . . , lk+1,k . . . , ln,k
e e

33
où
](k−1)
a
li,k = ik
e i = k + 1, . . . , n
](k−1)
akk
En résumé, l'élimination de Gauss avec une stratégie de pivot
partiel normalisé consiste à construire des transformations de Gauss
M1, . . . , Mn−1 et des permutations P1, . . . , Pn−1 telles que
I n×n
Mn−1Pn−1 . . . M1P1A = U ∈ K
où U est une matrice triangulaire supérieure.
Posons P = Pn−1 . . . P1 alors
L = P (Mn−1Pn−1 . . . M1P1)−1
est une matrice triangulaire inférieure est l'on a la factorisation
P A = LU
Remarque :

Au lieu de se lancer dans des permutations eectives et laborieuses

des lignes, dans la pratique, on se contentera d'associer aux équa-
tions, un vecteur d'index ` = [`1, . . . , `n]T ∈ NI n que l'on initialisera
à ` = [1, 2, . . . , n]T . Les permutations de lignes se résumeront à des
permutations des éléments de `.
Par ailleurs en pratique les vecteurs de normalisation s ne change
pas sensiblement.
On a alors l'algorithme :
real array (aij )n×n,(bi)n , (`i)n
real array (si )n

integer i, k, n

real r , rmax, smax , xmult

34
for i=1 to n do

begin

`←i
smax ← 0
for j = 1 to n do

begin

smax ← max (smax, |aij |)

end

si ← smax
end

for k=1 to n−1 do

begin

rmax ← 0
for i = k to n do

begin

r ← |a`i,k /s`i |
if (r > rmax) then

begin

rmax ← r
j←i
end

end

`j ↔ `k
for i = k + 1 to n do

begin

xmult ← a`i,k /a`k ,k

a`i,k ← xmult
for j = k + 1 to n do

begin

a`ij ← a`ij − (xmult)a`k j

35
 end

b`i ← b`i − (xmult)b`k

end

end

Théorème 3.2.5 Dans l'algorithme de Gauss appliqué à une ma-

trice A ∈ K
I n×n, la phase de triangularisation a un coût de l'ordre
de n3
3

Remarque :

Il faut observer que dans le processus de triangularisation, plus

un pivot a(k−1)
kk est petit en valeur absolue, plus les coecients mul-
(k−1)
aik
tiplicateurs (k−1) sont grands ce qui amplira l'eet des erreurs
akk
d'arrondi inhérentes à l'utilisation d'une arithmétique à précision
nie.
Ainsi l'objectif d'une stratégie de pivot est de réduire le plus
(k−1)
aik
possible les facteurs d'amplication que sont les coecients (k−1) .
akk
Il faudra tout au moins assurer la condition
(k−1)
aik
(k−1)
≤ 1pour i > k
akk
Cette condition est réalisée avec la stratégie de pivot partiel nor-
malisé.
On peut aussi envisager une stratégie dite de pivot total normal-
isé qui est plus stable du point de vu numérique.

36
3.2.6 Système linéaires spéciaux : factorisations de Cholesky, LDM t
...

Factorisation L-D-M T

Théorème 3.2.6 Si toutes les sous matrices principale de la ma-

trice régulière A ∈ K
I n×n sont régulières alors il existe deux matrices
triangulaires inférieures à diagonale unité L , M et une matrice di-
agonale D = diag(d1, . . . , dn) telles que A = LDM T
Exercice 7 Ecrire l'algorithme de la décomposition A = LDM T
en surchargeant aij par lij si i > j et par mij si i < j .
Exercice 8 1) On cosidère la matrice
 
2 −1 2
A =  2 −3 3 
6 −1 8
a) Trouver la factorisation A = LDU 0 où L est triangulaire
inférieure, D diagonale et U 0 triangulaire supérieure.
b) Utiliser cette décomposition pour résoudre Ax = b avec b =
[−2, −5, 0]T
2) Répeter la question 1) pour
   
−2 1 −2 1
A =  −4 3 −3  b=4
2 2 4 4

Factorisation L-D-L T

Théorème 3.2.7 Si A = LDM T est la décomposition LDM T

d'une matrice régulière symétrique A alors L = M et on a A =
LDLT .

37
 L'algorithme de décomposition de Cholesky pour une matrice
réelle, symétrique, denie positive A sera donné ici en calculant la
matrice triangulaire inférieure puis en surchargeant les termes aij
par gij pour i ≥ k
On a :
real array (aij )n×n
integer i, j, k, n

real sum

for k = 1 to n do

begin

for j=1 to k−1 do

begin

akk ← akk − akj ∗ akj

end
√
akk ← akk
for i = k + 1 to n do

begin

for j=1 to k−1 do

begin

aik ← aik − aij ∗ akj

end

aik ← aik /akk

end

end

38
3.3 Méthodes itératives de résolution de systèmes linéaires

Au = b

Nous aborderons la résolution du système linéaire Ax = b avec un

esprit complètement diérent de celui qui a prévalu lors de la mise
en ÷uvre des méthodes dites directes.
Il s'agit ici de construire une suite de vecteurs x(0), x(1), . . . , x(k), . . .
approchant la solution x du système. Le processus numérique est
conçu de manière à assurer la convergence de la suite vers la so-
lution. Ce processus peut s'arrête dès qu'une précision susante
préalablement xée est atteinte.
De façon générale un algorithme itératif pour le système linéaire
Ax = b est conçu en se donnant une matrice régulière Q ∈ Kn×n,
puis, partant d'un vecteur arbitraire x(0), d'engendrer la suite de
vecteurs x(1), . . . , x(k), . . . par l'équation reccurente :
Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .
En supposant que la suite x(k) k converge vers x∗, du fait de la

continuité des applications linéaires, le passage à la limite quand

k → ∞ dans l'équation
Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .
donne Qx∗ = (Q − A)x∗ + b ⇔ Ax∗ = b. Cette limite est donc la
solution du système initial.
L'algorithme général d'une méthode itérative est le suivant :
array (A)n×n ,(b)n,(c)n, (y)n, (x)n
integer k

begin

x ← x(0)
for k = 1 to kmax do

begin

39
 y←x
c ← (Q − A)x + b
Resolution de Qx = c
Ecrire k, x
if kx − yk <  then

begin

Ecrire "Convergence"
stop

end

end

Ecrire "Nombre maximal d'iterations atteint"

end

Remarque :

Le choix de la matrice régulière Q doit prendre en compte le fait

que le système reccurent
Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .
doit être facile à résoudre en x(k) lorsque x(k−1) est connu. Ce choix
doit être fait pour assurer une convergence éventuellement rapide
du processus.
Théorème 3.3.1 Soit b ∈ RI n et A, Q ∈ RI n×n deux matrices
régulières. Si le rayon spectrale de I − Q−1A satisfait à la con-
dition ρ(I − Q−1A) < 1 alors la suite de vecteurs x(k) dénie par
Qx(k) = (Q − A)x(k−1) + b converge vers x = A−1b pour toute
valeur initiale x(0).
Par dénition la matrice I − Q−1A est dite la matrice itérative
de la méthode.
40
3.3.1 Méthodes itératives standard : Jacobi, Gauss Seidel, méth-

odes de relaxation

Considérons le système Ax = b sous la forme explicite

n
X
aij xj = bi 1≤i≤n
j=1

Nous allons construire explicitement la matrice Q en décom-

posant la matrice A de façons simples dans les cas standards.

Méthode de Jacobi

Lorsque pour tout i on a aii 6= 0, en résolvant l'équation numero i

pour l'inconnue xi, on établit la relation de reccurence :
n
(k) (k−1)
X
aiixi =− aij xj + bi 1≤i≤n
j=1
j6=i

ou
 
n
(k) 1 X (k−1) 
xi = −  aij xj + bi  1≤i≤n
aii j=1

j6=i

On a ainsi formulé point par point la méthode de Jacobi

Dans la suite de cette section on décomposera la matrice A sous
la forme A = D − CL − CU où D est diagonale, CL triangulaire in-
férieure à diagonale nulle et CU triangulaire supérieure à diagonale
nulle avec
D = diag(A) = (aii)
CL = (−aij )i>j
CU = (−aij )i<j

41
 En formulation matricielle, la méthode de Jacobi s'écrit
Dx(k) = (CL + CU )x(k−1) + b
Ainsi Q = D et la matrice itérative de la méthode de Jacobi est
G = I − D−1A = D−1(CL + CU )
L'algorithme de la méthode de Jacobi peut s'écrire :
real array (A)n×n ,(b)n,(c)n, (y)n, (x)n
real sum, diag

integer i, j, k, km ax, n

begin

n ← taille(A)
x ← x(0)
for k = 1 to km ax do

begin

y←x
for i = 1 to n do

begin

sum ← bi
diag ← Aii
if |diag| < δ then

begin

Ecrire "Element Diagonal trop petit"

return

end

for j=1 to n do

begin

if j 6= i then

begin

sum ← sum − Aij yj

42
 end

end

xi ← sum/diag
end

Ecrire k, x
if kx − yk <  then

begin

Ecrire "Convergence"
stop

end

end

Ecrire "Nombre maximal d'iterations atteint"

end

Méthode de Gauss-Seidel

Intuitement on espère améliorer l'approximation en utilisant les in-

connues x1, . . . , xj−1 déjà mis à jour lors de l'évaluation de la j ieme
inconnue dans le processus de Jacobi. On établit alors la relation
de reccurence :
i−1 n
(k) (k) (k−1)
X X
aiixi =− aij xj − aij xj + bi 1≤i≤n
j=1 j=i−1
ou
 
i−1 n
(k) 1 X (k)
X (k−1)
xi = − aij xj − aij xj + bi  1≤i≤n
aii j=1 j=i+1

On a ainsi la méthode de Gauss-Seidel dont la formulation ma-

tricielle est :
Dx(k) = CLx(k) + CU x(k−1) + b
43
c'est à dire
(D − CL)x(k) = CU x(k−1) + b
et la matrice Q = D − CL est la matrice triangulaire inférieure de
A.
L'algorithme de la méthode de Gauss-Seidel s'écrire alors :
real array (A)n×n ,(b)n,(c)n, (y)n, (x)n
real sum, diag

integer i, j, k, km ax, n

begin

n ← taille(A)
x ← x(0)
for k = 1 to km ax do

begin

y←x
for i = 1 to n do

begin

sum ← bi
diag ← Aii
if |diag| < δ then

begin

Ecrire "Element Diagonal trop petit"

return

end

for j=1 to i−1 do

begin

sum ← sum − Aij xj

end

for j =i+1 to n do

begin

sum ← sum − Aij xj

44
 end

xi ← sum/diag
end

Ecrire k, x
if kx − yk <  then

begin

Ecrire "Convergence"
stop

end

end

Ecrire "Nombre maximal d'iterations atteint"

end

Exercice 9 Comparer les rayons spectraux ρ(I − Q−1 J A) et ρ(I −

G A) respectifs des méthodes de Jacobi et de Gauss-Seidel pour la
Q−1
matrice  
4 −1 −1
A =  −1 4 −1 
−1 −1 4

45