Chapitre IV

Systèmes d’Equations Linéaires

Considérons un système d’équations linéaires ( 
  donnés)

   
!! !  
.. .. .. .. (0.1)
. . . .
"#$%&# 

et cherchons sa solution $&&' . Très souvent, il est commode d’utiliser la notation matricielle
( )*+
(0.2)
(1*
0 2
Rappelons que le système (0.2) possède une solution unique si et seulement si ,.-'/ .

Bibliographie sur ce chapitre

Å. Björck (1996): Numerical Methods for Least Squares Problems. SIAM. [MA 65/387]
P.G. Ciarlet (1982): Introduction à l’analyse numérique matricielle et à l’optimisation, Masson.
J.J. Dongarra, C.B. Moler, J.R. Bunch & G.W. Stewart (1979): LINPACK Users’ Guide. SIAM.
D.K. Faddeev & V.N. Faddeeva (1963): Computational Methods of Linear Algebra. Freeman &
Co. [MA 65/271]
G.H. Golub & C.F. Van Loan (1989): Matrix Computations. Second edition. John Hopkins Univ.
Press. [MA 65/214]
N.J. Higham (1996): Accuracy and Stability of Numerical Algorithms. SIAM. [MA 65/379]
A.S. Householder (1964): The Theory of Matrices in Numerical Analysis. Blaisdell Publ. Comp.
[MA 65/262]
G.W. Stewart (1973): Introduction to Matrix Computations. Academic Press.
L.N. Trefethen & D. Bau (1997): Numerical Linear Algebra. SIAM. [MA 65/388]
J.H. Wilkinson (1969): Rundungsfehler. Springer-Verlag.
J.H. Wilkinson & C. Reinsch (1971): Handbook for Automatic Computation, Volume II, Linear
Algebra. Springer-Verlag.
88 Systèmes d’Equations Linéaires

IV.1 Elimination de Gauss

L’élimination dite “de Gauss” a été pratiquée pendant des siècles sans grand tam-tam, notamment
par Newton et par Lagrange (en 1781 dans ses calculs astronomiques). Toutefois, Gauss ayant le
souci de prouver l’existence des solutions pour son principium nostrum des moindres carrés (voir
notices historiques du cours d’Algèbre Linéaire, p. 17) décrit l’algorithme explicitement :
( 3
0 2 0 2
3
Soit donné le système (0.1) et supposons que ,.-'/ . Si , on peut éliminer la

variable dans les équations 4 à 5 à l’aide de l’équation 6 , c.-à-d., on calcule
7
89 
8
pour :  4 && 5 (1.1)

et on remplace la ligne : par 7
8>=

ligne :<; ligne 6 
De cette manière, on obtient le système équivalent
@ @ @ @
?  ? %A ? *?
@ @ @
 ?  %A  ?  *  ?
.. .. .. (1.2)
.@ @. . @
#?  %A &?  *  ?
où @ @
? B #
?   
 ; 7 8
$ 
 ?
@
*# 
?
@
 
; 7
# pour :  4 && 5 (1.3)
0
(si %C2 , on échange la première (DB
0 2
ligne de (0.1) avec une autre ligne pour arriver à *2 ; ceci
est toujours possible car ,.- / ).
Le système (1.2) contient un sous-système de dimension 5E;F6 sur lequel on peut répéter
la procédure @ @ pour éliminer  dans les équations G à 5 . On multiplie la ligne 4 de (1.2) par
7
HI
H? KJ  ?  et on la soustrait de la ligne : . Après 5L;M6 étapes
N ( PORQ N ( ? @ &  ? @ OSQ N ( ? @   ? @ ORQ TQ N ( ? UV @   ? "UV @ O9KW NX Y+O

Z
on obtient un système triangulaire Z Z
[ A *Y#
Z Z
 A  *Y 
.. .. .. (1.4)
. Z . .
&!*Y'
qui se résoud facilement par “back substitution”

Z Z Z
 Y' J & \ 
[ N Y'

_+O J 


; ]"
8^V pour :  5L;M6 && 6  (1.5)
(DC
0 2
Théorème 1.1 Soit ,.- / . L’élimination de Gauss donne
` ( *a X (1.6)
`
où est une matrice de permutation et Z Z Z
6 Z   Z 
7    
ab 6  X  
.. .. .. .. . (1.7)
. . . . Z ..
7 "  7 dc "UV #
6
(
La formule (1.6) s’appelle décomposition LR (left - right) de la matrice .
Systèmes d’Equations Linéaires 89

`
Remarque. Les colonnes et les lignes d’une matrice de permutation sont des vecteurs unité. On
`  e
a ,.- / 6 . Un exemple est
2 6 2 ligne 6 ligne 4
`  2 2  `  
6 ligne 4 ligne 6
2 2 6 ligne G ligne G

Démonstration. Supposons que toutes les permutations nécessaires soient déjà( faites avant que
` (
l’on commence l’élimination des variables (par abus de notation, nous écrivons au lieu de
dans cette démonstration). En utilisant les matrices

6 6
7  2
; 6 6
a9 7f 2  ag 2 7f   
; 6 ; 6 (1.8)
.. .. .. .. .. .. .. ..
. . . . . . . .
7 7
; " 2  2 6 2 ;   2 6
(
le premier pas de l’élimination de Gauss correspond à une multiplication de avec a% , le deuxième
avec a , etc.,
a% (  ( ? @ @
a ( ?  ( ? @    a["UV ( ?
U@
 ( ? "UV @  X 

Par conséquent,
X  N a["UVhaiUkjd&jda%$O[j ( (  N a["UV$a["Uj"&"jda%$O UV j X 
et
N UV
Il reste à montrer que la matrice a de (1.7) est égale à a[UVaiUIj.lj.a%$O . Pour ceci, nous
appliquons la même procédure à la matrice a . La multiplication de a avec a% élimine les éléments
de la première colonne en-dessous de la diagonale,N puis la multiplication avec a élimine ceux de
la deuxième colonne, etc. Finalement, on obtient a["UVa["U%jdjda$O[jdam*no identité, ce qu’il
fallait démontrer.
j ,p-'/  ,.- / abj ` (
Calcul
X du déterminant d’une matrice. La formule (1.6) implique que ,.- /
`  N _
O r
,.- / . Comme ,p-'/ ;q6 , où s est le nombre de permutations dans l’élimination de Gauss,
on obtient Z Z
(  N O r j jd&"j &\
,.-'/ ;t6 (1.9)

Résolution de systèmes linéaires.

(
En pratique,
(
on( rencontre souvent la situation où il faut résoudre
une suite de systèmes linéaires u  , v  v , v vwB v v , etc., possédant tous la même matrice.
Très souvent, on connaı̂t  v seulement après la résolution du premier système.
C’est la raison pour laquelle on écrit, en général, le programme pour l’élimination de Gauss en
deux sous-programmes :
DEC – calculer la décomposition LR (voir (1.6)) de la matrice;
(
SOL – résoudre le système <  . D’abord on calcule le vecteur Y (voir (1.4)), défini par
axYy `  , puis on résoud le système triangulaire X ) Y .

Pour le problème ci-dessus, on appelle une fois le sous-programme DEC et puis, pour chaque
système linéaire, le sous-programme SOL.
90 Systèmes d’Equations Linéaires

( ( @
Coût de l’élimination de Gauss. Pour le passage de à ? , on a besoin de
5L
N ;M6 O  divisions (voir (1.1)) et de
5L;M6 multiplications et additions (voirN (1.3)). 
( ? @
Le calcul de nécessite 5z;E4 divisions et 5z;E4 O multiplications et additions, etc. Comme le
travail dû aux divisions est ici négligeable, le coût total de la décomposition LR s’élève à environ
 f
N  N   g{ i} 5
5L;M6 O  5L;4 O  4  6 |
u opérations
G
(opération  multiplication
@
 addition).
 ? nécessite 5b;~6 opérations (voir (1.3)). Par conséquent, on obtient Y avec
{ N Le calcul de { 
5 
O 
 
 
 
  4  6 5 J 4 opérations. Similairement, la résolution du système (1.4) se fait
 ;€6
en 5 J 4 opérations. { f
En résumé, l’appel{ au  sous-programme DEC nécessite 5 J G opérations, tandis que SOL a
seulement besoin de 5 opérations (sur des ordinateurs sériels). A{ titre
 de comparaison, la
formule habituelle pour le déterminant d’une matrice 5b‚L5 contient 5„ƒ 5 "
J termes.

IV.2 Le choix du pivot

0
Dans l’élimination de Gauss, il faut au début choisir une équation (avec 
8  2 ; cet élément
s’appelle “le pivot”) à l’aide de laquelle on élimine  dans les autres équations. Le choix de cette
équation (choix du pivot) peut-il influencer la précision du résultat numérique, si l’on fait le calcul
sur ordinateur en virgule flottante?

Exemple 2.1 (Forsythe) Considérons le système

U‡
6 †2"2Ij 6 2 j&z 6 ˆ2"2Ij&‰ 6 †2"2
(2.1)
6 ˆ2"2Ij&z 6 ˆ2"2Ij&‰ 4 †2"2
avec pour solution exacte
6 2"ˆŠ"ŠŠ"‹
%  6 †2"2"2 6 2"2"2 6 + g *2d†Š"Š"Š"‹Š"Š"Š"‹%+ (2.2)
2"ˆŠ"ŠŠ"Š 2"ˆŠ"ŠŠ"Š

Appliquons l’élimination de Gauss et simulons un calcul en virgule flottante avec G chiffres signi-
ficatifs (en base 6 2 ).
U‡ 7 g! J Ž ‡ @
a) Si l’on prend g 6 †2"2Œj 6 2 comme pivot, on obtient 6 †2"2j 6 2 , l?  
‡ ‡ @ ‡ ‡
6 †2"2 ;M6 @ †2"2Žj @ 6 2  ;q6 †2"2Žj 6 2 et  ?  4 †2"2 ;16 ˆ22gj 6 2  ;q6 †2"2Žj 6 2 . Par conséquent,
I  ? J  ?   6 † 2"2 (exacte!), mais pour  nous obtenons
U‡
9 N  ; !O J 9 N 6 †2"2 ;M6 †2"2I= 6 ˆ22"O J N 6 †2"2j 6 2 O%Cd 2 

Le résultat numérique, obtenu pour  , est faux.

†2"2 et l’élimination de Gauss donne:
b) Si l’on échange Ules
7 9 ‡ deux @ équations de (2.1), U‡ le pivot est 6 @ U‡
6 †2"2„j 6 2 ,  ?   6 †2"2 ;‘6 @ ˆ22„j @ 6 2  6 †2"2 et  ?  6 †2"2 ;u4 ˆ2"2= 6 ˆ2"2j 6 2  6 †2"2 .
De nouveau, on obtient Ž* ? J  ?   6 †2"2 . Mais cette fois le résultat pour  est
% N  !O J l 9 N † 2"2
; 4 ;M6 ˆ22’= 6 †2"2"O J 6 †2"2 6 †2"2d

Les deux valeurs numériques (pour  et aussi pour  ) sont correctes.

Systèmes d’Equations Linéaires 91

Pour mieux comprendre dans quelle partie de l’élimination de Gauss on a perdu une informa-
tion essentielle, considérons les sous-problèmes (addition, soustraction, multiplication, division)
séparément et étudions leur “condition”.
X  X
Condition d’unN [O
problème. ConsidéronsN une application “ Wn Q n , le problème consistant
à calculer “ pour les données N<    
 &
 &
 h
  \ O . Il est intéressant d’étudier l’influence de
perturbations dans sur le résultat “ [O .

Définition 2.2 La condition ” d’un problème “ est le plus petit nombre tel que
• 

• • N iO N •
; — “ ; “ [O
˜
• 
• – eps • N [O • – ” j eps (2.3)
“
On dit que le problème “ est bien conditionné, si ” n’est pas trop grand. Sinon, il est mal condi-
tionné.

Dans cette définition, eps représente un petit nombre. Si eps est la précision de l’ordinateur
(voir le paragraphe II.5) alors, 
peut être interprété comme l’arrondi de 
. Remarquons en-
core que la condition ” dépend des données 
et du problème “ , mais qu’elle ne dépend pas de
N [O
l’algorithme avec lequel on calcule “ .

Exemple 2.3 (multiplication de deux N hnombres réels) Soient donnés les nombres  et  , con-
sidérons le problème de calculer “ !Ox‰ij& . Pour les deux valeurs perturbées

%E N 6 ™P$O$ I N 6 ™!O • ™_
•

– eps (2.4)
on a [j  ; ij&
 N N
ij 6 ™P$O 6 ™šO ; 6 *™P™™Pijd™
• • • ™š •
Comme eps est un petit nombre, le produit P™ j™š est négligeable par rapport à ™  et on
obtient [j  ; [j&
ij& – 4 j eps  (2.5)

On a donc ”  4 et ce problème est bien conditionné.

N
Exemple 2.4 (soustraction) Pour le problème “ h!O9 ;  , un calcul analogue donne
N N •  •  •  •
; O ;  ; !O  ™P ; ™š
 ;   ;  – • • j eps (2.6)
; 
”
Si sign9 ; sign (ceci correspond à une addition et non pas à une soustraction) on a ” 
 6 ; le
{
problème est bien conditionné.
Par contre, si   la condition ”  N•  •  •  • O J •  ;  • devient très grande et on est
confronté à un problème qui est extrêmement mal conditionné. Pour mieux illustrer l’effet de cette
grande condition, considérons l’exemple numérique
›
9 6  g 6 { 4 J 2
› › pour lequel ” N J › 2"O   6 2"2d
6 4 6
UV
En faisantUVle calcul avec G chiffres significatifs U f (en base 6 2 ), on obtient Œ<2" 6 ŠœŒj 6 2 , 
2d 6 Š 4 j 6 2 et ;  g *
 "
2 ˆ
 "
ž 
2 k
2 j 6 2 . Comme les deux premiers chiffres sont les mêmes pour 
et  , la soustraction les fait disparaı̂tre
N › j › O9*2d U f et on n’a plus qu’un chiffre qui est significatif (le résultat
exact est 6 J 6 4 G"Ÿ"Ÿ 6 j 2 ). On parle d’extinction de chiffres.
92 Systèmes d’Equations Linéaires

7
Explication du choix du pivot. Si  est très grand (ce qui est le cas dans la situation (a) de
l’exemple de Forsythe) alors,
@ { @
  ; 7 $ 
?   ;
7   # ? { 
@ 7 { 7 Ž @  (2.7)
 ? 
  ; $ 
;  l? 

Cette valeur, obtenue pour  , est en général correcte. Mais le calcul de 

9 N  ; !O J
{
nécessite une soustraction qui est très mal conditionnée car   et, à cause de l’extinction
de chiffres, on perd de la précision. 7
La conclusion de cette étude est qu’il faut éviter des 
 trop grands.
Recherche0 partielle de pivot. L’idée est de ne pas se contenter d’un pivot qui soit différent de
zéro (l D2 ), mais d’échanger les équations de (0.1) afin que soit le plus• grand élément (en
( 7
 •
valeur absolue) de la première colonne de . De cette manière, on a toujours – 6 . La même
stratégie est répétée pour les sous-systèmes apparaissant dans l’élimination de Gauss.
› ›"›
Expérience numérique. Pour chaque 5  &œ" nous choisissons 4 2"22 matrices aléatoires
avec coefficients 
  uniformément distribués dans ¡;q6 6¢ et des solutions 
uniformément dis-

tribuées dans £;t6  6#¢ . Alors on calcule en double précision les   pour cette solution exacte. Ensuite

10−1 errmax 10−1 errmax

sans recherche de pivot avec recherche de pivot
10−2 10−2

10−3 10−3

10−4 10−4

10−5 10−5

10−6 10−6

10−7 10−7
matrice quelconque matrice quelconque
10−8 5 10−8 5
10 20 30 40 50 10 20 30 40 50
10−1 errmax 10−1 errmax
sans recherche de pivot avec recherche de pivot
10−2 10−2

10−3 10−3

10−4 10−4

10−5 10−5

10−6 10−6

10−7 10−7
matrice orthogonale matrice orthogonale
10−8 5 10−8 5
10 20 30 40 50 10 20 30 40 50
› › ›› ›"›
F IG . IV.1: Erreurs pour 1 million de systèmes linéaires de dimensions ‚ à ‚
Systèmes d’Equations Linéaires 93

on applique l’algorithme de Gauss, une fois• sans recherche de pivot, et une fois avec recherche de
pivot, en simple précision. L’erreur ¤¥#¦

num ;
ex • de chaque résultat est représentée par un
petit point dans les dessins supérieurs de la figure IV.1. Bien que nous ne soyons pas surpris par les
nombreuses erreurs sans recherche de pivot, quelques cas demeurent inacceptables à droite ; bon
nombre de résultats restent cependant bons !
Faisons une deuxième expérience : une matrice avec 
 uniformément
N ( UVdistribués
N (
dans £;t6  6#¢
pour §u¨€: est complétée par A
y ; 
 , pour assurer que ©  n ; O nK O soit orthog-
onale. Cette matrice est calculée en double précision, le reste de l’expérience continue comme
auparavant (voir les résultats au bas de la figure IV.1). Cette fois-ci il n’y a pas d’exception dans la
bonne performance de l’algorithme de Gauss avec recherche de pivot. Nous allons démontrer ces
observations dans les paragraphes suivantes.

IV.3 La condition d’une matrice

En principe, un problème avec ª données et 5 solutions possède ª«‚)5 coefficients décrivant la
sensibilité de la 5 -ème solution par rapport à la ª -ème donnée. Devant cette myriade de valeurs,
il est parfois préférable d’exprimer la condition par un seul nombre. On réussira cela à l’aide de
normes de vecteurs et de matrices (recherche initiée par A. Turing 1948).
Rappel sur la norme d’une matrice. Pour une matrice à ª lignes et 5 colonnes, on définit
¬ ( ¬
¬ ( ¬  ¬ ( ¬ 
­¤
®­ ¥#]V¦ ¤¯
®°] ¥| ¦ ¬ ¬  (3.1)
( ¬ ( ¬
c.-à-d., la norme de est le plus petit nombre qui possède la propriété
¬ ( ¬ ¬ ( ¬ j ¬ ¬ 
– pour tout ²±bn X  (3.2)
¬ ( ¬ 
Evidemment, dépend des¬ ( normes choisies dans nX et n
X³
. Il y a des situations où l’on connaı̂t
¬
des formules explicites pour . Par exemple, si l’on prend la même norme dans les deux espaces
alors, ¬ ¬ 
pour x
8]V • 
• , on a
³
¬ ( ¬ 9 • 
 • ¶
]V
¤¯¥cµ´µ¦´µ´ c 
]V (3.3)

¬ ¬ g N   ·_
pour
]V • 
• O , on a
¬ ( ¬ I (¸( ¶
plus grande valeur propre de (3.4)
¬ ¬¹  • •
pour ¤¥#¦
8]V cµ´µ´µ´ c  
, on a

¬ ( ¬¹  • 
 • 

8]V¤cµ¥#´µ´µ¦ ´ c ³ ]V (3.5)

¬ ( ¬
La norme
¬ n ¬ 
d’une matrice satisfait toutes les¬ ( propriétés
¬ ( jdº ¬ ¬ j ¬ º ¬
d’une norme. En plus, elle vérifie
6 pour la matrice d’identité et – .
(
Après ce rappel sur la norme d’une matrice, essayons( d’estimer la condition du problème )
 . Pour ceci, considérons un deuxième système linéaire u  avec des données perturbées
A
Hk 
 N 6 ™š
O$ • ™_
 •
– ™_»i

 
N 6 ™_
šO$ • ™š
•
– ™_¼ (3.6)
94 Systèmes d’Equations Linéaires

où ™_» et ™š¼ spécifient la précision des données (par exemple ™_» – eps, ™_¼ – eps où eps est la
¬ j ¬
précision
¬ j ¬¹ de l’ordinateur). Les hypothèses (3.6) impliquent (au moins pour les normes et
) que ¬ ( ( ¬ ¬ ( ¬  ¬  ¬ ¬ ¬
; – ™_»½j ;  – ™š¼ij   (3.7)
¬ ¬
Notre premier résultat donne une estimation de ; , en supposant que (3.7) soit vrai.
( ‘< ( ‘  où (
Théorème 3.1 Considérons les deux systèmes N ( O%¾
linéaires et est une matrice
inversible. Si (3.7) est vérifié et si ™š»zj ” 6 , alors on a
¬ ¬ N( O
¬
;
¬ –
” j N ™_»½™_¼$O
69; ™_»½j ” N ( O (3.8)
N ( OxW¿ ¬ ( ¬ j ¬ ( UV ¬ N( O (
où ” . Le nombre ” s’appelle condition de la matrice .
( ( ) N( ( O t ( N
Démonstration. De  ;  ; ; ; [O , nous déduisons que

u ( UV N( ( O t N PO 
; ; ; ; (3.9)
¬ ¬
Maintenant, prenons la norme de (3.9), utilisons l’inégalité du triangle, les estimations (3.7), –
¬ ¬  ¬ ¬ ¬ ¬ ¬ ( ¬ ¬ ( ¬ j ¬ ¬
; et   – . Nous obtenons ainsi
¬ ¬ ¬ ( UV ¬ ™_»zj ¬ ( ¬ j N$¬ ¬  ¬ ¬ O™_¼j ¬ ( ¬ j ¬ ¬ 
; – ;

Ceci donne l’estimation (3.8).

La formule (3.8) montre que pour ™_»Sj ” N ( OtÀ 6 , l’amplification maximale de l’erreur des
N( O
données sur le résultat est de ” .
N( O (
Propriétés de ” . Soit une matrice inversible. Alors,
N( % O Á 6 (
a) ” pour toute ,
NÂ ( Ox N( O Â 0
b) ” ” pour C2 ,
N ( Ox ¬ (%Ä ¬ ¬ (gÈ ¬
c) ” ­¤
Ã&­ ¥#]V¦ ­¤
Ç­ Å¿]VÆ .

O aux matrices de dimension ª ˜ N(

La 0 propriété (c) permet d’étendre la définition de ” ‚ 5 avec
ª  5 .
( ( UV ¬
¬ n ¬  ¬ Ž ¬ ( ¬ j ¬ ( VU ¬
Démonstration. La propriété (a) est une conséquence de 6  – . La
propriété (b) est évidente. Pour montrer (c), nous utilisons
¬ ( VU ¬ ¬ È ¬ ¬ Ž
( È ¬ UV
¬ ( UV ¬  
¤
®°] ¥#| ¦ ¬ ¬ ¤
Ç°] ¥#| ¦ ¬ (gÈ ¬  Ç#°] ň|Æ ¬ È ¬
¤ 

Exemples de matrices ayant une grande condition. Considérons les matrices É  (matrice de
Hilbert) et Ê  (matrice de Vandermonde) définies par ( Y'k § J 5 )

8UV 
6
 É o

Ëc ]V  Ê 
  Y 
8c ]V 
: § ;*6

¬ j ¬¹
Leur condition pour la norme est donnée dans le tableau IV.1.
¸
Exemples de matrices ayant une petite condition. Une matrice Ì V
U est orthogonale si Ì V
U Ì € n .
¬ ¬ Ž
Pour la norme euclidienne, sa condition vaut 6 car Ì
¬
6 et Ì
¬ Ž
6 (l’inverse Ì  Ì ¸
est aussi orthogonale).
Systèmes d’Equations Linéaires 95

TAB . IV.1: Condition de matrices de Hilbert et Vandermonde

5 2 4 6 8 10 f 12
N ‡ | › Î
” É \
 O 4"Ÿ 4 †‹gj 6 2  4 †ŠIj 6 2\‡ Í G ˆžIj 6 2 Î G  j 6 2 G † ‹Ij 6 2 |
N  O ‹ › †œgj 2
” Ê 6 G Ÿ j 6 2 4 †žIj 6 2 6 †œ’j 6 2Ï 6 †2Ij 6 2

Concernant l’interpolation avec des fonctions splines, nous avons rencontré la matrice (voir le
paragraphe II.10, cas équidistant)
ž 6
(  6 6 ž
Ð . 6. .. 5 (3.10)
6 ..
.
..
.
. .
Ð O . Posons alors N(
 6 . Avec la formule (3.5), on vérifie Ð
Le facteur 6 J ¬ n’influence
( ¬¹ ќ
pas ” ¬ ( UV ¬¹ ( ( ў N nŒÒO
facilement que . Pour estimer ,¬ écrivons

¬ ¹
sous la forme ( UV où n
est l’identité et Ò contient le reste. On voit que Ò  6 J 4 . En exprimant par une série
géométrique, on obtient
f
¬ ( UV ¬¹ 6 ¬ Ò ¬¹ ¬ Ò ¬ ¹ ¬ Ò ¬ ¹ 6 
– ž 6    & –
4
¹ÓN ( O
Par conséquent, ” – G indépendamment de la dimension du système.

IV.4 La stabilité d’un algorithme

Le but de ce paragraphe est d’étudier l’influence des erreurs d’arrondi sur le résultat pour l’élimi-
nation de Gauss. Commençons par la définition de la stabilité d’un algorithme et avec quelques
exemples simples. N
Un algorithme pour résoudre le problème “ [O est une suite d’opérations élémentaires Ô & ,
Ô  (addition, soustraction, multiplication, division, évaluation d’une racine, d’une fonction élé-
mentaire, & ) telle que N N N N N
“ [O Ô  Ô UV & Ô  Ô [O$O&ˆOhO (4.1)
N
En général, il existe beaucoup d’algorithmes différents pour résoudre le même problème “ [O .

N
O
L’amplification de l’erreur, en faisant l’opération Ô , est décrite par la condition ” Ô (voir la
définition dans le paragraphe IV.2). L’estimation
N N $O[j N O[jd"j N \O
” “ O – ” Ô ” Ô ” Ô (4.2)

est une conséquence simple de la définition de la condition d’un problème.

Définition 4.1 Un algorithme est numériquement stable (au sens de “forward analysis”) si
N $O[j N O[jdj N \O N
” Ô ” Ô ” Ô – Const j ” “ O (4.3)
N
où Const n’est pas trop grand (par exemple, Const *Õ 5 O ).
N
La formule (4.3) exprime le fait que l’influence des erreurs d’arrondi durant le calcul de “ [O
n’est pas beaucoup plus grande que l’influence d’erreurs dans les données (qui sont inévitables).
96 Systèmes d’Equations Linéaires

‡ N N T[O$O
Exemple 4.2 Soit ‘ 6 2 et considérons le problème de calculer 6 J 6 . Examinons les
deux algorithmes suivants : ×
Ö× N ½ 6
a) 6 O ; Q N t .
½ 6 Ö 6 O
Toutes ces opérations sont très bien conditionnées (voir le paragraphe IV.3). Ainsi, cet algorithme
est numériquement stable. ×
Ö× 6 J 6 6  6
b) N t ; t N ½ .
6t 6; Q
J 6 O Ö 6 6 O
Dans cet algorithme, seules les trois premières { opérations
N
sont bien conditionnées. La soustraction,
à la fin, est très mal conditionnée car 6 J 6 J z 6 O . Ainsi, cet algorithme est numériquement
instable.
La vérification, si un algorithme (non-trivial) est stable (au sens de “forward analysis”), est
souvent très complexe et difficile. Pour cette raison, Wilkinson (1961, J. Ass. Comp. Mach. 8) a
introduit une autre définition de la stabilité d’un algorithme.
N
Définition 4.3 Un algorithme pour résoudre le problème “ [O est numériquement stable (au sens
Ä
de “backward analysis”) si le résultat numérique
Ä  N iO
peut être interprété comme un résultat exact
pour des données perturbées (c.-à-d., “ ) et si
• 
•
; 

• 
• – Const j eps (4.4)

où Const n’est pas trop grand et eps est la précision de l’ordinateur.
Remarque. Pour l’étude de cette stabilité, il ne faut pas connaı̂tre la condition du problème.
f j‡
Exemple 4.4 Considérons le problème de calculer le produit scalaire ij&kb . On utilise
l’algorithme × ij&
N h!h f '‡!O Ö× [j&R f j&‡
f j&‡ (4.5)
Ö
Le résultat numérique (sous l’influence des erreurs d’arrondi) est
 N 6 ™$O[j& N 6 ™!O N 6 RØ$O˜ f N 6  ™ f O[j&‡ N 6 ™š‡!O N
6 bØO
N f
6 ˜Ø O
• • • • f j ‡
où š™
 Ød – eps. Ce résultat est égal à ij % si l’on pose
N N
x 6 ™PhO 6 ˜ØV$O f ‰ f N 6 ™ f O N 6 RØ"O$
g N 6 ™šO N 6 ˜Ø f O ‡Ž‰‡ N 6 ™‡!O N 6 RØ f O$

Ainsi, (4.4) est vérifié pour Const  4 (on néglige les produits ™_
jHØd ). En conséquence, l’algorithme
(4.5) est toujours numériquement stable (au sens de “backward analysis”).
Cet exemple montre bien qu’un algorithme peut être stable, même si le problème est mal
conditionné. Ainsi, il faut bien distinguer les notions “stabilité numérique” et “condition d’un
problème”.
(
B2 ) ayant la décompo- (€0
La stabilité de l’élimination de Gauss. Soit donnée une matrice (,.-'/
( Da X
sition (on suppose que les permutations nécessairesX sont déjà effectuées). En appliquant
l’élimination de Gauss, nous obtenons deux matrices a et , qui représentent la décomposition
( Wˆ aj X
exacte de la matrice . Pour montrer la stabilité numérique (au sens de • “backward
analysis”) de l’élimination de Gauss, on a besoin de trouver une estimation de la forme 
 ; 
 • –
• 
 • j
Const j eps. En tous cas, il faut estimer la différence 
 ; 
 .
Systèmes d’Equations Linéaires 97

( X
Théorème 4.5 (Wilkinson) Soit une matrice • 7
 •
inversible et a , le résultat numérique de l’élimination
de Gauss (avec recherche de pivot, c.-à-d. – 6 pour tout :  § ). Alors,
• 
 
 • – N
; 4 jdŽj ¤ÅˆÆ :i;M6  § Oj eps (4.6)
• Ù!@ •
où o ¤¯¥¦ 8
c  c Ù l?
H .
Ù!@ Ù!UV @
Démonstration. Lors de la Ú ème étape de l’élimination de Gauss, on calcule
?  à partir de
?  .
Si l’on prend en considération les erreurs d’arrondi, on obtient
Ù!@
N 
?  Ù!UV @ Ù!UV @ N

?   ;
ِj
7 Ù ?  j 6 ™_
HšÙO$O N 6 ˜Ød
HšÙO
Ù!UV @ 7 Ù!UV @ (4.7)

?  ;
ِj Ù?  R[Û
šÙ
N 
où (en négligeant les termes Õ eps O )
• Û[
HšÙ • •
?  Ù@ • j • Ød
Ù •  • 7 Ù!UV @ • •

Ù • j • Ù ?  j ™_
Ù • – 4 jdgj eps 
– (4.8)

Ëc @ 7 Z
Ëc $@ 7 Ù!UV @
(
Par définition de , on a 
Œ ÜÞ
Ù!]Vݵß?
َj Ùg ÜÞ
Ù!]Vݵß!?
ÙIj Ù ?  et, en utilisant la formule
(4.7), on obtient pour :¨E§
 

y N
? Ù!UV @
? 
Ù!@
RÛ[
šÙOx*
„ Û[
šÙ
Ù!]V ; Ù!]V (4.9a)
@
car l?
H B2 dans cette situation. Pour : – § , on a

8UV
UV

y N
? Ù!UV @ 
? 
Ù!@
˜Û
šÙ!O 7

j 
? 

8UV @
*
„ Û[
Ù
Ù!]V ; Ù!]V (4.9b)
7
car H

Þ 6 . Les formules (4.9) ensemble avec l’estimation (4.8) démontrent l’estimation (4.6).

Conséquence. L’élimination de Gauss est stable (au sens de “backward analysis”) si le quotient
• H
?  Ù!@ • • 
 •
¤

Ëc  ¥#c Ù ¦ ¤¯
Ëc¥ ¦ (4.10)

n’est pas"UVtrop grand. En fait, il existe des matrices

(
pathologiques pour lesquelles ce quotient
atteint 4 , où 5 est la dimension de la matrice (voir exercice 12). Mais heureusement cette
constante est en général beaucoup plus petite. Pour illustrer ceci, en suivant une idée de Trefethen
& Schreiber (1990, SIAM J. Matrix Anal. Appl. 11) nous avons pris un grand nombre de matrices
de dimensions allant de 4 à 4 ž dont les éléments sont des nombres aléatoires dans £;t6  6#¢ . Nous
avons dessiné dans la figure IV.2 le quotient (4.10) en fonction de la dimension de la matrice
(chaque point représente la moyenne de G 2 échantillons). En échelle doublement logarithmique, le
résultat semble mystérieusement devenir une droite.

ž quotient (4.10)
G

5
100 101
F IG . IV.2: Stabilité numérique de l’élimination de Gauss
98 Systèmes d’Equations Linéaires

IV.5 L’algorithme de Cholesky

Etudions l’élimination de Gauss pour le cas important où
( ( ¸  (
est symétrique ( ¸ ( ) et 2
( 0
est définie positive ( ¨ pour B2 ).
Le théorème suivant montre que, dans cette situation particulière, il n’est pas nécessaire d’effectuer
une recherche de pivot.
(
Théorème 5.1 Soit une matrice symétrique et définie positive.
a) L’élimination de Gauss
(
est
X
faisable sans recherche de pivot.
b) La décomposition *a satisfait
Z Z
X  à½a ¸ avec àá ,.ň¥"â
N # & dO$ (5.1)
¸ ( N & 2d&&2"O ¸ (
Démonstration. a) On a y ¨ 2 (avec y 6
est définie ) car la matrice
positive. Alors, on peut choisir comme pivot dans la première étape de l’élimination de Gauss.
Ceci donne ¸ l ¸
(  ( ? @ 
ã ; Q 2 ã ?
@ (5.2)
@ 
8
où Y#
?  BA
H ; j"  pour :  §  4 # 5 , ce qui est équivalent à

@
ã ?  ã ;
6 jdŽjd ¸  (5.3)
l
@
La matrice ã ? Ä est symétrique
X "UV ÄRB
(trivial). Montrons qu’elle
0 2 Ä ¸
est
@ Ä aussi définie positive. Pour ceci,
nous( prenons un ±mn , . Il faut montrer que ã ? ¨ 2 . La partition de (5.2) et le fait
que soit définie positive impliquent
¸  
N  Ä ¸ O Bl   Ä ¸ g Ä ¸ ã Ä
ã Ä 4 ij ¨ 2d (5.4)

Ä ¸
En posant x ; J dans (5.4), on obtient de (5.3) que

Ä ¸ ã ? @Ä  Ä ¸ ã Ä 6 N Ä ¸ dO  2d
; ¨

Par récurrence, on voit que la deuxième et aussi les autres étapes de l’élimination de Gauss sont
faisables sans recherche de pivot.
b) La formule (5.1) est une conséquence X ä de l’unicité
¸ de l’élimination de Gauss pour des matri-
(  ( ¸  X ¸ a ¸  a N à½a ¸ O
ces inversibles. En effet, on peut écrire à a et on obtient ,
d’où ˜
a *
 a .
( åa X  a X
Pour montrerUV l’unicité
XX V
U de l’élimination de Gauss, supposons et considérons
l’identité a 
a  . Le produit de deux matrices triangulaires inférieures reste une matrice
triangulaire inférieure; UV de même pour les matrices triangulaires supérieures. Comme les éléments
de la diagonale de a a sont tous égaux à 6 , on a
UV XX UV *n 
a ab (5.5)

ce qui implique l’unicité de l’élimination de Gauss.

Systèmes d’Equations Linéaires 99

La décomposition
( *axà½a ¸ (5.6)
Z
(
s’appelle décomposition rationnelle ·_ Z
deNçæ Cholesky .æ Comme
H
¨ 2 ( est définie positive),
1 Z
on
( 
peut considérer la racine à  ,.Å¿¥â h  
 
 #
  #
  O , et la décomposition (5.6) devient
N axà ·_ O N à ·_ a ¸ O¯ N axà ·_ O N axà ·_ O ¸ ·_
. Par abus de notation, en écrivant a pour axà , nous
obtenons la décomposition de Cholesky 7 2
(  a a ¸
x où ab .. ..  (5.7)
. .
7  & 7 #

Une comparaison des coefficients dans l’identité axa ¸ donne pour

( 

7 Ù  7 Ù$ è 7 ÙÙ

:  Ú W Ùٌ
7
8 7 Ù!Þ&è 7
8c Ù!UV 7 Ùc Ù!UVé 7
HÙ 7 Ù Ù
:¨*Ú W 
ُ

et on en déduit l’algorithme suivant:

Algorithme de Cholesky.
for 7 Ú W¿ 6 Nto 5 do Ù!UV 7  ·_
ÙٌWˆ ÙÙ ; $]V Ù O ;
for 7 : ˆ
W  Ú  6 to 5 Ù!doUV

ٌWˆ N A
HÙ ; $]V 7
 7 Ù PO J 7 Ù Ù .

Le coût de cet algorithme. En négligeant les 5 racines, le nombre d’opérations nécessaires est
  f
N { N } 5
O[j Ú [Oç ) 
Ù!]V 5L;êÚ | 5L; œ

Ceci correspond à la moitié du coût de la décomposition LR.

(
Pour résoudre le système Ñ« , on calcule d’abord la décomposition
¸ de Cholesky (5.7).
Puis, on résoud successivement les deux systèmes axY)ë et a DëY , dont les matrices sont
triangulaires.
Comme pour l’élimination de Gauss, on peut étudier la stabilité de l’algorithme de Cholesky.
(
Théorème 5.2 Soit une matrice symétrique et définie positive. Notons
(  a²j a ¸ , où a est la
matrice triangulaire obtenue par l’algorithme de Cholesky. Alors,
• 
 
 • – | j ¤ÅˆÆ N :  § O[j eps
; (5.8)
• • • •
où |  ¤¥#¦
8c  
  ¤¥#¦
A
H
.

Ce résultat démontre que l’algorithme de Cholesky ( est toujours numériquement stable. Il n’est
donc pas nécessaire de faire une recherche de pivot, si est symétrique et définie positive.
1
Le “Commandant Cholesky” (1875–1918) entra à l’École Polytechnique à l’âge de vingt ans et en sortit dans
l’arme de l’Artillerie. Affecté à la Section de Géodésie du Service géographique, en juin 1905, il s’y fit remarquer
de suite par une intelligence hors ligne, une grande facilité pour les travaux mathématiques, un esprit chercheur,
des idées originales, parfois même paradoxales, mais toujours empreintes d’une grande élévation de sentiments et
qu’il soutenait avec une extrême chaleur. (...) Cholesky aborda ce problème en apportant dans ses solutions, ... une
originalité marquée. Il imagina pour la résolution des équations de condition par la méthode des moindres carrés un
procédé de calcul très ingénieux ... (copié du Bulletin géodésique No. 1, 1922).
100 Systèmes d’Equations Linéaires

IV.6 Systèmes surdéterminés – méthode des moindres carrés

Considérons un système d’équations linéaires
ìältí9í AB
!ìä tí9í! AB 
.. .. .. .. (6.1)
. . . .
³  ³ +%&è ³ AB ³
( X  X ³ (
où ª Á 5 (matriciellement:   avec ±în et R±Fn est une matrice ª ‚m5 ).
;
Evidemment, le système (6.1) ne possède, en général, pas de solution. L’idée est de chercher un
vecteur tel que ¬ ( ¬
;  ŽQ ¤Å¿Æ (6.2)
pour la norme euclidienne. Une justification probabiliste de cette condition sera donnée dans le
paragraphe IV.8. Le nom “méthode des moindres carrés” indique le choix de la norme dans (6.2)
(la somme des carrés des erreurs doit être minimale).
( Á X ³
Théorème 6.1 Soit une matrice ªï‚L5 (avec ª 5 ) et soit ±bn . Le vecteur est solution
de (6.2) si et seulement si ( ¸ ( ) ( ¸ +
(6.3)
Les équations du système (6.3) s’appellent “équations normales”.
Démonstration. Les minima de la fonction quadratique

N
Ô [O%Wˆ
¬ (
;  ¬  N(
; PO ¸ N ( ; PO%E ¸ ( ¸ ( ;4
¸ ( ¸  ¸ 

N N ¸ ( ¸ ( ¸ ( O.
sont donnés par 2ð Ô v [Ox 4 ;
˜±‘n €ò ( • X  ó X ³
Interprétation géométrique.
Xk³ L’ensemble ñ est un sous-espace
(
linéaire de n .
Pour un ²
 í
± n arbitraire, est une solution de (6.2) si et seulement si  est la projection
orthogonale de  sur ñ . Ceci signifie que ( ; Œô (gÈ pour tout È ±án X . On en déduit que
( ¸ N(
; POC2 et on a ainsi établi une deuxième démonstration de (6.3).

Exemple 6.2 Pour étudier le phénomène de la thermo-électricité, on fait l’expérience suivante. On

soude un fil de cuivre avec un fil de constantan{ de manière à obtenir une boucle fermée. Un point
de soudure est maintenu à température fixe (õ | 4 ždö C), alors que l’on fait varier la température õ
de l’autre. Ceci génère une tension Ì , laquelle est mesurée en fonction de õ (voir le tableau IV.2
et la figure IV.3). Les données du tableau IV.2 sont prises du livre de P.R. Bevington2.
On suppose que cette dépendance obéit à la loi

Ì  g õ Y õ (6.4)
et on cherche à déterminer les paramètres \ et Y . Les données du tableau IV.2 nous conduisent
au système surdéterminé (5  G , ª  4V6 )

Ì
[*g õ
VY õ
 :  6 && 46  (6.5)
›
En résolvant les équations normales (6.3) pour ce problème, on obtient z ; 2dˆ‹‹"œ , ÷2dˆ2 G 4
Ó
Y í
 d
2  › Š"‹Kj 2 U‡
et 6 . Avec ces paramètres, la fonction (6.4) est dessinée dans la figure IV.3. On
observe une très bonne concordance avec les données.
Remarque. Les équations normales( ¸ (6.3)
(
possèdent toujours au moins une ¸ ( solution
¸ ( ) (la projection
¬ ( ¬  Á*2
sur ñ existe toujours). La matrice est symétrique et non-négative ( ).
2
P.R. Bevington (1969): Data reduction and error analysis for the physical sciences. McGraw-Hill.
Systèmes d’Equations Linéaires 101

TAB . IV.2: Tensions mesurées en fonction de la température ø

: õ
ö C Ì
: õ
ö C Ì
: õ
ö C Ì

2 2 ˆ‹"Š
" ‹ › 2 †ž 4
d ›
6 ; G 6 Ÿ 2 6 ˆ‹"‹
› 2"ˆœ"Š Š ž2 2d†œ 6 ›
4 ; 6 œ Ÿ 4 6 2
G 6 2 ; 2" › G 6 2 ž › 2d†‹ 4 #6 Ÿ ‹"2 4  G6
ž › 2" G ž › 2 ‹ › ›
6 ; 6 6
" 6 †2 G 6 ‹ 4  ž
› 2" 6 › ››
4 2 ; 6#4 6  4"4 6 Š Š"2 4 Ÿ ‹
œ › 2d†2 4 œ2 › Š ›
4 6#G 6 †ž 4 2 G ˆ2"2
Ÿ G 2 2d 4 2 6 ž œ › 6 †œ"‹ 46 6 2"2 G  4"4

3
Ì

1
õ | õ
0

õ
0 50 100
F IG . IV.3: Tension en fonction de la température et schéma de l’expérience

*2 pour *2 ). ( ( 0 0

Elle est définie positive si les colonnes de sont linéairement indépendantes (
Dans cette situation, on peut appliquer l’algorithme de Cholesky pour résoudre le système (6.3).
Mais, souvent, il est préférable de calculer la solution directement de (6.2) sans passer par les
équations normales (6.3).

IV.7 Décomposition QR d’une matrice

(
Dans l’élimination de Gauss, on a multiplié l’équation ‘~ parX la matrice triangulaire a["UVùj
lj.aIj.a% . De cette manière, on a réduit le problème original à äY où X est une matrice
(
triangulaire supérieure. Malheureusement, la multiplication de ;  avec a[
ne conserve pas la
norme du vecteur.
Pour résoudre (6.2), nous cherchons une matrice orthogonale © telle que
X v Yv
¸ N( +Ox X Yy
© ; ; 2 ; Y vv (7.1)
X Y v Y v v O est la partition N ¸
où v (une
¸  matrice carrée deX  dimension X5 ³ )U"est triangulaire supérieure et
de K
Y  © telle que Y v ê
± n et Y v v ê
± n . Comme le produit par une matrice orthogonale ne
change pas la norme du vecteur, on a
¬ (  ¸ N(    ¬ Y v v ¬  
;  ¬   ¬
© ; PO ¬   ¬&X
; Y ¬   ¬&X v
; Yv ¬   (7.2)
On obtient alors la solution de (6.2) en résolvant le système
X v u*Y v 
(7.3)
102 Systèmes d’Equations Linéaires

¸ ïn ) et une matrice

Le problème consiste àX calculer une ¸ matrice
(  X
orthogonale © (c.-à-d., © ©
triangulaire supérieure telles que © ou de façon équivalente
(  X 
© (7.4)
(
Cette factorisation s’appelle la “décomposition QR” de la matrice . Pour arriver à ce but, on peut
se servir des rotations de Givens (voir exercice 11 du chapitre V) ou des réflexions de Householder.
Réflexions de Householder (1958). Une matrice de la forme

É *n ;4&úú ¸ où ú

¸
ú  6 (7.5)

a les propriétés suivantes : É

ó
û É est une réflexion à l’hyper-plan &ò • ú ¸ uB2
N ¸ ú
car É ) ;bú j 4#ú [O et É tR‘ô ú .
û É est symétrique.
û É est orthogonale, car
¸ Nn ¸ ¸ N ¸ ž úú ¸ ž úú ¸ úú ¸ 
É É  ;4#úú O n ;4&úú O9 n ; n\
( (
En multipliant avec des matrices de Householder, nous allons essayer de transformer en une
matrice de forme triangulaire.
L’algorithme de Householder
¸ -XBusinger
³ ¸ 9- Golub. Dans une première étape, on cherche une
matrice É 9 * n ;4#ú ú (ú %± b n et ú ú 6 ) telle que
 j!j!j
‚ ‚
( 2 ‚ j!j!j ‚
É   (7.6)
.. .. ..
. . .
2 ‚ j!j!j ‚
( ( ( y   N &2d&&2"O ¸
Si l’on dénote
•¡Â •  par
¬ ¬lagpremière colonne de , il faut que É et on
obtient É
(  ¬ ( ¬  . La forme particulière de É implique que
( % ( ¸ ( x Â
É ;4&ú j ú 
¸ (
L’expression ú est un scalaire. Par conséquent,
( Â
ú xBãj&ü. où ü.% ; (7.7)
¬ ¬ K
et la constante ã est déterminée par ú 6 . Comme on a encore la liberté de choisir le signe
de , posons  9 N ¬ ( ¬ 
; sign l $O[j (7.8)
( Â
pour éviter une soustraction mal conditionnée dans le calcul de ü.% ; .
( ( N ( O_ ( (
Calcul de É . Notons par  et É les § èmes colonnes de et É respectivement. Alors,
on a
N ( O_ (  ¸ ( y (  ¸ ( ùj#ü. 4 
É ;4&ú ú ;˜ý j&ü où ý 
ü ¸ .
ü (7.9)

Le facteur ý peut être calculé à l’aide de

UV ü ¸ ü. 6 ( ¸ ( Â $  Â  Â N Â hO$
ý   ;4  ; ; (7.10)
4 4
Systèmes d’Equations Linéaires 103

Dans une deuxième étape, on applique la procédure précédente UV à la sous-matrice de dimension

N N
ªþ;Ñ6 O ‚ 5²;Ḡ6 O de (7.6). Ceci donne un vecteur ÿú t±*n X ³ et une matrice de Householder
N 2d &O ¸
ÿÉ ˜3n ;14Vÿú  ÿú  . En posant ú R ÿú , une multiplication de (7.6) par la matrice É ˜
n ;4#ú  ú ¸ donne

   j!j!j
j!j!j j!j!j ‚
‚ 
‚ ‚ ‚ ‚ ‚
2 2 2 ‚ j!j!j ‚
(
É  É  É  .. ã
 ..  2 2 j!j!j 
. . É ÿ !ã W W ‚W ‚W
2 2 2 2 j!j!j
‚ ‚

En continuant cette procédure, on obtient après 5 étapes (après 5;‰6 étapes si ª  5 ) une matrice
triangulaire X
(  X  v
É Œjd&j É  É 2

¸
©
¸ 
Ceci donne la décomposition (7.4) avec © É Œj"j É  É .
 {
Coût de la décomposition QR. ¸ La première étape exige le calcul de par la formule (7.8)N ( ( ª
opérations), le calcul de 4 J ü üV par la formule { N (7.10) (travail négligeable) et le calcul de É Oš
pour §  4 # 5 par la formule (7.9) ( 5q;6 Oj 4 j ª opérations). En tout, cette étape nécessite
environ 4&ªz5 opérations. f
Pour la{ décomposition QR, on a alors besoin de
N  N 
4 5  5L;M6 O è 6 O { 4&5 J G opérations si ª  5 (matrice carrée);
N N
4&ª 5  5L;M6 Oè 6 O ªz5 opérations si ª 5 . { 
J 4


En comparant encore ce travail avec { celuif de la résolution des équations normales ( z
ª 5
( ¸ ( J œ
opérations pour le calcul de et 5 opérations pour la décomposition de Cholesky de
( ¸ (
), on voit que la décomposition QR coûte au pire le double.
( Â
Remarque. Si les colonnes de la matrice sont linéairement indépendantes, tous les
sont
non nuls et l’algorithme de( Householder–Businger–Golub est applicable. Une petite modification
(échange des colonnes de ) permet de traiter aussi le cas général.
Concernant la programmation, il estÂ

important de ne calculer ni les matrices É
, ni la matrice
© . On retient simplement les valeurs et les vecteurs ü
(pour :  6 &# 5 ) qui contiennent déjà
toutes les informations nécessaires pour la décomposition. Comme pour l’élimination de (
Gauss,
on écrit
 deux sous-programmes.
X DECQR fournit la décomposition ¸ QR de la matrice (c.-à-d.
les
, ü"
et X la matrice ). Le sous-programme ¸ SOLQR calcule ©  et la solution du système
triangulaire v TáY v (voir (7.3)). Le calcul de © Œ É ðj  j É  É h se fait avec une formule
analogue à (7.9).

(
Exemple 7.1 Si les colonnes de sont “presque” linéairement dépendantes, la résolution du
problème (6.2) à l’aide de la décomposition QR est préférable à celle des équations normales.
Considérons, par exemple,
6 6 6
(  ™ 2  y 2
2 ™ 2

où ™ est une petite constante, disons ™ ¾ eps. Avec un calcul exact, on obtient

( ¸ (  6 ™ 6   ( ¸ k 6
6 6 
 ™ 6
104 Systèmes d’Equations Linéaires

et la solution est donnée par


6
 
6   Õ N ™ O
xg
4 
 ™ 4
 ( ¸ (
Un calcul en virgule flottante fait disparaı̂tre le ™ dans et cette matrice devient singulière. On
n’obtient pas de solution. 
Par contre, l’algorithme de Householder–Businger–Golub donne (en négligeant ™ )  ’ ;t6 ,
ü.x N 4 ™#2"O ¸ , & et à la fin

;t6 æ ;q6 ¸ k ;qæ 6

X  2 4 jd™  © ™ J 4
æ

2 2 ; ™ J 4

La résolution de (7.3) donne une bonne approximation de la solution exacte.

IV.8 Etude de l’erreur de la méthode des moindres carrés

Comme c’est le cas dans l’exemple du paragraphe IV.6, nous cherchons les paramètres #h
d’une loi 

$]V
Y' N O"yB  (8.1)
    $UV
   
N N
qui relie les variables et  (les fonctions hY  O sont données, p. ex. Yh Oy ). Supposons que
pour plusieurs valeurs de (disons  & 
 
 
  ³ , ª ï5 ) l’on puisse mesurer les quantités # ³ .
On obtient ainsi le système surdéterminé


 yB
_ :  6  ª (8.2)
]V


N
où 
˜ Y'
O . En pratique, les 
sont des mesures légèrement erronées et il est naturel de
les considérer comme des valeurs plus ou moins aléatoires. L’étude de l’erreur de la solution ,
obtenue par la méthode des moindres carrés, se fait alors dans le cadre de la théorie des probabilités.


Rappel sur la théorie des probabilités


Considérons desX variables
X
aléatoires (dites “continues”) qui sont spécifiées par une fonction
de densité Ô ù
W n Q n , c.-à-d., la probabilité de l’événement que la valeur de se trouve dans
l’intervalle &+O est donnée par
¼

¹
` N
–  ¾BPOx
 N
Ô [O
}
(8.3)
N [O%Á*2 X ¹ N }
pour u±mn U Ô [O )
avec Ô

6 . et
On appelle espérance (mathématique) de la variable aléatoire le nombre réel


¹
Û  ñ
N
 O9
U ¹
N }
Ô [O i (8.4)

et variance la valeur

¹ ¹
s  
 Var
N
 O9
U ¹
N
;
}
Û ’O Ô N [O )

U ¹
 N
Ô [O
}
; Û

  (8.5)
Systèmes d’Equations Linéaires 105

Exemple 8.1 Si une variable aléatoire satisfait (8.3) avec (voir la figure IV.4)
Û 
N [O% 6 j -'¦ 6 ;
Ô æ ; (8.6)
4 j s 4 s

alors on dit que laN variable

 aléatoire satisfait la loi normale ou la loi de Gauss
 – Laplace que l’on
symbolise par Ò 9
Û  O Û
s . On vérifie facilement que est l’espérance et s la variance de cette
variable aléatoire.

La loi normale est parmi les plus importantes en probabilités. Une raison est due au “théorème
de la limite centrale” qui implique que les observations pour la plupart des expériences physiques
obéissent à cette loi.

 %

Û ;4ds Û ;s Û Û s ÛL 4ds

F IG . IV.4: Fonction de densité pour la loi normale

Rappelons aussi que 5 variables aléatoires  &  

 sont indépendantes si, pour tout 
š
,

 
on a
` N 

[¾C'
š :  ` N 

– 6 && 5 O9
8]V
–
[¾*
çO (8.7)

  
Lemme 8.2 N
Soient
NÄ O
et deux variables  aléatoires indépendantes avec comme fonctions de
Ô [O et  ý ±mn X avec  *
0 2

  
densité
Â
respectivement et soient . Alors, les variables aléatoires
 ý et  possèdent les fonctions de densité
¹
6
•¡Ây• Ô
;bý
 et
N
Ô =  O N È Ox
U ¹ Ô
NÈ
;
Ä O
 NÄ O } Ä 
(8.8)

Leur espérance mathématique est

ñ
NÂ
  ý O9
Â
ñ
N
 Oé ý  ñ
N
   Ox ñ
N
 O ñ
N
 O (8.9)

et leur variance satisfait

Var
NÂ
  ý O9
 
Var
N
 O Var
N
   Ox Var
N
 O Var
N
 O (8.10)

     ý découle de (pour  ¨ 2 )

Démonstration. La fonction de densité pour la variable aléatoire

` N
–
Â
  ý ¾*POx ` ;bý
 –  ¾
 ;˜ý
  ?
¼ U "@£·

  
U "@£· Ô
N [O } )
¼
 UV
Ô
;bý
Â
}

 
?

 
Les propriétés (8.9) et (8.10) pour  ý en sont une conséquence directe.
È
Comme et sont supposées indépendantes, on obtient (en posant Et )
Ä
Â

¼ ¹
` N
  –  ¾*PO%
  ® ^ Ã ¼ Ô
N [O
 NÄ O } } Ä

 U ¹ Ô
NÈ
;
Ä O
 NÄ O } Ä } È
106 Systèmes d’Equations Linéaires

et on trouve la fonction
¹ ¹ de densité pour    . Un
¹ calcul
¹ direct donne

ñ
N
   Ox
U ¹
È
U ¹ Ô
NÈ
;
Ä O
 NÄ O } Ä } È 
U ¹ U ¹
N Ó Ä O N [O
Ô  N Ä O } Ä } )
ñ
N
 OV ñ
N
 O

et, de façon similaire, on¹ obtient¹

Var
N
   Ox
U ¹¹
È 
¹ U
¹ Ô
NÈ
;
Ä O NÄ O } Ä } È
; 

Û y^  

U ¹ U ¹
N ½ Ä O  N [O N Ä O } Ä }
Ô ; 
N Û b˜Û ùO     Var
N
 OÞ Var
N
 O$

Remarque. Si et sont
les variables aléatoires
   
 ý et  obéissent aussi à cette loi (exercice 16).  
deux variables aléatoires indépendantes qui obéissent à la loi normale,

Revenons maintenant au problème (8.2). Pour pouvoir estimer l’erreur du résultat numérique
, faisons les hypothèses suivantes :
H1: La valeur 
est la réalisation d’une épreuve pour une variable aléatoire ºy
. On suppose que
les ºk
soient indépendantes
 et qu’elles obéissent à la loi de Gauss–Laplace
 avec ý
comme
espérance et s
comme variance (les ý
sont inconnus, mais les s
sont supposés connus).
les '
par¸ les
 
H2: Le système surdéterminé (8.2) possède une solution X  unique
(
si l’on remplace
N
nombres ý
, c.-à-d. qu’il existe un vecteur ±Tn tel que  ý où ý  ý & ý ³ O .

Motivation de la méthode } des moindres

} carrés. Par l’hypothèse H1, la probabilité que ºy
soit
dans l’intervalle 
_
. 
šO avec 
(infinitésimalement) petit est

` N
} { 6 6 
;bý
 }
– ºy
[¾*
 '
çO æ j j 
š
4 4 j s
-'¦ ;
s

Comme les y
º
sont indépendants, la formule (8.7) implique que
` N '
} { ³ 6 6 
;bý

 }
– y
º [

*
¾  .

  _

 :  6 
 
 
 
  ª O æ j '
- ¦ ; j 
(8.11)

8]V 4 j s
4 s

 ãj -'¦ ;
6
³ '
;bý

*
 
ã j -'¦ ;
6
³ '
; ]V 
  
 
4
8]V s
4
8]V s

Selon une idée de Gauss (1812), la “meilleure” réponse 
pour les
(inconnus) est celle pour
laquelle la probabilité (8.11) est maximale (“maximum likelihood”). Alors, on calcule #h

de façon à ce que ³  



; j# Q ¤ÅˆÆ  (8.12)

8]V s
$]V s

Si l’on remplace '
J s
par 
et 
 J s
par 
 , la condition (8.12) est équivalente à (6.2). Par la suite,
nous supposerons que cette normalisation soit déjà effectuée (donc, s
 6 pour :  6 &# 5 ).


Estimation de l’erreur
N ( ¸ ( O UV ( ¸ 
La solution de (8.12) est donnée par F
N ( ¸ ( O UV ( ¸
. La solution théorique satisfait 
ý . Alors,
³

;  N ( ¸ ( O UV ( ¸ N 
;bý O ou 
; 

$]V
Â
 N  
;bý PO

 N( ¸ ( UV ( ¸
où
 est l’élément (:  § ) de la matrice O . L’idée est de considérer la valeur 
comme
la réalisation d’une variable aléatoire
définie par
³
 ³

Þ
$]V
Â
ºy
ou  
;
[
]V
Â
 N ºy
;˜ý +O (8.13)
Systèmes d’Equations Linéaires 107

Théorème 8.3 Soient ºŒ$#º ³ des variables aléatoires indépendantes avec ý
comme espérance
et s
[ 6 comme variance. Alors, la variable aléatoire
, définie par (8.13), satisfait 
ñ
N

Ox 
et Var
N

çO9 ™_

(8.14)
N( ¸ ( O V
U
où ™_

est le : ème élément de la diagonale de .

Remarque. Les autres éléments de

N ( ¸ ( O UV
sont les covariances de
et  .  
Démonstration. La formule N ºy
çOx (8.9) donne ñ
N
O)
. Pour calculer la variance de
, nous
6 et la formule (8.10). Ceci donne avec


N 2d&&&2d &2d &+&2"O ¸
6
 
utilisons le fait que Var
que 
s 
³
$]V
 
  ¬ ¸
N ( ¸ ( O UV ( ¸ ¬   ¸ N ( ¸ ( O UV ( ¸ ( N ( ¸ ( O UV



¸ N ( ¸ ( O UV
[B™_

_

Exemple 8.4 Pour l’expérience sur la thermo-électricité (voir le paragraphe IV.6), on a supposé
que les mesures 
ont été faites f
avec une précision correspondant à s
’ä2d†2 6 . Pour le système
surdéterminé (on écrit   h  ! '
 pour \+Y et 
pour Ì
)

6 j&[ õ
#j % õ
#j f  

 :  6 # 46
s
s
s
s

N ( ¸ ( O UV
la matrice devient
"!
"! "#
U‡ U U
2  G › œ’j 6 2
" ; d 2  #6 G ŠIj 6 2 2d 6"6#G j 6 2 Í

"#
N( ¸ ( O UV U U U
 ; 2d 6G ŠIj 6 2 U 2d Ÿ œ › j 6 2 Í ; d 2  Ÿ6#G j 6 2 (8.15)
U U V
2" 6"6G j 6 2 Í ; 2d Ÿ6#G j 6 2 2d Ÿ6#G j 6 2

et on obtient
$„*2d†œ"2’j 6 2
U
 &%
9B2" 4 ‹Ij 6 2
U f
 (' xB2d 4Ÿ j 6 2
"!
U

)
s s s

Ceci implique qu’avec une probabilité de Š › , la solution exacte (si elle existe) satisfait
U‡ U‡
ð ; 2dˆ‹‹"œIe2d†2 6#4  kB2d†2 G › 4 e2dˆ22"2"œd YyB2d › Š"‹’j 6 2 e2dˆ2 › ž’j 6 2 

Test de confiance du modèle

N 
Etudions encore si les données
_ ç
O sont compatibles avec la loi (8.1). Ceci revient à justifier
l’hypothèse H2. X ( (
En utilisant la décomposition © de la matrice , le problème surdéterminé S1 se trans-
forme en (voir (7.1))
X v Yv Yv
²  ¸ +
2 Y vv où Y vv © (8.16)
¬ ¬ 
La grandeur de Y v v  est une mesure de la qualité du résultat numérique. Théoriquement, si l’on a

ý à la place de  et à la place de , cette valeur est nulle.
de la matrice © par
 . Alors, les éléments ³ du vecteur Y> © ¸  sont *
* *
Notons les éléments
³ N
donnés par Yh
 ]V š
š  et ceux du vecteur Y'v v satisfont aussi Y'
 $]V š
  ;uý PO . Il est alors
naturel de considérer les variables aléatoires
³
ã9

$]V * š
N ºy ;bý PO$ :  5  6  ª  (8.17)

³
8]""^V ã

Le but est d’étudier la fonction de densité de .
108 Systèmes d’Equations Linéaires

Lemme
N
8.5 Soient º$&&º ³ des variables aléatoires indépendantes satisfaisant la loi normale
Ò ý
_ 6 O . Alors, les variables aléatoires ã9"^V$&&&ã ³ , définies par (8.17), sont indépendantes et
satisfont aussi la loi normale avec
N ã9
O9B2d N ã
Ox
ñ Var 9 6  (8.18)
Démonstration. Pour voir que les 9 ã
sont indépendants, calculons la probabilité ` N 
– ã9
g¾

N ¸ N ¸
et par ã et º les vecteurs ã’$&&ã ³ O et º#º ³ O . Alors, on a
+

_ :  5  6 & ª O . Notons par l’ensemble  ò Ä ±Tn X ³ • 
– Ä
[¾ 
_ :  5  6 & & ª
ó
+
` N 

– ã9
¾*
š :  5  6 # ª O9
` N ã*±
+ O9 ` N
©
¸ Nº
;Rý % O ± O +
 ` N º ;bý ± © N O$O (a) + , 6 6
³
Ä
 } Ä . } Ä
.- /
 ³
@ Næ
0
O³ -'¦ ; (8.19)
? 4 4 8
]V
³ ³ ¼ È

 /
È
 } È } } È
6 6 V È ³ 6
-

(b) 
š
Næ O³ -'¦ ; æ -'¦ ;
4 4
8]V
8]""^V 4 4
des ºy
et (b) découle de la transformation
+
L’identité
Ä  È
(a) est une conséquence
Ä
 
de l’indépendance
È
 Ä
©• , car ,.-'/„© ó  6 et

(la matrice © est orthogonale). En utilisant
x÷ò ±
nX ³ A
– Ä
¾*
, on déduit de la même manière que
0 ¼ È

` N 

– ã9
[¾B
šO9
4
` N ãB±
+
Ox*&.
 æ
6
4
-'¦ ;
} È

 (8.20)

Une comparaison de (8.19) avec (8.20) démontre l’indépendance

N 2d O de ã9^V$#ã ³ (voir la définition
9
ã

(8.7)). Le fait que les satisfont la loi normale Ò 6 est une conséquence de (8.20).

Théorème 8.6 (Pearson)

N
à la loi normale Ò 2d 6 O . Alors, la fonction de densité de la variable aléatoire

Soient $&&  des variables aléatoires indépendantes qui obéissent 
  


  è  

(8.21)
est donnée par (voir figure IV.5)

1
6 # ·_UVN U ® ·_
&·_ j N O
j& Ô  [O%
j (8.22)
4 5 J 4
N
pour ¨ 2 et par Ô  [Ox*2 pour – 2 (“loi de

à 5 degrés de liberté”). L’espérance de cette 2
variable aléatoire vaut 5 et sa variance 4&5 .

5  G
.2
*‹
5
 Š › %
.1 5 6 ‹

0 10 20 30 40
F IG . IV.5: Fonction de densité (8.22)

Démonstration. Considérons d’abord le cas 5  6 . Pour 2 – ¾C , on a

` N
– ¾*POx ` N æ
¼
–

x¾ O ` N ; æ Á
¼
¨*;
æ
+ O
}3 
æ
%  4 
 4 æ
6 j
U ® ·_ }
) æ
6 j
U ·_
j æ 3   

N
ce qui démontre (8.22) pour 5  6 car 6 J 4 Où
4
æ
.
4
1
Systèmes d’Equations Linéaires 109

   
Pour le cas général, nous procédons par
 récurrence.
 Nous utilisons le résultat du Lemme
 8.2
qui affirme que  la fonction de densité de qç "^V est la convolution de celle de q 

avec celle de ^V . Le calcul
®
 4  4 
1 51
N Œ N 6 N UV·_ U ® U Ë@£·_ &·_UV U ·_ }
Ô  = Ô $O [Ox æ N J Oj #·_ j N J O | ; O ?
4 j 6 4 4 5 4

 æ
1 N J 
U ® ·_
#·_ j N 51 |
®
N
;
  
O
UV·_ &·_UV }
4 j 6 4 O j 4 5 J 4 O
"^V £@ ·_UV U ® _· 
 æ
4 j1
?
N J Oj #·_ j N J O
6 4 4 51
5 4
|
N
6x; 76 6 O
·_ &·_UV }
6  N
Ô ^V [O

nous permet de conclure.

Pour les variables aléatoires ã9
de (8.17), ce théorème montre que

³ 
ã
(8.23)

8]""^V
N
est une variable aléatoire ayant comme fonction de densité Ô ³ U" [O (on rappelle qu’après normal-
isation, on a s
 6 pour les variables aléatoires ºy
).
Appliquons
¬ Y v v ¬   ce résultat à l’exemple du paragraphe IV.6 (voir la formulation (8.12)). Dans ce
› 
cas, on a 4 4 et ªä;b5  6 ‹ degrés de liberté. La figure IV.5 montre que cette valeur de
¬ Y v v ¬ 
est suffisamment petite pour être probable.
Si l’on avait travaillé avec le modèle plus simple

Ì  g õ (8.24)
¬ ¬ › œd 
(à la place de (6.4)) on aurait trouvé Y v v   4 G et ªF;‘5  6 Š . Cette valeur est trop grande
pour être probable. La conclusion est que, pour les données du tableau IV.2, la loi (8.24) est à
rejeter sur la base de ces mesures.

IV.9 Exercices

exige environ 8
f
1. Pour calculer l’inverse d’une matrice dont la dimension
opérations. Donner cet algorithme.
8 est très grande, il existe un algorithme qui

2. Supposons que la décomposition 9:

de la matrice est à disposition. Pour ; <>=@?A=CB des vecteurs

D ;FEG<H? IKJGL B
donnés, trouver un algorithme efficace pour résoudre le système
¸

; B ; < UV UV

? J
qui utilise uniquement la résolution des systèmes et . Cet algorithme est connu sous la
formule de Sherman - Morrison - Woodbury. ¸
Indication. Calculer d’abord une formule pour .

3. Considérons le problème de calculer le produit scalaire


M J =@NPO L J N RQ
]V

Quelle est sa condition?

110 Systèmes d’Equations Linéaires

; 8TSVU
W; W L W; W W ; W L W ; W XQ
4. Soit une matrice à coefficients réels. Montrer que
¸ ¹  ¸ 
et

\ ^L ] _ `badceY5Z af`b_gceY5[
5. Considérons une matrice-bande avec une largeur inférieure et une largeur supérieure
dire,
H
(c’est-à- Y"Z 9 : Y5[
5Y9 Z : L Y [ Lih
si et si ). Montrer que les matrices et de la décomposition
avec et sans la recherche de pivot ont aussi une structure de bande. Pour le cas tridiagonal,
, donner les largeurs des bandes apparaissants dans les décompositions et estimer le coût
en opérations des algorithmes.

6. Pour résoudre le système linéaire


J L B = _ Lkh = QCQCQ =@8
j

UV
 
(9.1)
]V

(matrice du type Vandermonde), dériver un algorithme qui nécessite seulement l D8 I


opérations.
Indications.
(a) Le système (9.1) est équivalent à

Y D j IKJ L B D Y I YVmT8n` h =

  pour deg
]V

B D Y I L KD q Ipo B Y D _ I L Ho h _ q DKq I DKq I DKq I DKq I DKq I


$]V   et

$]V 
UV

= ` j = ` j ` j = ` j ` j ` j = QCQCQ
où
Y
.
  f
(b) Choisir pout les éléments de la base

h ]
; L ` r h rh
7. (a) Pour la matrice

W; W W; W W; W  et ¹
calculer , .
 W; W m W; W
B LFsCt"u Dv ` Dv ` hwI I
(b) Démontrer que pour des matrices symétriques nous avons toujours .
·f ·f
8. Les valeurs de la suite Ù

B LFLFsCsCt"t"uu DDvv I` xpsCDt"v u ` DyDvhwI hwI I = II

peuvent être calculées par les formules:
· f · f
Ù
Ù
·f
B LFsCt"u DDzwv hwIx Dv ` Dv{Dv = hwII Dv hwI II Q
·f

Ù B ‡·f
` ·f ‡·f
E ` E `
Lih Q ]wh|rwr~}w}ww}|r z r z h€]ww]€‚wƒ„}w‚€ z wh€‚w‚w}€]w] v i
L w
h w
] w
] w
] w
] ]
B
Calculer à l’aide d’une calculatrice la valeur pour
||||| ) et pour
v (le résultat est
grand, quelle formule est préférable pour

J ` z Y J `† Li]
un calcul en virgule flottante?

9. Les racines du polynôme peuvent être calculées par

‡ L YˆE Y E = ‡ L Y‰` Y E† Q
  

Montrer que pour

l’aide de la relation
YG‡ c ‡ ] L `f Šc ]

(grand) et (très petit) cet algorithme est numériquement instable. A
, trouver un algorithme qui est numériquement stable.
Systèmes d’Equations Linéaires 111

10. Soient donnés  J = J = CQ QCQ = J  . Une estimation de la variance peut être calculée par chacune des deux

‹ L 8n` h h J `Œ8P ‹ L 8n` h h J G` 

formules suivantes:
  
   




]V
8]V

L  J
8 L z J LiŽw]w]wh J LiŽw]w]wŽ
où 
est l’espérance. Quelle formule est la plus stable?

8]V

r
a) Appliquer les deux algorithmes à l’exemple , ,  et simuler un calcul en

8 L z J J
virgule flottante avec chiffres.
b) Etudier l’influence des erreurs d’arrondi pour les deux algorithmes, si mais que et

sont arbitraires.
¸
; L 9(9
L;
11. a) Calculer la décomposition de Cholesky
h
pour la matrice de Hilbert
LiŽ = } = ‚ = h z = hw Q
_ EGag` h =

Ëc $]V cµ´µ´µ´ c 8
 L‘ zwDv ` h“v ’I”„D aX’ D ` hwI”„v ’ D aghw` I” hwI” Q
b) Comparer le résultat numérique avec les valeurs exactes

Ù
aX` agE ` (9.2)

9 ;F` 9 9
Combien de chiffres sont exacts?
¸

;F`†9(9 9
c) Si dénote le résultat numérique, calculer le résidu .
¸

L; D \ I D\ I
Calculer aussi le résidu pour la matrice , donnée par (9.2).
Ù!@
12. Pour une matrice
H notons par
?  les matrices des étapes intermédiaires de l’élimination

•—– Œt ˜ \ ˜ m z ’ •—– t ˜ \ ˜ Q
de Gauss. Montrer qu’avec une recherche de pivot partielle, on a
Ù!@ UV


? 
8c 


Ëc  c Ù (9.3)

Pour la matrice suivante, on a égalité dans la formule (9.3):

hh h hh
; L `` h ` h h h Q
` h ` h ` h QCQCQ ` h h
.. .. .. ..
. . . .

; U 8 I œL tUš™FJ8 x
› D ; ž•—€ž– ˜˜ ; ˜˜ ž •—/ž Ÿ ˜˜ ; N ˜˜ Q
13. Soit une matrice à lignes et colonnes ( ). On définit pour les matrices non-carrées,

® ]V Ã ]V

Pour la norme Euclidienne, montrer que 

¸ › D ; ; I L D › D ; II 

‡ ™ ‡ ™ CQ QCQ ™ ‡ ™ ]
¸
; ; .
 

•¡–% t W ; J W L ‡ = •¡% Ÿ W ; J W L ‡ Q
Indication. Transformer la matrice symétrique sous forme diagonale
et montrer que  
  
­®­ ]V ­®­ ]V

¢ 
^¢ D £ ¥¤“I ¦ ] Q ‚w‚w‚ww}w ] Q ‚w‚w‚w‚w‚ z ] Q ‚w‚whw‚w] ƒ z  ] Q ‚w‚whw‚w h z } ] Q ‚w‚wz  ] z Ž z
]
14. Voici quelques valeurs pour la densité de l’eau en fonction de sa température ø .

ø ö
ø
112 Systèmes d’Equations Linéaires

(a) Approcher ces valeurs par un polynôme de degré 2 (méthode des moindres carrés).
(b) Pour quelle valeur de ø , la densité est-elle maximale et quelle est cette valeur maximale?
J§L xw
¨ D ILih `†¢ D I
Indication. Si vous préférez calculer avec des nombres plus petits, faites la transformation
ø ø .
ø ,

;: 8
15. Soit une matrice inversible de dimension . Montrer que la décomposition
nale et triangulaire supérieure) est unique, si l’on suppose que

ª c ] f© _ : Lih = QCQC© Q =@8 I
(où est orthogo-

D = ‹
pour .

­ D  = ‹« I ¬
16. Soient
 
et
®¯«°E†± ]
²® c «°E³¬
deux variables aléatoires indépendantes obéissant à la loi normale
respectivement. Montrer que (pour ) et
­
et
obéissent aussi à cette

¨ D J>I
loi.

17. Soit « ´
une variable aléatoire qui obéit à la loi

avec 8 degrés de liberté (c.-à-d.,  de (8.22)

µ D« I L 8 D « I(L z 8 Q
est sa fonction de densité). Montrer que

et Var

W j W L W ; J `†B W
18. Effectuer une étude complète de l’erreur du modèle trouvé à l’exercice 14. Pour cela, trouver les

écarts types des coefficients du polynôme et effectuer un test de confiance du modèle.
Indication. v v
¤ L D; ; I
.
¸ UV
19. Les éléments de la diagonale de

;
jouent un rôle important pour l’étude de l’erreur de
la méthode des moindres carrés. Supposons que nous avons à disposition la décomposition de la ©f:
¤ L D: : I
matrice .
¸ UV
(a) Démontrer que .
¤ 8 xw}
f
8
;
(b) Trouver un algorithme pour calculer la diagonale de en
colonnes de ; 1 opération = 1 multiplication + 1 addition).
opérations ( = nombre de