0% ont trouvé ce document utile (0 vote)

65 vues132 pages

MNE Cours

Transféré par

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

65 vues132 pages

MNE Cours

Transféré par

Hope2904

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

METHODES NUMERIQUES

Manfred GILLI

Dpartement dconomtrie e e e Universit de Gen`ve e e

Version : 25 mars 2006

Table des mati`res e

1 Arithmtique en prcision nie e e 1.1 Reprsentation en virgule ottante e 1.1.1 Normalisation . . . . . . . . 1.1.2 Underow et Overow . . . 1.2 Prcision, erreurs darrondi . . . . . e 1.2.1 Prcision machine . . . . . . e 1.2.2 Digits signicatifs . . . . . . 1.3 Mesures de lerreur . . . . . . . . . 1.3.1 Erreur relative de oat(x) . 2 Instabilit numrique et condition e e 2.1 Condition dun probl`me . . . . . e 2.2 Condition dune matrice . . . . . 2.3 Remarques . . . . . . . . . . . . . 3 Complexit des algorithmes e 3.1 Taille du probl`me . . . . . . . e 3.2 Crit`res de comparaison . . . . e 3.3 La fonction O() . . . . . . . . . 3.4 Oprations lmentaires (ops) e ee 3.5 Classication des algorithmes . 1 2 2 4 5 6 8 8 9

. . . . . . . .

dun probl`me e 11 . . . . . . . . . . . . . . . . . . . . 11 . . . . . . . . . . . . . . . . . . . . 12 . . . . . . . . . . . . . . . . . . . . 13 17 17 17 19 20 20 23

. . . . .

4 Resolution numrique de syst`mes linaires e e e

5 Syst`mes triangulaires e 27 5.1 Forward substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.2 Back substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.3 Proprits des matrices triangulaires unitaires . . . . . . . . . . . . . 28 ee 6 Factorisation LU 6.1 Formalisation de llimination de Gauss e 6.2 Matrices de permutation . . . . . . . . 6.3 Pivotage partiel . . . . . . . . . . . . . 6.4 Considrations pratiques . . . . . . . . e 6.5 Elimination de Gauss-Jordan . . . . . 3 31 32 38 41 43 45

. . . . .

7 Matrices particuli`res e 7.1 Factorisation LDM . . . . . . . . . . 7.2 Factorisation LDL . . . . . . . . . . . 7.3 Matrices symtriques dnies positives e e 7.4 Factorisation de Cholesky . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . .

47 47 50 51 52 55 55 56 56 58 59 59 60 63 67 68 69 71 73 74 75 79

8 Matrices creuses 8.1 Matrices par bande . . . . . . . . . . . . . . . . . . . . 8.2 Matrices creuses irrguli`res . . . . . . . . . . . . . . . e e 8.2.1 Reprsentation de matrices creuses . . . . . . . e 8.2.2 Conversion entre reprsentation pleine et creuse e 8.2.3 Initialisation de matrices creuses . . . . . . . . . 8.2.4 Principales fonctions pour matrices creuses . . . 8.3 Proprits structurelles des matrices creuses . . . . . . ee 8.3.1 Exemples de matrices creuses . . . . . . . . . . 9 Mthodes itratives stationnaires e e 9.1 Mthodes de Jacobi et de Gauss-Seidel e 9.2 Interprtation gomtrique . . . . . . . e e e 9.3 Convergence de Jacobi et Gauss-Seidel 9.4 Mthodes de relaxation (SOR) . . . . . e 9.5 Structure des algorithmes itratifs . . . e 9.6 Mthodes itratives par bloc . . . . . . e e 10 Mthodes itratives non stationnaires e e 11 Equations non-linaires e 11.1 Equation ` une variable . . . . . . . . . . . . . . a 11.1.1 Recoupement par intervalles (bracketing) . 11.1.2 Itrations de point xe . . . . . . . . . . . e 11.1.3 Bisection . . . . . . . . . . . . . . . . . . . 11.1.4 Mthode de Newton . . . . . . . . . . . . e 11.2 Syst`mes dquations . . . . . . . . . . . . . . . . e e 11.2.1 Mthodes itratives . . . . . . . . . . . . . e e 11.2.2 Formalisation de la mthode du point xe e 11.2.3 Mthode de Newton . . . . . . . . . . . . e 11.2.4 Quasi-Newton . . . . . . . . . . . . . . . . 11.2.5 Newton amorti (damped ) . . . . . . . . . . 11.2.6 Solution par minimisation . . . . . . . . . 11.3 Tableau synoptique des mthodes . . . . . . . . . e 12 Dcompositions orthogonales e 12.1 La factorisation QR . . . . . . . . . 12.1.1 Matrices de Givens . . . . . 12.1.2 Rexion de Householder . . e 12.1.3 Algorithme de Householder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

81 82 86 88 89 91 92 92 92 95 99 100 101 101 105 105 106 106 108

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Businger - Golub

12.2 Dcomposition singuli`re SVD . . . . . . . . . . . . . . . . . . . . . e e 12.2.1 Approximation dune matrice par une somme de matrices de rang un . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.2.2 Calcul numrique du rang dune matrice . . . . . . . . . . . e 12.2.3 Interprtation des valeurs singuli`res . . . . . . . . . . . . . e e 12.2.4 Complexit . . . . . . . . . . . . . . . . . . . . . . . . . . . e 12.2.5 Calcul du pseudoinverse . . . . . . . . . . . . . . . . . . . . 13 Moindres carrs e 13.1 Moindres carrs linaires . . . . . . . . . . . . . e e 13.1.1 Mthode des quations normales . . . . . e e 13.1.2 Solution ` partir de la factorisation QR . a 13.1.3 Dcomposition ` valeurs singuli`res SVD e a e 13.1.4 Comparaison des mthodes . . . . . . . e 13.2 Moindres carrs non-linaires . . . . . . . . . . e e 13.2.1 Mthode de Gauss-Newton . . . . . . . . e 13.2.2 Mthode de Levenberg-Marquardt . . . . e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 108 . . . . . . . . . . . . . 109 109 109 110 111 113 114 114 117 118 119 120 122 123

Chapitre 1 Arithmtique en prcision nie e e

Sources dapproximation dans le calcul numrique e Les deux sources derreur qui interviennent systmatiquement dans le calcul numrique e e sont : e Les erreurs de troncature ou de discrtisation qui proviennent de simplications du mod`le mathmatique comme par exemple le remplacement dune drive par e e e e une dirence nie, le dveloppement en srie de Taylor limit, etc. e e e e Les erreurs darrondi qui proviennent du fait quil nest pas possible de reprsenter e (tous) les rels exactement dans un ordinateur. e Dans ce chapitre on expliquera dabord quelle est lapproximation utilise pour e reprsenter les nombres rels dans un ordinateur et ensuite on discutera les consquences e e e de cette reprsentation inexacte. e Rappelons que dans un ordinateur le support pour toute information est un mot constitu par une squence de bits (en gnral 32), un bit prenant la valeur 0 ou 1. e e e e
231 230 23 22 21 20

0 0 0 0 1 0 1 0

Cette squence de bits est alors interprte comme la reprsentation en base 2 dun e ee e nombre entier. Ainsi les nombres entiers peuvent tre reprsents exactement dans un ordinateur. e e e Si tous les calculs peuvent se faire en nombres entiers, on parle darithmtique en e nombres entiers. Ces calculs sont alors exacts.
1 e e e Une division 3 ncessite un nombre inni de digits pour reprsenter le rsultat exactement. Dans la pratique, on recourra ` larithmtique en virgule ottante pour a e reprsenter une approximation du rsultat. Ainsi appara le probl`me des erreurs e e t e darrondis.

Remarque : Il existe une catgorie de logiciels (e.g. Maple) qui permettent de faire e des calculs symboliques en traitant des expressions algbriques. Ces mmes logiciels e e 1

eectuent des calculs, impliquant des nombres rels, avec une prcision quelconque, e e limite uniquement par la performance de lordinateur. Ceci ne permet videmment e e pas de contourner les probl`mes dans la plupart des calculs numriques. e e

1.1

Reprsentation en virgule ottante e

Dans un ordinateur, une variable relle x = 0 est reprsente de la faon suivante : e e e c x = n be o` n est la mantisse (ou partie fractionnelle), b la base et e est lexposant. Dans un u ordinateur, la base est toujours gale ` 2. e a Exemple 1.1
12.153 = .75956 24 ou en base 10 on a .12153 102

Pratiquement on partitionne le mot en deux parties, lune contenant e et lautre contenant n. (Le premier bit ` gauche indique le signe). a

Quelle que soit la taille du mot on ne disposera donc que dun nombre limit de bits e pour reprsenter les entiers n et e. e Soit t le nombre de bits disponibles pour coder la mantisse. En base 2, on peut alors coder les entiers allant de 0 ` t1 2i = 2t 1. a i=0 Exemple 1.2
t=3 0 0 0 =0 1 1 1 = 22 + 21 + 20 = 23 1 = 7

1.1.1

Normalisation

An de maximiser le nombre de digits signicatifs, on normalise la mantisse, cest`-dire on limine les digits nuls ` gauche. a e a Exemple 1.3
nombre comme x = 0 0 1 2 3 105 ou 1 2 3 4 5 107 Soit x = 0.00123456. En base 10 et avec 5 digits on peut reprsenter ce e

La deuxi`me reprsentation exploite au maximum le nombre de digits disponibles. e e

Ainsi, pour un mot en base 2, on fait varier la mantisse n dans lintervalle

1 t 1 t

1 0 0
2t1

1 1 1
2t 1

En rajoutant 1 ` droite (ce qui change le signe en <) et en multipliant la mantisse a par 2t , la valeur normalise de la mantisse n variera dans lintervalle e
1 2

n 2t < 1 .

Lensemble des nombres en virgule ottante F que lon peut ainsi reprsenter constie tue un sous-ensemble de R. Si lon reprsente la mantisse avec un mot de t bits, les e lments f F sont dnis par ee e f = 1 d2 dt 2t 2e = n 2et

avec di = 0 ou 1 pour i = 2, . . . , t et d1 vaut toujours 1 a cause de la normalisation. Si lexposant peut prendre les valeurs enti`res de lintervalle [L, U] on a que e m |f | M Exemple 1.4 avec m = 2L1 et M = 2U (1 2t ) .

Construction de lensemble F pour t = 3 et e [1, 2]. f = n 2et n

1 1 1 1 0 0 1 1
2 0 =2 =4 2 0 1 =2 +2 =5 2 1 0 =2 +2 =6 2 1 0 1 =2 +2 +2 =7

2et

e = 1 e = 0 e = 1 e = 2 1/16 1/8 1/4 1/2 1/4 5/16 3/8 7/16 1/2 5/8 3/4 7/8 1 5/4 3/2 7/4 2 5/2 3 7/2

En reportant ces nombres sur une droite on observe que les nombres en virgule ottante ne sont pas galement espacs. e e
.................................................................................................................... ....................................................................................................................

Pour dcrire lensemble des nombres rels qui trouvent une reprsentation dans F e e e on dnit lensemble e G = {x R tel que m |x| M} {0}

et loprateur oat : G F qui fait correspondre ` un lment de G un lment e a ee ee de F comme plus proche f F ` x qui satisfait |f | |x| (chopping) a oat (x) = plus proche f F ` x (perfect rounding). a Exemple 1.5
Soit lensemble F dni pour t = 3 et e [1, 2] dans lexemple 1.4, e alors avec chopping on a oat (3.4) = 3 et oat (0.94) = 0.875 et avec perfect rounding le rsultat est oat (3.4) = 3.5 et oat (0.94) = 1. e

1.1.2

Underow et Overow

Soit op lune des 4 oprations arithmtiques + et a, b deux nombres rels ; e e e alors, si |a op b| G, on est dans une situation derreur (overow si |a op b| > M et underow si |a op b| < m). Matlab utilise des mots en double prcision de 64 bits avec t = 52 le nombre de bits e de la mantisse et e [1023, 1024] les valeurs possibles pour lexposant. On a alors m 1.11 10308 et M 1.67 10308 . La fonction Matlab realmin et realmax produit m respectivement M. Les 11 bits rservs ` lexposant permettent de reprsenter des entiers (positifs) de e e a e 11 0 ` 2 1 = 2047. En soustrayant ` cet entier la constante 1023, appele oset, et a a e s1 qui est dnie comme o = 2 1 avec s le nombre de bits rservs pour lexposant e e e on peut reprsenter les entiers dans lintervalle e ce qui correspond ` 0 1023 = 1023 et 2047 1023 = 1024 comme indiqu plus a e haut. Dans le standard IEEE1 les situations doverow et underow ne provoquent pas larrt des calculs. Un overow produit le symbole Inf qui se propage dans la suite e des calculs et un underow peut produire le rsultat zro. e e Ainsi une situation doverow est fatale pour la poursuite des calculs. Dans la pratique on peut cependant souvent transformer une situation doverow en une situation dunderow sans consquences pour la prcision des calculs. e e Exemple 1.6
Evaluons avec Matlab lexpression c= a2 + b2

(0 o) e (2s 1 o)

avec a = 10200 et b = 1. Le rsultat attendu est c = 10200 mais avec Matlab a2 produit un e overow do` c = Inf. On peut viter ce probl`me en procdant a une normalisation u e e e ` c=s
1

a s

b s

Institute of Electrical and Electronics Engineers. Organisation professionnelle avec plus de 100 000 membres qui chapeaute 35 socits techniques dont la Computer Society. ee

avec s = max{|a|, |b|} = 10200 , do` u c = 10

200

b 10200

= 10200

car ( 101 )2 produit un underow avec la valeur zro. Ceci ne g`ne cependant pas la e e 200 prcision de calculs tant donne que la quantit 10400 est insigniante lorsque lon ade e e e ditionne a 1. `

1.2

Prcision, erreurs darrondi e

Avant de prciser la dnition de prcision dun syst`me en virgule ottante, illuse e e e trons avec un exemple les probl`mes qui peuvent surgir lorsque les calculs se font en e prcision nie. e

Exemple 1.7

Rsolution numrique dune quation du deuxi`me degr ax2 + bx + c = e e e e e 0. Les solutions analytiques sont : x1 = b b2 4ac 2a x2 = b + b2 4ac . 2a

Algorithme 1 Transcription de la solution analytique.

1: = b2 4ac 2: x1 = (b )/(2a) 3: x2 = (b + )/(2a)

Pour a = 1, c = 2 et une arithmtique a 5 digits signicatifs lalgorithme 1 est appliqu e ` e ci-apr`s pour des valeurs direntes de b. e e

b oat () oat (x2 ) x2 5.2123 4.378135 4.3781 0.41708 0.4170822 121.23 121.197 121.20 0.01500 0.0164998 1212.3 1212.2967 1212.3 0 0.001649757 Attention : oat (x2 ) = (b+ oat ())/(2a)

erreur relative % 0.004 9.1 Catastrophe ! !

Ci-apr`s on propose un autre algorithme qui exploite le thor`me de Vi`te et qui vite le e e e e e calcul de la dirence de deux nombres qui est a lorigine de la perte de prcision si ces e ` e nombres sont tr`s proches. e

Algorithme 2 Exploite le thor`me de Vi`te x1 x2 = c/a . e e e

1: 2: 3: 4: 5: 6: 7: = b2 4ac if b < 0 then x1 = (b + )/(2a) else x1 = (b )/(2a) end if x2 = c/(a x1 )

b 1212.3

oat () 1212.3

oat (x1 ) 1212.3

oat (x2 ) 0.0016498

x2 0.001649757

Pour caractriser la prcision dun syst`me de reprsentation des nombres rels en e e e e e virgule ottante on prsente souvent les mesures suivantes : e Prcision machine (eps) e Digits signicatifs Erreur darrondi (u) Erreur relative ()

Bien que ces nombres di`rent lg`rement entre eux et dpendent du schma dare e e e e rondi (chopping ou perfect rounding), tous donnent une mesure qui caractrise la e granularit du syst`me en virgule ottante. e e

1.2.1

Prcision machine e

La prcision machine dune arithmtique en virgule ottante est dnie comme le e e e plus petit nombre positif eps tel que oat (1 + eps) > 1 . Regardons quelle est la valeur de eps pour une machine avec des mots dont la mantisse comporte t bits et lorsquon arrondit selon la technique dite chopping. Reprsentons le nombre 1 et le plus proche nombre suivant : e
1 t t 1 t

1 0 0 2
2t1

2 =1

1 0 1 2t 21 = 1 + 21t
2t1 +1

La distance qui spare le nombre 1 du plus proche nombre suivant est gale ` 21t . e e a Ainsi pour une mantisse comportant t bits on a eps = 21t . Vrions ce rsultat sur notre mantisse avec t = 3 introduite avant : e e

1 eps 2 ....................................................................................................................................................................... ......................................................................................................................................................................

7 16 5 8 6 8 7 8 5 4

1 2

On voit aussi, que dans une situation o` on utilise perfect rounding, lon obtient u 1 eps = 2 21t = 2t . Avec Matlab, qui utilise t = 52 bits pour la mantisse et perfect rounding, la prcision e machine est alors eps = 252 2.22 1016 . Exemple 1.8
Montrons, que dans certains cas, lapproximation que fournit le rsultat e dun calcul en virgule ottante peut tre surprenant. Considrons lexpression e e
k=1

1 = . k

En calculant cette somme avec un ordinateur on obtient videmment un rsultat qui est e e ni. Ce qui peut cependant surprendre est que cette somme sera certainement infrieure e a 100. Le probl`me ne vient pas dun underow de 1/k ou dun overow de la somme ` e e partielle n 1/k mais du fait que pour n qui vrie k=1 1/n
n1 k=1

< eps

1 k

la somme reste constante. Pour le montrer considrons quil existe n pour lequel e oat oat
n1 1 k=1 k

1 n
1 k=1 k

= =

n1 1 k=1 k

+ 1/n n1

n1 1 On peut facilement exprimenter que k=1 k < 100 pour des valeurs de n qui peuvent e tre envisages pour le calcul pratique et comme eps 2 1016 on tablit que e e e

1/n = 2 1016 100 do` lon trouve que n est dordre 1014 . u

1014 2

Considrant un ordinateur dont la performance est de 1 Gops (109 ops) et quil faut trois e oprations lmentaires par tape, la somme cessera de cro apr`s (3 1014 )/(2 109 ) = e ee e tre e 1.5 105 secondes, cest-`-dire apr`s quelques 42 heures de calcul sans avoir dpass la a e e e valeur de 100.

1.2.2

Digits signicatifs

La prcision machine dnit le nombre de digits signicatifs dun nombre rel dans sa e e e 2 reprsentation en virgule ottante. Dans un mot de 32 bits , on rserve en gnral 23 e e e e bits pour la mantisse ce qui donne eps = 222 2.38107 et 1+eps = 1.0000002 38 ce qui donne 8 digits signicatifs. Dans ce cas, il est inutile de lire ou imprimer des rels de plus que 8 digits. e Pour les illustrations, on choisira souvent la base 10 pour la dnition de lensemble e F . On aura alors f = d1 d2 . . . dt 10t 10e avec d1 = 0 et 0 di 9, i = 2, . . . t. On utilise ensuite la reprsentation en point e xe dans laquelle seuls les premiers t digits sont signicatifs. Voici ` titre dexemple, a pour une base de 10 et t = 3, les digits signicatifs : 2.37 , 139 00, 0.00 293 7.

1.3

Mesures de lerreur

On peut envisager plusieurs faons pour mesurer lerreur e entre une valeur apc proche x et une valeur exacte x. e Erreur absolue Elle est dnie comme e | x|. x Dans une situation avec x = 3 et x = 2 lerreur absolue vaut un, ce qui dans ce cas ne peut tre considr comme petit. Par contre la mme erreur absolue avec e ee e 9 9 x = 10 + 1 et x = 10 peut certainement tre considre comme relativement e ee petite par rapport ` x. a Erreur relative La remarque prcdente nous conduit ` la dnition de lerreur relative e e a e | x| x |x| qui est dni si x = 0 et pas dnie si x = 0 (dans ce cas lerreur absolue ferait e e bien laaire). Pour les exemples prcdents lerreur relative est respectivement 0.5 e e et 109 ce qui indique une petite erreur pour le deuxi`me cas. e
2

Ceci correspond a une prcision simple, qui est en gnral le dfaut pour la longueur dun mot. ` e e e e

Combinaison entre erreur absolue et erreur relative Lutilisation de lerreur relative pose des probl`mes lorsque x prends des valeurs e qui sont proches de zero. Pour pallier ` cet inconvnient on utilise souvent dans la a e pratique la mesure suivante : | x| x . |x| + 1

Cette mesure a les caractristiques de lerreur relative si |x| 1, et les caractristiques e e de lerreur absolue si |x| 1.

1.3.1

Erreur relative de oat(x)

Considrons maintenant lincrment de la mantisse d au dernier bit ` lextrme e e u a e droite. Cette valeur est appele roundo error, que lon note e u = 2t . Lerreur relative due ` larrondissement est a | oat (x) x| =. |x| Un rsultat fondamental de larithmtique en virgule ottante est alors que quelque e e soit x G, lerreur relative commise par lapproximation oat (x) est borne comme e || u lorsquon utilise perfect rounding et || 2 u pour chopping. Dmonstration 1.1 Soit une mantisse avec t bits. La plus grande erreur relative se e produit lorsquon reprsente 1/2. On consid`re x = 2k = 2t1 2t 2k+1 et le nombre e e t 1) 2t 2k . On a alors en virgule ottante qui le prc`de (2 e e
1 t 1 t

oat (x) x = et

1 1 1
2t 1

t k+1 1 0 0 2 2 2t1

2t+k

| oat (x) x| = 2t . |x|

La relation qui dnit lerreur relative peut aussi scrire e e oat (x) = x(1 + ) avec || u .

Cette relation est ` la base de toute tude des erreurs darrondi. a e Considrons maintenant les erreurs relatives associes aux oprations suivantes : e e e oat (a b) = (a b)(1 + 1 ) oat (a b) = (a b)(1 + 2 ) oat (a b) = (a b)(1 + 3 ) On peut alors montrer que |1 | u et |2 | u

alors que 3 nest pas garanti dtre petit. Cest le cas notamment, lorsque laddition e de deux nombres tr`s proches en valeur absolue (dirence) donne un rsultat tr`s e e e e petit. Ceci est connu sous le nom de catastrophic cancellation.
Soit lopration x = .123456 .123465 = .000009 = 9 106 . Si lon e choisit une reprsentation dcimale avec t = 5 et perfect rounding on a e e

Exemple 1.9
oat (x) =

105

105 =

105

3 =

Lerreur absolue est et vrie e

| oat (x) x| | 10 106 + 9 106 | 1 = = = .1111 6 | |x| | 9 10 9 | oat (x) x| = 106 | oat (x) x| u = 105

mais lerreur relative 3 est tr`s grande. e

Un fait important est que larithmtique en virgule ottante nest pas toujours e associative. Exemple 1.10 Pour le montrer, considrons une base de 10, t = 3 et lexpression e e e 103 + 1 1. On vrie alors aisment que
oat ( oat (103 + 1) 1) = 0 et oat (103 + oat (1 1)) = 103 . En eet pour laddition les exposants des deux chires doivent tre identiques ce qui e implique labandon de la normalisation pour le chire le plus petit oat (103 + 1) =
0 0 0

103 101 +

103 101 = 1 .

Cest une consquence du catastrophic cancellation. e Calcul de ea pour a > 0 (Golub and Van Loan, 1989, p. 62). Golden mean (Press et al., 1986, p. 18), rcurrence instable, etc. e

Exemple 1.11

Chapitre 2 Instabilit numrique et condition e e dun probl`me e

Dans la section prcdente, on a montr que lon peut seulement reprsenter un souse e e e ensemble ni des rels dans un ordinateur, ce qui conduit ` des erreurs darrondi. e a Comme on la vu, ces derniers peuvent, dans certains cas, provoquer linstabilit de e la solution ou mme la destruction totale des rsultats du calcul. e e Si la prcision des rsultats nest pas acceptable, il est important de distinguer les e e deux situations suivantes : Les erreurs darrondi sont considrablement grossies par la mthode de calcul. e e Dans ce cas, on parle dune mthode ou dun algorithme numriquement instable. e e Une petite perturbation des donnes peut engendrer un grand changement de la e solution. Dans ce cas, on parle de probl`me mal conditionn . e e Les exemples ` la n de la section prcdente illustrent des cas dalgorithmes numriquement a e e e instables et proposent galement les modications ` apporter an que les erreurs e a darrondi ne jouent plus le mme rle destructeur pour la prcision des calculs. e o e

2.1

Condition dun probl`me e

Pour valuer la sensibilit dun probl`me ` une perturbation des donnes, on dnit e e e a e e la condition du probl`me. Soit un probl`me qui consiste ` valuer f (x) et une pere e ae turbation des donnes x + x. On consid`re alors le rapport e e
|f (x+x)f (x)| |f (x)| | x | x

|f (x + x) f (x)| |x| |x| |f (x)| |x f (x)| |f (x)|

et pour x 0 on a

cond(f (x)) 11

la condition du probl`me f (x), ce qui correspond ` la valeur absolue de llasticit. e a e e Dans le cadre dun calcul numrique, il sagit de llasticit de la solution de f (x) e e e par rapport aux donnes. Si llasticit est grande alors le probl`me est dit mal e e e e conditionn. Pour des probl`mes pratiques il est le plus souvent tr`s dicile dvaluer e e e e cette lasticit. On peut distinguer 3 situations : e e f (x) tr`s grand et x, f (x) normales e x tr`s grand et f (x), f (x) normales e f (x) tr`s petit et x, f (x) normales e Considrons les deux matrices A et C et le vecteur b e A= 1 1 2 1 2 + d 1 C= 2 4 1 2+e b= 1 1 .

Si d est tr`s petit, le syst`me linaire Ax = b est un probl`me mal conditionn. Si e e e e e par contre on consid`re une mthode qui chercherait la solution de x en rsolvant e e e 3 le syst`me transform CAx = Cb, avec d = 2 et e tr`s petit, on a a faire ` une e e e a mthode numriquement instable. La gure 2.1 illustre ces situations. e e
2 1
.. ... .... .. .... ....... ....... ....... ...... . ...... ...... . . .... .. . .... ..... . .. ..... . . ... . ..... . .... . .... .... ... . . . .... .... . . .... .... . . ..... . ..... . .... .... . .... . ...... . ...... ...... . ...... . ...... . . ......... .. . ........ .. . . ... . .... .. ... .. . . . . . . . . . . . . . . . . . . . ............................................... . . . ................................................. . .

2 1

|d| 1 2

3 d= 2 |e| 1

Fig. 2.1 Illustration dun probl`me mal conditionn et dune mthode e e e numriquement instable. e

2.2

Condition dune matrice

Nous discuterons ici comment valuer la sensibilit dun probl`me ` une perturbation e e e a des donnes lorsquil sagit du probl`me de la rsolution dun syst`me linaire Ax = e e e e e b. Exemple 2.1
Considrons un syst`me linaire Ax = b, avec e e e A= .780 .563 .913 .659 b= .217 .254

o` lon vrie facilement que la solution exacte est x = [1 1] . On rappelle que, dans la u e pratique, on ne conna presque jamais la solution exacte. Avec Matlab, on obtient : t x = A\b = .99999999991008 .99999999987542 .

Ceci est un probl`me tr`s mal conditionn. En eet considrons la matrice e e e e E= .001 .001 .002 .001 5.0000 7.3085

et rsolvons le syst`me perturb (A + E)xE = b. On obtient alors e e e xE = (A + E)\b = .

Dnition de la condition de A Soit un syst`me linaire Ax = b et la perture e e bation E. On peut alors crire e (A E)(x + d) = b do` lon tire la perturbation d engendre pour x u e d = (A E)1 Ex . En appliquant M v M v pour une norme matricielle subordonne ` une e a norme vectorielle, on crit e d (A E)1 et, en divisant par x et en multipliant par
d x E A A E

E x , on obtient

(A E)1 A

ce qui rappelle lexpression de llasticit. Pour E 0 on a e e (A) = A1 A avec (A) la condition de la matrice A. Pour le moment, on ne discutera pas comment obtenir numriquement ces normes. e La commande cond de Matlab donne la condition dune matrice. Pour notre exemple, on a cond(A) = 2.2 106 . Si cond(A) > 1/sqrt(eps), il faut se mer pour la suite des calculs. Pour eps = 2.2 1016, on a 1/sqrt(eps) = 6.7 108. e (On perd la moiti des digits signicatifs.) e

2.3

Remarques

La condition dune matrice constitue une borne suprieure pour la mesure des e probl`mes que lon peut rencontrer lors de la rsolution dun syst`me linaire. A e e e e titre dexemple soit le syst`me linaire Ax = b avec e e A= 1010 0 0 1010

et b quelconque. La matrice A vrie cond(A) = 1020 mais la solution xi = bi /Aii e ne pose aucun probl`me numrique. e e Soit un syst`me linaire Ax = b et la solution calcule xc . On appelle alors erreur e e e rsiduelle la quantit r = Axc b. On imagine souvent que cette erreur rsiduelle e e e est un indicateur pour la prcision des calculs. Ceci nest le cas que si le probl`me e e est bien conditionn. Illustrons ce fait avec lexemple prcdent. Soit deux candidats e e e xc pour la solution et lerreur rsiduelle associe : e e xc = .999 1.001 r= .0013 .0016 et xc = .341 .087 r= 106 0

On constate que le candidat nettement moins bon produit une erreur rsiduelle e nettement infrieure. Donc, pour un probl`me mal conditionn, lerreur rsiduelle e e e e nest pas un indicateur utilisable. Montrons avec deux exemples quune matrice peut tre bien conditionne pour la e e rsolution du syst`me linaire associ et mal conditionne pour la recherche des e e e e e valeurs propres et vice versa. Exemple 2.2 Soit A = triu(ones(20)). A est mal conditionne pour le calcul des e valeurs propres. Avec la fonction eig on obtient
1 , . . . , 20 = 1 et en posant A20,1 = .001 on obtient 1 = 2.77, . . . , 20 = .57 .08i .

A est bien conditionne pour la rsolution du syst`me linaire, car cond(A) = 26.03. e e e e

Exemple 2.3

Soit la matrice A= 1 1 1 1+

A est bien conditionne pour le calcul des valeurs propres. Pour = 0 on obtient les e valeurs propres 1 = 0 et 2 = 2 . Pour = .001 on obtient les valeurs propres 1 = .0005 et 2 = 2.0005 .

A est mal conditionne pour la rsolution du syst`me linaire. Pour = 0 A est singuli`re, e e e e e avec cond(A) = . Pour = .001, on a cond(A) = 4002.

Pour conclure, on peut rappeler que si lon veut obtenir des rsultats numriquement e e prcis avec un ordinateur, il faut : e

un probl`me bien conditionn, e e un algorithme qui est numriquement stable lorsque lon excute avec une prcision e e e arithmtique nie, e un logiciel qui soit une bonne transcription de lalgorithme. Un algorithme numriquement stable ne pourra pas rsoudre un probl`me mal e e e conditionn avec plus de prcision que celle contenue dans les donnes. Cependant, e e e un algorithme numriquement instable produira de mauvaises solutions pour des e probl`mes bien conditionns. e e

Chapitre 3 Complexit des algorithmes e

Ce chapitre aborde de mani`re succincte la problmatique de la comparaison des e e algorithmes. Beaucoup de probl`mes sont de nature ` pouvoir tre rsolus par une e a e e varit dalgorithmes dirents. Il se pose alors le probl`me du choix du meilleur ee e e algorithme. La notion de complexit dun algorithme est utilise pour comparer les e e performances des algorithmes. Une comparaison approfondie est souvent une tche a ardue et nous nintroduirons que des notions tr`s gnrales. e e e

3.1

Taille du probl`me e

La notion de base en complexit des algorithmes est la dnition de la taille du e e probl`me. Il sagit naturellement de la quantit de donnes qui doivent tre traites. e e e e e Pratiquement, la taille dun probl`me est mesure par la longueur du plus petit codage e e ncessaire ` la reprsentation des donnes du probl`me (nombre dlments quil faut e a e e e ee traiter). Par la suite, on tentera dexprimer le temps ncessaire au droulement de e e lalgorithme comme une fonction de la taille du probl`me. e Exemple 3.1 Soit le probl`me qui consiste a rsoudre un syst`me linaire Ax = b. La e ` e e e
taille du probl`me est alors caractrise par lordre n de la matrice A si la matrice est e e e dense. Dans le cas dune matrice rectangulaire, la taille est dnie par le nombre de lignes e n et le nombre de colonnes m. Sil sagit danalyser un graphe, les grandeurs caractrisant la taille du probl`me sont en e e gnral le nombre de sommets n et le nombre darcs m. e e

3.2

Crit`res de comparaison e

Une comparaison dtaille de la performance de deux algorithmes est souvent tr`s e e e dicile. Beaucoup de crit`res peuvent tre utiliss. Le principal crit`re est le temps e e e e 17

dexcution de lalgorithme. An dobtenir des rsultats qui soient indpendants de e e e la performance dun ordinateur particulier (ou compilateur particulier), on mesure le temps en comptant le nombre doprations lmentaires excutes. e ee e e On appelle complexit de lalgorithme le nombre doprations ncessaires pour rsoudre e e e e un probl`me de taille n. Cette complexit est donc proportionnelle au temps dexcution e e e eectif. Exemple 3.2 Soit un algorithme qui recherche de faon squentielle un mot particulier c e
dans une liste de n mots. La complexit est alors donne par e e C(n) k1 n + k2 avec k1 0 et k2 0, deux constantes indpendantes de n qui caractrisent une installae e tion particuli`re. e

On retient deux mesures pour compter le nombre doprations lmentaires pour e ee caractriser la performance dun algorithme, le pire des cas et le cas moyen. e Nombre doprations dans le pire des cas Il sagit dune borne suprieure e e du nombre doprations ncessaires pour rsoudre un probl`me de taille n dans sa e e e e conguration la plus dfavorable pour lalgorithme tudi. e e e Exemple 3.3 Pour lalgorithme squentiel de lexemple prcdent, cest la situation ou e e e
le mot recherch se trouve en derni`re position dans la liste. On a C(n) = k1 n + k2 . e e

Nombre moyen doprations Souvent des algorithmes qui sont tr`s mauvais e e pour le pire des cas se comportent tr`s bien en moyenne1 . La dtermination du e e nombre moyen doprations ncessite souvent un analyse mathmatique sophistique e e e e et, la plupart du temps, on recourt simplement ` lexprimentation. a e Autres crit`res Dautres crit`res sont la place mmoire ncessaire ` la rsolution e e e e a e du probl`me (space complexity) et la simplicit de la mise en oeuvre dun algorithme. e e Pour des ordinateurs ` architecture parall`le, la quantit et la frquence des coma e e e munications entre processeurs sont des facteurs extrmement importants. Le temps e ncessaire pour communiquer des donnes est souvent de plusieurs ordres de grane e deur plus lev que le temps ncessaire pour eectuer des oprations lmentaires e e e e ee dans un processeur.
Lalgorithme du simplexe pour les probl`mes de programmation linaire constitue un exemple e e frappant tant donn que dans le pire des cas il ncessite 2m itrations, avec m le nombre de e e e e contraintes, et que lon observe que dans la pratique ce nombre nexc`de en gnral pas 3m/2. e e e
1

3.3

La fonction O()

Soient deux algorithmes A1 et A2 de complexit CA1 (n) = 1 n2 et CA2 (n) = 5n, e 2 alors lalgorithme A2 est plus rapide que A1 pour n > 10. En fait quels que soient les coecients de n2 et n, il existe toujours un n ` partir duquel A2 est plus rapide a que A1 .
C(n) 50
.. ... . .. . ... . . .. A1.... ............ .. .. . .. ... . . ... . . . .. ..... . ........ . . A2 . . . ...... ..... . . .... ... ...... ...... .. . . .. ... .... . . . ... . . . ..... . ...... . ....... . .. ... . . ... .. ... .. . . ... ... .... ... . . .... ..... . ... ..... . ... ... . . ... .. ... . ...... ....... . . ...... ..... ............... .. .. ..... ............................................................... . . . .. . ........................................................................ .

Dnition 3.1 Une fonction g(n) est dite O(f (n)) sil existent des constantes co e et no telles que g(n) est infrieure ` co f (n) pour tout n > no . e a Ainsi lalgorithme A1 est dordre O(n2 ) et lalgorithme A2 est dordre O(n). La notion dordre est donc proportionnelle au temps dexcution sans tre encombr e e e de caractristiques particuli`res dune machine tout en restant indpendante dun e e e input particulier.
e Exemple 3.4 Un algorithme qui ncessite complexit dordre e O(n3 ).
n3 3

+ n2 + 2 n oprations lmentaires a une e ee 3

Exemple 3.5 Soit lalgorithme rcursif maxmin qui recherche llment maximum et e ee
llment minimum dun vecteur de 2n lments. ee ee
function [a,b]=maxmin(s) if length(s)==2 if s(1) > s(2), a = s(1); b=s(2); else, a = s(2); b=s(1); end else k=length(s); [a1,b1]=maxmin(s(1:k/2)); [a2,b2]=maxmin(s((k/2)+1:k)); if a1 > a2, a=a1; else, a=a2; end if b1 < b2, b=b1; else, b=b2; end end

On dsire tablir sa complexit en comptant le nombre de comparaisons eectues. Le e e e e nombre de comparaisons en fonction de n est donn par la rcurrence suivante : e e C2n = 1 pour n = 1 2 C2n1 + 2 pour n > 1 .

Rsolvons alors cette rcurrence : e e C21 C22 C23 C24 C2n = 1 = 2 C21 + 2 = 2 + 2 = 2 C22 + 2 = 22 + 22 + 2 = 2 C23 + 2 = 23 + 23 + 22 + 2 . . . = 2 C2n1 + 2 = 2n1 + 2n1 + 2n2 + + 21
n1 k=1

= 2n1 + = 3 n 2 2 . 2

1 2k = 2n1 + 2n 2 = 2n (1 + ) 2 2

En posant N = 2n le nombre de comparaisons devient 3 N 2 et la complexit de lalgoe 2 rithme est O(N ).

3.4

Oprations lmentaires (ops) e ee

Pour les algorithmes rsolvant des probl`mes du domaine de lalg`bre linaire, on e e e e utilise souvent le ops (oating point operations) comme mesure de dnombrement e des oprations lmentaires. Chacune des quatre oprations lmentaires addition, e ee e ee soustraction, multiplication et division compte pour un op2 . Soient les vecteurs x, y Rn , z Rm et les matrices A Rmr et B Rrn on a alors : Opration ops e x+y n xx 2n (n > 1) xz nm AB 2 m r n (r > 1)

3.5

Classication des algorithmes

La notion de complexit est tr`s importante, car elle permet une classication des e e algorithmes et des probl`mes. On distingue notamment deux classes, suivant que le e nombre doprations lmentaires est une fonction : e ee polynmiale de la taille du probl`me, o e non-polynmiale de la taille du probl`me. o e Seuls les algorithmes appartenant ` la premi`re classe peuvent tre qualis defa e e e caces. Pour les probl`mes de la deuxi`me classe, on ne peut garantir dobtenir e e
2

Dans certains ouvrages un op est dni comme zi = xi + yi . e

la solution dans un intervalle de temps raisonnable. Ceci est illustr dans les tae bleaux 3.1 et 3.2. Dans le tableau 3.1 on donne le nombre doprations lmentaires pour rsoudre des e ee e probl`mes de taille variable en fonction de direntes complexits couramment rene e e contres. On trouve aussi linformation sur la relation entre oprations lmentaires e e ee et temps de calcul qui correspondent ` une puissance de calcul de 106 oprations a e lmentaires par seconde (1 Mops). On voit notamment laccroissement extraordiee naire du nombre doprations lmentaires (et du temps de calcul) lorsque la come ee plexit du probl`me devient non-polynmiale. e e o Tab. 3.1 Flops en fonction de la complexit. e
Complexit e log2 n n n log2 n n2 n3 2n n! 2 1 2 2 22 23 22 2 n (Taille du probl`me) e 8 128 1024 3 7 10 23 27 210 3 23 7 27 10 210 6 14 2 2 220 9 21 2 2 230 28 5 213 2128 5 2714 21024 7 28766

1 Mops 226 ops/min, 232 ops/h, 236 ops/jour, 245 ops/anne, e 252 ops/siecle.

Le tableau 3.2 montre quun accroissement de la performance des ordinateurs a un eet ngligeable sur le gain de la taille des probl`mes qui peuvent tre envisags pour e e e e une solution. En eet les accroissements de performance que lon peut sattendre sont dordre polynomial et ne modient pas lordre des temps de calcul. Tab. 3.2 Eet dun accroissement de vitesse sur la taille du probl`me. e
Complexit e n n2 2n 8n Taille maximale N1 N2 N3 N4 Accroissement de vitesse 8 128 1024 8 N1 128 N1 1024 N1 2.8 N2 11.3 N2 32 N2 N3 + 3 N3 + 7 N3 + 10 N4 + 1 N4 + 2.3 N4 + 3.3

A titre dexemple, la performance de quelques ordinateurs et la date de leur mise sur le march. On remarque quen 1985 le rapport de performance entre un ordinae teur personnel et un ordinateur de haute performance (Cray) tait de mille. Aussi e voit-on que la performance des ordinateurs personnels a t multiplie par mille ee e approximativement en quelque 15 ans.

10 10

Machine 8086/87 a 8 MHz ` Cray XMP Pentium 133 MHz Pentium II 233 MHz Pentium III 500 MHz Pentium IV 2.8 GHz Earth Simulator

Mops 0.03 33 12 28 81 1000 30 TFlops

Anne e 1985 1985 1995 1997 1999 2003 2003

10 10 10 10 10

1984 1987 1990 1993 1995 1998 2001 2004

Fig. 3.1 Evolution de la performance des PC.

Chapitre 4 Resolution numrique de syst`mes e e linaires e

Le probl`me consiste ` trouver la solution du syst`me dquations e a e e a11 x1 + a12 x2 + + a1n xn = b1 a21 x1 + a22 x2 + + a2n xn = b2 . . . . . . . . . . . . . an1 x1 + an2 x2 + + ann xn = bn En notation matricielle ce syst`me scrit e e Ax = b et on rappelle quil existe une solution unique si et seulement si det(A) = 0. Rappel Soit Ax = b avec x Rn , r = rang(A) et B = [A b] avec r = rang(B). Alors solution de est possible si r = r est unique si r = r = n Ax = b comporte une innit de solutions si r = r < n . e x = A1 b ce qui sugg`re que lon a besoin de linverse de la matrice A pour calculer x. e Ceci constitue une mthode numrique tr`s inecace ( 3 fois plus de calculs que e e e ncessaires et 2 fois plus despace mmoire pour moins de prcision). e e e En gnral on na quasiment jamais besoin de linverse dans les probl`mes dalg`bre e e e e 1 linaire. Il nest pas non plus indiqu de calculer A pour rsoudre des syst`mes e e e e x(i) = A1 b(i) 23 i = 1, . . . , k

La solution dun syst`me linaire se note e e

mme si k est tr`s grand. e e En r`gle gnrale, si lon na pas besoin danalyser les lments de A1 il est tr`s e e e ee e 1 vraisemblable quon puisse et doive viter le calcul de A dans le probl`me dalg`bre e e e linaire considr. e ee Le choix dune mthode particuli`re pour la rsolution dun syst`me linaire dpendra e e e e e e galement du type et de la structure de la matrice A intervenant dans le probl`me. e e La matrice A peut tre caractrise par une structure : e e e creuse, diagonale, avec concentration des lments non-nuls autour de la diagonale ee (aij = 0 si |i j| > b, b est appele la largeur de la bande), e matrices par blocks, tridiagonale (cas particulier dune matrice par bande avec b = 1), triangulaire, ou par une quantication particuli`re, comme par exemple les matrices : e symtriques (A = A ), e dnies positives (x Ax > 0, x = 0), e semi dnies positives (x Ax 0, x = 0), e matrices de Tplitz (T Rnn est une matrice de Tplitz sil existe des scalaires o o rn+1 , . . . , r1 , r0 , r1 , . . . , rn1 , tels que Tij = rji quelque soit i et j), matrices Hessenberg (une matrice de Hessenberg est une matrice triangulaire plus une diagonale immdiatement adjacente non nulle), e Hankel, Vandermonde, . . . etc. Le dterminant de A joue un rle central dans les considrations thoriques, mais e o e e nest daucune utilit pour la rsolution numrique (sauf rares exceptions). La r`gle e e e e de Cramer xi = det(Ai ) det(A)

est incroyablement inecace pour rsoudre un syst`me linaire si lon calcule le e e e dterminant selon le schma de lexpansion des mineurs. e e Souvent on entend dire que lon ne peut rsoudre un syst`me donn car le dterminant e e e e est nul. En utilisant des mthodes ecaces de calcul on dcouvre que le syst`me ne e e e peut tre rsolu avant davoir la valeur du dterminant. e e e La valeur du dterminant ne contient aucune information sur les probl`mes numriques e e e que lon peut rencontrer en rsolvant Ax = b. e

Exemple 4.1 Soit deux matrices A et B dordre 40 et vriant la structure ci-apr`s : e e

A= .1 .1 .. . .1 .1 1 1 1 1 1 1 . .. .. . B= . . . 1 1 1

On a det(A) = 1040 et det(B) = 1. La matrice A ne pose videmment aucun probl`me, e e alors que la matrice B est tr`s mal conditionne (cond(B) 1013 ). e e Multiplier les m quations du syst`me Ax = b par 100 modie la valeur du dterminant e e e m . Ceci change la valeur du dterminant sans modier les calculs numriques. de 100 e e

Dans les chapitres qui suivent on discutera des techniques qui consistent ` transfora mer le probl`me original en un probl`me ayant une forme particuli`re, pour laquelle e e e la solution dvient triviale ou tr`s simple. Par exemple, on transforme le syst`me e e e Ax = b en Ux = c avec U triangulaire (ou diagonale) et ainsi x est facile ` calculer. a Dans dautres situations on peut manipuler Ax = b de la sorte ` obtenir le syst`me a e Rx = c avec R orthogonale, ce qui se rsoud facilement par rapport ` x. e a Avertissement Les algorithmes prsents dans les chapitres qui suivent ne sont pas spciques ` e e e a un language ou une machine prcise. Ils sont crits dans un pseudo-language de e e programmation que lon pourra traduire dans nimporte quel language de programmation formel. Ce pseudo-language est tr`s proche du language Matlab sans toutefois e respecter exactement la syntaxe de Matlab (p. ex. la produit AB est cod A*B dans e Matlab, ou encore A = B est cod A~=B). e

Chapitre 5 Syst`mes triangulaires e

Beaucoup de mthodes de rsolution transforment le syst`me original en un syst`me e e e e triangulaire qui poss`de la mme solution. De ce fait on est souvent amen ` manie e ea puler des syst`mes triangulaires dans les algorithmes. Ainsi on discutera dabord la e rsolution des syst`mes triangulaires. e e Considrons le syst`me triangulaire infrieur suivant : e e e 11 0 21 22 x1 x2 = b1 b2 .

Si les lments 11 , 22 = 0, alors les inconnues peuvent tre dtermins de faon ee e e e c squentielle, cest-`-dire e a x1 = b1 /11 x2 = (b2 21 x1 )/22 .

5.1

Forward substitution

Pour un syst`me Lx = b, avec L triangulaire infrieure, la solution de la i-`me e e e quation scrit : e e

i1
... . . . . . . . . . . . . . . . . . . . ... . . . . .. . ... . . . .. . . ... . . .. . .. . . . .. .. . . . . . . ............... ... ii . . ................. ... . . . . . . . . . . . . .. . . . . .. .... . ........ .. .. . ................. . . . . . . . . . . . . . . . .. . . . ...... .. .. . . . .. .. .. . ... . . .. . .. . . .. . .. . . . .. .. . . . .. . . .. . . .. . .. . . .. . .. . . .. . . .. . . .. . ................................... . . .....................................

xi =

ij xj
j=1

. .. . .. . .. . .. . .. . .. . .. . .. . .. . . .. .. .. ..... ... .. ... ... .

xi =

. ...... ..... .. . . .... .. ..

Cette procdure sappelle forward substitution. On remarque que bi nest utilis e e que pour le calcul de xi (voir aussi la gure), ainsi on peut remplacer bi par xi . Lalgorithme 3 rsoud le syst`me triangulaire infrieur Lx = b dordre n par forward e e e substitution et remplace b avec la solution x. Lalgorithme 3 ncessite n2 + 3n 4 ops. Sa complexit est donc dordre O(n2 ). e e 27

Algorithme 3 fsub
1: b1 = b1 /L1,1 2: for i = 2 : n do 3: bi = (bi Li,1:i1 b1:i1 )/Lii 4: end for

Dmonstration 5.1 1 op pour la premi`re instruction. Pour i = 2 on a 5 ops (1 op e e

pour le produit scalaire si les vecteurs sont de longueur 1, 2 ops pour calculer les indices et 2 ops pour la soustraction et la division). Pour i > 2 on a n (2(i1)+4) = n2 +3n10 i=3 ops (2(i 1) ops pour le produit scalaire si les vecteurs sont de longueur > 1). Le total est alors n2 + 3n 4 ops pour n 3.

5.2

Back substitution

On peut rsoudre un syst`me triangulaire suprieur Ux = b de faon analogue ` e e e c a celle prsente avant, soit e e
n
..................................... . .................................... . .... . . . . . .... . . . ... . .. . . . .. .. . . . .. . .. . . . .. . . . . . . .... ..................... . ........ii......... . . .. .............. . ... ... . . . . . . . . . ....... . . ..... . ...................... . . . . . . . . ........... ... .. ... . .. . .. . . . . .. . . .. . . .. . . .. . .. . . .. . . .. . . .. . . .. . .. . . . .. .. . . .. . . .. . . .. . .. . . . .. . . . . . . . . . . . . . . . . . . ..

xi =

uij xj
j=i+1

uii

. . ..... .... . . ... .. . . .. .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . ..

xi =

. ..... ..... .. . . .... .. ..

et que lon appelle back-substitution. On peut ` nouveau remplacer bi par xi . a Lalgorithme 4 rsoud le syst`me triangulaire suprieur Ux = b dordre n par back e e e substitution et remplace b avec la solution x. Algorithme 4 bsub
1: bn = bn /Un,n 2: for i = n 1 : 1 : 1 do 3: bi = (bi Ui,i+1:n bi+1:n )/Uii 4: end for

La complexit de lalgorithme 4 est dordre O(n2 ). e

5.3

Proprits des matrices triangulaires unitaires e e

Dnition 5.1 On appelle une matrice triangulaire avec une diagonale unitaire une e matrice triangulaire unitaire. On rappelle ici quelques proprits du produit et de linverse des matrices trianguee laires. Linverse dune matrice triangulaire infrieure (suprieure) est triangulaire infrieur e e e (suprieur). e

Le produit de deux matrices triangulaires infrieures (suprieures) est triangulaire e e infrieur (suprieure). e e Linverse dune matrice triangulaire infrieure (suprieure) unitaire est triangue e laire infrieure (suprieure) unitaire. e e Le produit de deux matrices triangulaires infrieures (suprieures) unitaires est e e triangulaire infrieure (suprieure) unitaire. e e

Chapitre 6 Factorisation LU
Cest la mthode de prdilection pour la rsolution de syst`mes linaires comportant e e e e e une matrice A dense, sans structure et sans quantication particuli`re. e On vient de voir quil est tr`s simple de rsoudre des syst`mes triangulaires. La e e e transformation du syst`me original en un syst`me triangulaire se fait en choisissant e e des combinaisons linaires appropries des quations du syst`me. e e e e Exemple 6.1 Soit le syst`me e
3x1 + 5x2 = 9 6x1 + 7x2 = 4 si lon multiplie la premi`re quation par 2 et que lon la soustrait de la deuxi`me on e e e obtient 9 3x1 + 5x2 = . 3x2 = 14

Cette procdure sappelle limination de Gauss. Dans la suite on donnera ` ce e e a procd une formulation matricielle, notamment en termes de factorisation dune e e matrice. On prsentera un algorithme qui factorise une matrice donne A en une e e matrice triangulaire infrieure L et une matrice triangulaire suprieure U, tel que e e A = LU. Exemple 6.2 Pour la matrice A correspondant a lexemple prcdent la factorisation ` e e
A = LU est : 3 5 6 7
A

1 0 2 1
L

3 5 0 3
U

La solution du syst`me original Ax = b est alors obtenue en rsolvant deux syst`mes e e e triangulaires successivement. On remplace A par sa factorisation Ax = L Ux = Ly = b
y

et lon cherche la solution y ` partir du syst`me triangulaire infrieur Ly = b, puis a e e la solution de x ` partir du syst`me triangulaire suprieur Ux = y. a e e

6.1

Formalisation de llimination de Gauss e

Il sagit de formaliser une procdure qui transforme une matrice en une matrice e triangulaire suprieure en liminant, colonne apr`s colonne, les lments non nuls en e e e ee dessous de la diagonale comme illustr dans lexemple 6.3. e Exemple 6.3 Soit la matrice
1 4 7 2 5 8 3 6 10

Lors de la deuxi`me tape on transforme cette nouvelle matrice en soustrayant 2 fois la e e deuxi`me ligne de la troisi`me ligne pour obtenir la forme triangulaire recherche. e e e 1 4 7 0 3 6 . 0 0 1

que lon veut transformer en une matrice triangulaire suprieure. On obtient cette forme en e deux tapes. Lors de la premi`re tape on soustrait 2 fois la premi`re ligne de la deuxi`me e e e e e ligne et 3 fois la premi`re ligne de la troisi`me ligne pour obtenir la matrice e e 1 4 7 0 3 6 . 0 6 11

Matrice dlimination M e Formulons la transformation en matrice triangulaire comme une transformation nonsinguli`re dun syst`me linaire. Considrons dabord le vecteur ` deux lments e e e e a ee x = [x1 x2 ], alors si x1 = 0, on peut dnir la transformation e 1 0 x2 /x1 1 x1 x2 = x1 0 .

De mani`re gnrale, tant donn un vecteur x avec xk = 0, on peut annuler tout e e e e e les lments xi , i > k, avec la transformation ee
1 . . .. . . 0 Mk x = 0 . . .. . . 0 0 . . . 1 (k) k+1 . . . n
(k)

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

0 x1 . . . . . . 0 xk 0 xk+1 . . . . . . xn 1

x1 . . . xk = 0 . . . 0

o` i = xi /xk , i = k + 1, . . . , n. Le diviseur xk est appel le pivot, la matrice u e (k) Mk est appel la matrice de transformation de Gauss et les i constituent les e multiplicateurs de Gauss. Proprits des matrices Mk e e Proprit 6.1 Les matrices Mk sont triangulaires infrieures avec diagonale unie e e taire et de ce fait elles sont non-singuli`res. e Proprit 6.2 Les matrices Mk peuvent scrire comme e e e
1

(k)

Mk = I

(k) ek

... ... .... ... . . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. (k) . .. . .. . .. . .. . ..

e k

e e avec (k) = [0, . . . , 0, k+1, . . . , n ] et ek la k`me colonne de la matrice identit. Proprit 6.3 Linverse de la matrice Mk scrit e e e
1 Mk = I + (k) ek

(k)

qui peut tre dduit directement de Mk en inversant simplement le signe des multiplie e cateurs de Gauss. On vrie (I (k) ek )(I+ (k) ek ) = I (k) ek + (k) ek (k) ek (k) ek = e I tant donn que ek (k) = 0. e e Proprit 6.4 Soient deux matrices dlimination Mk et Mi avec i > k, alors e e e Mk Mi = I (k) ek (i) ei + (k) ek (i) ei = I (k) ek (i) ei tant donn que ek (i) = 0. Ainsi le produit peut tre considr comme lunion e e e ee des deux matrices.
1.. .
.. .. .. .. .. .. .. .. .. .. . 1 .... . . . . .....1. . . . . ... ... ..... .... . .. . .. .. . .. ..... . ... . .. . . ... . .. ..... . ... . .. . .. .. .. . .. ..... . .. . .. . .. . ... . ..... . .. .. . .. . .. .. . .. ..... . .. . .. . . ... . .. .. ..... . . .. . .. .. .. . .. ..... . . .. . ..... . ..... . . .. . 1 . .. .

1.. .

.. .. .. .. .. .. .. .. .. .. . .1. . . . . . ..1. . .. . . .. . .. . .. ...... . .. . . . .. .. . . . .. ... ... . . .. . . .. . . . .. . . .. ... ... . .. . . .. . .. .. .. . .. ... . . .. . .. .. . .. . . . . .. ... .. . .. . .. . . .. .. .. . .. ... . .. . .. . . .. . . . .. .. . .. ... . . .. . .. . .. .. .. . .. ... .. . . .. . . .. .. . .. ... . . .. . 1 . .. .

1.. .

.. .. .. .. .. .. .. .. .. .. 1 .... . . . . .....1. . .. . . ... . . . ..... ...... . .. . . .. .. .. . . ..... ... ... . . .. . . ... . . . . . ..... ... ... .. . .. . . .. .. .. .. ..... ... .. . . .. . .. . ... . . . .. ..... ... . .. . .. . .. .. .. .. ..... ... . .. . .. . .. . ... . . . .. ..... ... .. . . .. . .. .. .. .. ..... ... . .. . . .. . ..... .. ..... ... . . .. . 1 . .. .

Mk+1

I+ (k) e + (k+1) e k k+1

1 La mme proprit est vrie pour le produit Mk Mi1 . e ee e e

Revenons ` la matrice A ` factoriser a a Considrons maintenant la matrice A(k1) dont les premi`res k 1 lignes ont dj` e e ea la forme triangulaire suprieure. Il sagit alors de transformer la colonne k de la e matrice, que nous dsignerons ici comme le vecteur x, de la sorte que xi = 0 pour e i = k + 1, . . . , n.
k 1:k1
....................................................... .. . .................................................... .. . . . . . . . . . . . . . . . . . . . . . . . . . ...... . .. . . . . . . . ... .. . . . . . . . . . . . . . ..... .... .. . . ............................................. .. ... .. .. . . . . . . . . ... . . . . . . . . . . ..... ... . . . . ... .. . . . . . . . .... . .... .. . . . . . . . ... . . . . . . ........ . . . . ... . . . ... .. . . . . ...... . . . . . ... . . . . . ... . . . .. ...... . . . . . ................................ ......... . ....... . . . .. . .............. ....................... . . .. . ... ... ... . . . . . . . . . . . .. .................................. .. . .. . . .. . .. ........... .............. .. . ... . .... ............. . . .............................................. ... ... . .... . .. .. .. .... ..... .... ........... .. . . . .. .............................................. .. .. ... .. . .. . ..... . ... .... .. . . . . . ...................... . ... . . . . . . . . . . . . ...................... . . . . ...................... . . .. . . . . . . . . . . . . . . . ...................... . . . . . . ...................... . ... . . . . . . . . .. . . . . . . . . . . . . . . ...................... . . .

vecteur x A(k) =

A(k1) =

(k1) A k:n,k:n

....................................................... .. . .................................................... .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . ... .. . . . . . . . . . . . . . .. . . ................................................ .. . . . . . . . . ... . . . . . . . . . . . . . ... . . . . ... .. . . . . . . . . . . . . . .. . . . . . . . ... . . . . . . . . . . . . . . . ... . . . ... .. . . . . . . . . . . . . . ... . . . . . ... . . . . . . . . . . . . . . .. ....................................... . . .. . . ........... ........................... . . .. . ... .. . . . . . . . . . . . . . .. .................................. .. . . .. . .. .......... .............. .. ... ....... ............. ............................................... . . ... ... ... ..... ... ... ... ... . . ............................. . . . .. ............................................ . ....................... . . .................. . . . . . ... . ... . . . . . .. . . ... . ... . . . . ....................... . ....................... . . ... ... . . ....................... . ................... .. . . ... . ... . . . . . ................. . . . . . ... . ... . . . . . . . . ................... .. . . ....................... . ....................... . . ....................... .. . . ... . ... . . . . . ................. . ................... . . . . . . ... . ... . . . .. . . ... . ... . . . . ....................... . ....................... . . . . . . ... . ... . . . . . ................. .. . . ... . ... . . . . ....................... . . ....................... . ................... .. . . ... . ... . . . . . ................. . . ................... . . ... . ... . . . . ....................... . . ................. . . .. .

Cette transformation sop`re en faisant le produit Mk A(k1) et seulement les colonnes e j > k sont aectes par la transformation. e Le produit de Mk par la matrice A(k1) sav`re tre particuli`rement simple. Seulee e e (k1) ment la matrice Ak:n,k:n intervient dans cette transformation. Comme le rsultat e de la transformation de Gauss est connu pour la colonne Ak+1:n,k (zros), seuls les e lments de Ak+1:n,k+1:n doivent tre calculs. ee e e En fait la transformation Mk A(k1) peut tre formule comme une transformation e e (k1) M1 applique ` la sous-matrice Ak:n,k:n. Lalgorithme 5 applique cette transformation e a M1 ` une matrice C. a Soit une matrice C Rmm , lalgorithme 5 remplace C2:m,1 par le vecteur des multiplicateurs (1) , et la sous-matrice C2:m,2:m par le produit M1 C. Les lments de ee C1,1:m restent inchangs. e Algorithme 5 (tgm1)
1: 2: 3: 4: 5: function C = tgm1(C) m = size(C, 1) if C11 = 0 then Pivot est nul, arrt e C2:m,1 = C2:m,1 /C11 C2:m,2:m = C2:m,2:m C2:m,1 C1,2:m
(k) (k)

Lalgorithme 5 ncessite 2(m1)2 +m1 = (2 m1)(m1) oprations lmentaires. e e ee Les multiplicateurs conservs seront utiliss dans une tape ultrieure. e e e e Dtaillons lalgorithme 5. La transformation scrit e e M1 C = (I (1) e1 ) C = C (1) e1 C

ce qui correspond ` a
. . . . . ............................. . . . . . . . . . . . . . . . . . . . . . . . . ...... ..... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... ... .... . . . . 1. . .

e 1 (1)

. . . . . ............................. . . . . . . . . . . . . . . . . . . . . . . . .

C1,1:m

et ce qui explique linstruction 5 de lalgorithme 5. Exemple 6.4 Soit la matrice C a laquelle on applique la transformation M1 C. `
2 1 3 C = 4 5 6 , 6 7 8 C2:3,1 = 2 3 alors 1 0 0 2 1 3 2 1 3 M1 C = 2 1 0 4 5 6 = 0 3 0 3 0 1 6 7 8 0 4 1 5 6 7 8 2 3 1 3 = 3 0 4 1

C2:3,2:3 =

En appliquant n 1 fois dans lordre M1 , . . . , Mn1 la transformation de Gauss ` a une matrice A dordre n Mn1 M2 M1 A = U on obtient une matrice U qui est triangulaire suprieure. Ce procd est appel e e e e limination de Gauss. e Exemple 6.5 Soit la matrice A de lexemple du dbut de la section prcdente e e e
1 2 3 A = 4 5 6 , 7 8 10 1 0 0 M1 = 4 1 0 , 7 0 1 1 2 3 M1 A = 0 3 6 , 0 6 11 A(1) 1 2 3 = 0 3 6 . 0 0 1

1 0 0 M2 = 0 1 0 , 0 2 1

M2 M1 A = M2 A(1)
A(2)

Rsumons ltape k de llimination de Gauss : e e e On a une matrice A(k1) = Mk1 M1 A qui a une forme triangulaire suprieure e des colonnes 1 ` k 1. a (k1) Les multiplicateurs dans la matrice Mk sont calculs ` partir du vecteur Ak+1:n,k . e a Pour quils soient dnis il faut que Ak,k e
(k1) Llment Ak,k ee (k1) (k1) Ai,k /Ak,k est (k1)

= 0.
(k)

qui intervient dans le calcul des multiplicateurs de Gauss i

appel le pivot. Le pivot doit tre dirent de zro et sa grandeur e e e e relative est dimportance pour la prcision des calculs. e

Comment identier la matrice L A partir de llimination de Gauss e Mn1 M1 A = U qui conduit ` une matrice triangulaire suprieure U on dduit, en crivant a e e e (Mn1 M1 )1 (Mn1 M1 ) A = (Mn1 M1 )1 U
I L

que la matrice L dans la factorisation A = LU est constitu par le produit e

1 1 L = (Mn1 M1 )1 = M1 Mn1 .

A partir des proprits 6.3 et 6.4 on vrie aisment que ee e e

L = (I +

(1) e1 ) (I

(n1) en1 )

=I+
k=1

(k) ek

et que les lments de la matrice L correspondent aux multiplicateurs de Gauss qui ee ont t conservs par lalgorithme 5 (instruction 4). ee e Existence de la factorisation LU Il appara dans lalgorithme 5 que dans le cas dun pivot nul, la factorisation nexiste t pas toujours. Le thor`me qui suit, montre quun pivot nul sidentie avec une souse e matrice diagonale singuli`re. e Thor`me 6.1 Existence : Une matrice A dordre n admet une factorisation LU si e e elle vrie det(A1:k,1:k ) = 0 pour k = 1, . . . , n 1. Unicit : Si la factorisation existe e e et A est non-singuli`re, alors elle est unique et det(A) = n uii . e i=1 Exemple 6.6 Matrice A non-singuli`re avec pivot nul. e
1 3 2 1 0 0 u11 u12 u13 3 9 5 = 21 1 0 0 u22 u23 2 5 4 31 32 1 0 0 u33

En dveloppant le produit de la matrice triangulaire infrieure avec la matrice triangulaire e e suprieure on dduit que : e e a11 = 1 u11 = 1 u11 = 1 a21 = 21 u11 = 3 21 = 3 a31 = 31 u11 = 2 31 = 2 a12 = 1 u12 = 3 u12 = 3

a22 = 21 u12 + 1 u22 = 9 u22 = 0

Cependant le produit de la troisi`me ligne avec la deuxi`me colonne e e a32 = 5 = 31 u12 + 32 u22 = 6 conduit a une contradiction. On a det(A) = 0, mais det(A1:2,1:2 ) = 0. `

Finalement pour une matrice A dordre n et vriant les conditions du thor`me 6.1 e e e la factorisation en une matrice triangulaire infrieure L et une matrice triangulaire e suprieure U est obtenu avec lalgorithme 6. Cet algorithme remplace A avec U et L e ` lexception des lments de la diagonale de L (les lments diagonaux de L valent a ee ee 1). Algorithme 6 (eg) Elimination de Gauss.
1: function A = eg(A) 2: for k = 1 : n 1 do 3: Ak:n,k:n = tgm1(Ak:n,k:n ) 4: end for
n1 Lalgorithme 6 comporte k=1 2(n+1k)1 (nk) = 2 n3 1 n2 1 n+1 oprations e 3 2 6 lmentaires. Il sagit de la formulation classique de llimination de Gauss. ee e

On va montrer par la suite que cet algorithme nest pas numriquement stable. e Choix du pivot Lors de la prsentation de llimination de Gauss on a dni le pivot. Il sagit de e e e llment akk de la matrice A au dbut de ltape k. Pour pouvoir calculer les multiee e e plicateurs de Gauss, le pivot doit tre dirent de zro. Etant donn que les calculs se e e e e font avec une arithmtique en prcision nie, la grandeur relative du pivot inuence e e la prcision des rsultats. e e Exemple 6.7 Soit le syst`me e
0.0001 x1 + 1 x2 = 1 1.0000 x1 + 1 x2 = 2 (1) (2)

pour lequel la solution exacte est x1 = 1.0001 . . . et x2 = 0.9998 . . .. Appliquons llimination e de Gauss avec une arithmtique en base 10 et t = 3. e On a a11 = 0.0001 le pivot et 1 = 1/.0001 = 10000 le multiplicateur de Gauss. Lquation e (2) devient alors 1 x2 = 2 1 x1 + 1 x1 + 10000 x2 = 10000 9999 x2 = 9998 (2) eq. (1) (1 )

Avec 3 digits signicatifs on obtient alors x2 = 1. On remplace x2 dans (1) et on obtient x1 = 0. (Catastrophe numrique ! !). e Permutons les quations : e 1.0000 x1 + 1 x2 = 2 0.0001 x1 + 1 x2 = 1 (2) (1)

On a a11 = 1 le pivot et 1 = .0001/1 = .0001 le multiplicateur de Gauss. Lquation (1) e devient maintenant 0.0001 x1 + 1.0000 x2 = 1.0000 0.0001 x1 0.0001 x2 = 0.0002 0.9999 x2 = 0.9998 (1) eq. (2) (1 )

Avec 3 digits signicatifs on obtient alors x2 = 1. On remplace x2 dans (2) et on obtient x1 + 1 = 2 do` x1 = 1. (Excellent pour une prcision de 3 digits ! !). u e

Les petits pivots sont ` lorigine de coecients relativement grands dans la matrice a triangularise. La factorisation de la matrice A de lexemple prcdent donne e e e A= .0001 1 1 1 = 1 0 10000 1 .0001 1 0 9999 = LU .

Lors de la rsolution du syst`me triangulaire on eectue alors des soustractions entre e e nombres relativement grands ce qui peut conduire ` une perte de prcision fatale. a e 0 1 En appliquant la permutation P = 1 0 on obtient PA = 1 1 .0001 1 = 1 0 .0001 1 1 1 0 .9999 = LU

ce qui modie lamplitude des lments dans la matrice U. ee On peut rencontrer des pivots tr`s petits sans que le probl`me soit mal condie e 1 tionn (e.g. pour A = 1 0 on a cond(A) = 1). Ainsi, llimination de Gauss, telle e e que prsent jusquici, est un algorithme numriquement instable. Il est donc absoe e e lument ncessaire dintroduire dans lalgorithme un mcanisme de permutation des e e lignes et/ou des colonnes an dviter des petits pivots. e

6.2

Matrices de permutation

La permutation des lignes dune matrice peut se formaliser en criture matricielle e avec une matrice de permutation qui nest rien dautre quune matrice didentit e avec des lignes rordonnes. La matrice suivante est un exemple de matrice de pere e mutation :

0 1 P = 0 0

0 0 0 1

0 0 1 0

1 0 . 0 0

Le codage dune matrice de permutation dordre n se fait ` laide dun vecteur p de a longueur n. On aura p(k) gal ` lindice de lunique colonne non nulle de la ligne k. e a La matrice de lexemple ci-dessus sera alors cod p = [4 1 3 2]. e Soit P une matrice de permutation, alors le produit P A permute les lignes de A, le produit AP permute les colonnes de A, P est orthogonale, P 1 = P , le produit de deux matrices de permutation engendre une matrice de permutation.

Dans ce qui suit on sera plus particuli`rement intress ` la permutation de deux e e e a lignes (ou colonnes) dune matrice. La matrice de permutation particuli`re qui core respond ` une opration de ce type est la matrice dchange. Il sagit de la matrice a e e didentit dans laquelle on a chang deux lignes, par exemple e e e
0 0 E= 0 1 0 1 0 0 0 0 1 0 1 0 0 0

Le produit EA change la ligne 1 avec la ligne 4 dans A. Le produit AE change la e e colonne 1 avec la colonne 4 dans A. Pour la suite nous sommes particuli`rement intresss a changer les lignes dans e e e ` e lordre, cest-`-dire on change dabord la premi`re ligne avec une autre ligne, puis a e e la deuxi`me ligne avec une autre ligne, puis la troisi`me ligne, etc. e e A titre dillustration considrons un vecteur x R4 pour lequel on change dabord e e la ligne 1 avec la ligne 3, puis la ligne 2 avec la ligne 3 et nalement la ligne 3 avec la ligne 4. Le produit des matrices qui correspond ` cet change est : a e
1 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 x1 0 x2 0 x3 1 x4 0 1 = 0 0 0 0 0 1 1 0 0 0 0 x1 0 x2 1 x3 0 x4 x3 x1 = x4 x2

ce que peut tre formalis comme e e E3 E2 E1 x ou les matrices Ek sont dnies ` laide du vecteur e(k), k = 1, 2, 3 comme : e a Ek est la matrice didentit ou la ligne k est chang avec la ligne e(k). e e e

Pour notre exemple e = [3 3 4]. La matrice de permutation P qui correspond ` ces a trois changes e E3 E2 E1 x = P x
P

est dnie par p = [3 1 4 2] . e La permutation dun vecteur x Rn qui correspond aux r changes successifs e Er E1 x et dnis par le vecteur e Rr peut tre obtenu par la procdure suivante e e e qui remplace x : Algorithme 7 (perm)
1: function x = perm(e, x) 2: for k = 1 : r do 3: x(k) x(e(k)) 4: end for
. . . . . . . . .

e(k)
. . . . . . . . . . . . . . . . . .

Le symbole signie changer le contenu : e

1: s = x(k) 2: x(k) = x(e(k)) 3: x(e(k)) = s

. . . . . . . .. . . . . . . . . . .. . . . .. . . ....... .. . . .......... . . . ...... ....... 2: ...... . . . ... . 1: ............... . . ......... . . . . . .. . .. . .. .. .... 3: .... ........................ s .......................

Le vecteur p qui dnit la matrice de permutation P qui correspond aux changes e e de lignes successifs dnis par un vecteur e peut tre obtenu avec la mme procdure e e e e dans laquelle le vecteur x est remplac par le vecteur 1 : n, soit p = perm(e, 1 : n). e Exemple 6.8 Soit le vecteur x = [19 20 13 55] et le vecteur e = [3 3 4]. En appliquant la procdure perm a x et a un vecteur [1 2 n] on obtient e ` `
Etat du vecteur x 19 20 13 13 20 19 13 19 20 13 19 55 Etat 1 3 3 3 du vecteur 1 : n 2 3 4 2 1 4 1 2 4 1 4 2

k 1 2 3

e(k) 3 3 4

55 55 55 20

k 1 2 3

e(k) 3 3 4

0 1 Px = 0 0

0 0 0 1

1 0 0 0

0 19 0 20 1 13 0 55

13 19 = 55 . 20

Si pour une permutation donne y = Er E1 x on veut revenir au vecteur original e on doit procder aux changes des lignes dans lordre inverse : e e x = (Er E1 )1 y .

Etant donn que les matrices Ek , k = 1, . . . , r sont orthogonales on vrie (Er E1 )1 = e e 1 E1 Er . La procdure qui remplace x par le produit (Er E1 ) x sera : e
for k = r : 1 : 1 do x(k) x(e(k)) end for

On peut facilement vrier cette procdure avec lexemple introduit plus haut. e e

6.3

Pivotage partiel

Montrons comment changer les lignes lors de la factorisation LU de la sorte ` ce e a quaucun multiplicateur ne devienne suprieur ` 1 en valeur absolue. e a Exemple 6.9 Soit la matrice A que lon dsire factoriser : e
3 17 10 A = 2 4 2 . 6 18 12

Pour obtenir les multiplicateurs les plus petits possibles a11 doit tre llment le plus e ee grand de la premi`re colonne. Pour changer la premi`re ligne avec la troisi`me ligne on e e e e utilise la la matrice dchange e 0 0 1 6 18 12 E1 = 0 1 0 et on forme le produit E1 A = 2 4 2 . 1 0 0 3 17 10 On calcule alors M1 (tous les multiplicateurs sont infrieurs a 1 en valeur absolue) e ` 1 0 0 6 18 12 M1 = 1 1 0 M1 E1 A = 0 2 2 . 3 1 0 1 0 8 16 2 0 0 1 0 1 et M2 = 0 1 0 0

Pour obtenir a ` 3. On a 1 E2 = 0 0

nouveau les plus petits multiplicateurs on change la ligne 2 avec la ligne e 0 0 6 18 12 1 0 M2 E2 M1 E1 A = 0 8 16 . 1 0 0 6 4 1

La stratgie prsente dans cet exemple est appele pivotage partiel . A chaque tape e e e e e on recherche le maximum dans la colonne Ak:n,k . Une alternative est le pivotage complet, ou lon recherche le maximum dans la matrice Ak:n,k:n.
.. . . . . . . . .. .. . . . . . . .. .. .. . . . . . . .. .. .. . . . . . .. .. . . . . . .. .. .. . . . . .. .................... . . ....... .......... . .... . . . . . . . . . . . .. . . . . . . . . . .. .... . . . . . ......... . . ...... . . . . . ..... . .... . . .. ......... . . . . . . .... .... . ......... . . . . .. ... ....... ..... ... . . . . . . . .. . . . . . . .... . . . . . . . . . . .... . . . . . . . . . . .. . . . . .. . . . . . . . .. .. . . . . . . .. .. .. . . . . . . .. .. .. . . . . . .. .. . . . . . .. .. . . . . .. .. .................... . . ................. . .... .............. . . . .. . . . . . . . . . . . . . . .... .............. .... . .. . . . . . . . . .... . . .. . . . . . . . . ...... .... . . . . . . . ...... .... .................. . . . ...... . . .... ............ . .. . . .. . . . . . ... .. . . . . . .. . . . . . . . . . .... .............. . . . . . . .. . . . . . . .

Ak:n,k

Ak:n,k:n

La stratgie du pivotage partiel peut se rsumer par les tapes suivantes : e e e

for k = 1 : n 1 do Dterminer Ek tel que si z est la colonne k de Ek A, on a |zk | = zk:n e Echanger les lignes Appliquer la transformation de Gauss end for

Lalgorithme transforme A en une matrice triangulaire suprieure U en eectuant e les produits : Mn1 En1 M1 E1 A = U . Le pivotage partiel a comme consquence quaucun multiplicateur nest suprieur ` e e a un en valeur absolue. Lalgorithme 8 ralise llimination de Gauss avec pivotage partiel. Tous les multie e plicateurs sont 1 en valeur absolue. A1:k,k est remplac par U1:k,k , k = 1 : n et e Ak+1:n,k par une permutation de Mk k+1:n,k , k = 1 : n 1. Le vecteur e1:n1 dnit e les matrices dchange Ek . Ek change la ligne k avec la ligne ek , k = 1 : n 1. e e Algorithme 8 (egpp) Elimination de Gauss avec pivotage partiel.
1: for k = 1 : n 1 do 2: Chercher p, k p n tel que |Apk | = Ak:n,k 3: Ak,1:n Ap,1:n 4: ek = p 5: Ak:n,k:n = tgm1(Ak:n,k:n ) 6: end for

Lalgorithme 8 remplace la partie triangulaire suprieure de A avec U. La partie e triangulaire infrieure contient les lments de L mais dans le dsordre tant donn e ee e e e les changes de lignes eectus (LU = A). e e Si llimination de Gauss avec pivotage partiel a t applique suivant lalgorithme 8 e ee e pour obtenir la forme triangulaire suprieure e Mn1 En1 Mn2 En2 M1 E1 A = U Mn1 Mn2 M1 En1 En2 E1 A = U
P

on a P A = LU avec P = En1 E1 et L une matrice triangulaire infrieure unitaire avec |ij | 1. e e e La k`me colonne de L correspond au k`me vecteur de la transformation de Gauss qui a t permut, cest-`-dire si Mk = I (k) ek alors Lk+1:n,k = gk+1:n avec ee e a g = En1 Ek+1 (k) .

Exemple 6.10 Si lon applique lalgorithme 8 a la matrice A de lexemple prcdant on ` e e

obtient 6 18 12 1/2 8 16 1/3 1/4 6 et e= 3 3 0 0 1 P =1 0 0 0 1 0 .

Avec Matlab la commande [L,U] = lu(A) retourne une matrice triangulaire infrieure L e permute. On a LU = A. La commande [L,U,P] = lu(A) retourne une matrice triangue laire L mais LU = P A. 1 6 18 12 , U = L = 1/2 1 8 16 et 1/3 1/4 1 6 0 0 1 3 17 10 6 18 12 P A = 1 0 0 2 4 2 = 3 17 10 . 0 1 0 6 18 12 2 4 2

6.4

Considrations pratiques e

Montrons comment rsoudre des syst`mes linaires ` partir dune factorisation de e e e a matrice. Soit la collection de syst`mes linaires AX = B avec A Rnn , none e singuli`re et B Rnq . En posant X = [x1 . . . xq ] et B = [b1 . . . bq ] les q solutions e sobtiennent avec la procdure : e
Calculer P A = LU for k = 1 : q do Rsoudre Ly = P bk e Rsoudre U xk = y e end for

Remarquons que lon consid`re le syst`me linaire permut P AX = P B et la face e e e torisation P A = LU. Cette procdure est programme dans lalgorithme 9. e e Algorithme 9 (rsl) Rsolution dun syst`me dquations linaires AX = B e e e e
1: 2: 3: 4: 5: 6: 7: function B = rsl(A, B) [n,q] = size(B) [A,e] = egpp(A) for k = 1 : q do B1:n,k = fsub1(A,perm(e,B1:n,k)) B1:n,k = bsub(A,B1:n,k) end for

Les fonctions egpp et bsub correspondent aux algorithmes 8 et 4. La fonction fsub1 correspond ` une variante de lalgorithme 3 dans laquelle la matrice L est une matrice a triangulaire unitaire. On remarque que la matrice A nest factorise quune fois. Si B = In alors la solution e correspond ` A1 . a

Considrons encore le syst`me linaire Ak x = b avec A Rnn , b Rn et k entier e e e positif. Un possibilit serait de calculer C = Ak puis de rsoudre Cx = b. Ceci e e 2 3 1 2 1 3 2 3 2 3 ncessite 2(k 1)n + 3 n 2 n + 2n = (2k 2 )n + 2 n ops. Pour montrer e que lon peut procder plus ecacement considrons la solution du syst`me linaire e e e e 3 A x = b. On peut alors crire e A AAx = b
b1

A Ax
b2

= b1

Ax = b2 ce qui se formalise avec lalgorithme suivant :

Calculer P A = LU for i = 1 : k do Remplacer b avec la solution de Ly = P b Remplacer b avec la solution de U x = b end for

Dans ce cas le nombre de ops est 2 n3 + k2n2 . 3 Finalement montrons avec un exemple comment viter le calcul de linverse lorsquil e 1 sagit dvaluer une expression du type s = q A r. On remarque que A1 r est la e solution du syst`me linaire Ax = r. On rsout alors ce syst`me comme montr e e e e e avant, puis on value s = q x. Cet exemple montre quil faut toujours raisonner en e termes dquations linaires ` rsoudre et non en termes dinverse dun matrice. e e a e Revenons ` la remarque faite ` la page 23 qui disait quil ne fallait jamais recourir a a ` linverse pour rsoudre un syst`me linaire. a e e e En ne considrant que les termes dordre le plus lev dans le comptage du nombre e e e doprations lmentaires, rappelons que la factorisation LU ncessite 2 n3 oprations e ee e e 3 2 lmentaires et la solution dun syst`me triangulaire n oprations lmentaires. ee e e ee Ainsi le nombre doprations lmentaires pour la rsolution dun syst`me linaire e ee e e e 2 3 avec LU est dordre 3 n .
2 Linverse A1 sobtient en rsolvant n syst`mes linaires do` elle ncessite 3 n3 + e e e u e 8 n2n2 = 3 n3 oprations lmentaires et le produit A1 b ncessite n2 oprations e ee e e supplmentaires. Donc cette approche est quelque trois fois plus co teuse quune e u solution avec LU. On vrie aussi que la situation o` lon a Ax = B avec un grand e u nombre de colonnes dans B ne modie pas les conclusions qui prc`dent. e e

Illustrons encore que la factorisation produit des rsultats plus prcis. Considrons e e e 12 le syst`me 3x = 12 compos dune seule quation et la solution x = 3 = 4 obtenue e e e par une transformation qui est quivalente ` ce que produit la factorisation. Si lon e a proc`de ` une inversion explicite en utilisant une arithmtique ` trois digits on a e a e a 1 x = 3 12 = .333 12 = 3.99. On voit que linversion explicite ncessite une e opration de plus pour aboutir ` un rsultat moins prcis. e a e e

6.5

Elimination de Gauss-Jordan

Llimination de Gauss visait la transformation du syst`me original en un syst`me e e e triangulaire an quil puisse tre rsolu plus facilement. e e Llimination de Gauss-Jordan est une variante qui transforme le syst`me original e e en un syst`me diagonal en annulant galement les lments qui se trouvent en dessus e e ee de la diagonale. La matrice de transformation pour un vecteur x donn est de la e forme
1 . . .. . . 0 0 0 . .. . . . 0 0 . . . 1 0 0 . . . 0 1 . . .
(k)

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

k1 1 (k) k+1 . . . n
(k)

(k)

0 x1 . . . . . . 0 xk1 0 xk 0 xk+1 . . . . . . xn 1

0 . . . 0 = xk 0 . . . 0

o` i u

(k)

= xi /xk , i = 1, . . . , n.

La transformation de Gauss-Jordan est environ 50% plus coteuse que lliminau e tion de Gauss. Cependant lorsquon applique la transformation aussi au vecteur b la solution du syst`me diagonal ncessite que n divisions. La mthode reste cepene e e dant globalement plus co teuse. Un autre dsavantage est quelle ne garantit pas la u e stabilit numrique, mme lorsquon pivote. e e e Un syst`me diagonal peut cependant tre facilement rsolu en parall`le ce qui ace e e e tuellement conf`re ` cette mthode un certain regain dintrt. e a e ee

Chapitre 7 Matrices particuli`res e

Lorsquon eectue des calculs numriques il est tr`s important de tenir compte de la e e structure du probl`me an de choisir lalgorithme le plus adapt pour sa resolution. e e On va notamment montrer que dans le cas de matrices symtriques, dnies positives e e ou matrices par bande, il existe des variantes de la factorisation LU qui, pour de telles matrices, sont plus ecaces. Dans ce chapitre on prsentera dabord une formalisation dirente de la factorisae e tion LU sous la forme dun produit de trois matrices LDM . Cette factorisation na pas dintrt pratique mais elle sera utile pour obtenir la factorisation LDL dune ee matrice symtrique. Finalement on discutera la situation tr`s importante dune mae e trice dnie positive pour laquelle il existe une factorisation de la forme GG appele e e factorisation de Cholesky.

7.1

Factorisation LDM

Montrons dabord quil est possible de factoriser la matrice A en un produit de trois matrices A = LDM avec D une matrice diagonale et L et M deux matrices triangulaires infrieures e unitaires. Thor`me 7.1 Si toutes les sous-matrices diagonales de A sont non-singuli`res, e e e alors il existe deux matrices triangulaires infrieures unitaires L et M et une matrice e diagonale D tel que A = LDM . Cette factorisation est unique. Dmonstration 1 A partir du thor`me 6.1 on sait que A admet la factorisation A = e e e
LU . Posons D = diag(d1 , . . . , dn ) avec di = uii , i = 1, . . . , n, alors D est non-singuli`re et e M = D1 U est triangulaire suprieure unitaire. Ainsi A = LU = LD(D1 U ) = LDM . e Lunicit dcoule de lunicit de la factorisation LU. e e e

La factorisation LDM peut donc tre obtenu ` partir de la factorisation LU. Il existe e a cependant un algorithme qui permet dobtenir L, D et M directement. Supposons que lon connaisse les j 1 premi`res colonnes de la factorisation A = LDM : e

1:j1

j+1:n

. . .. . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . .... . . . . . . . . . . ..... . .. . . . . . . . . ..... . . . . . . .. . . . . . . . . . . . . .. . . . . . .. . . . . . . . . . . . .. . . . . . . .

1:j1

=
j+1:n

... . . . . . ... . . . . . . . . . . . . . . . .. . . ... .. . . . . . .. . . ... . . .. . . . ... . . . .... . . . .. . . . .. . . .. . . . . . ... . . ... 1 . . . .. . . . . . . ... . .... . . . . . . .. . . ................... . . .. . . . ....... ........ . ................. . ................... . . ... . . .. . . . . . . . . . . . . . ..... ... .. . . . . . . . . . . . . . . . ..... ... . .. . .. .. . . . . . . . . . . . ..... . .. . .. . . . . . . .... .. . .. .. . . . . . .. . . .. . . . . . ..... .. . ... . .. . . .. .. .. . . . . . . . ....... . . . .............. .. ................... . .. . ...................................... . . . .. .. j ....... Lj+1:n,j

. . . . . . . . . . . . . . . ...... dj .. ....... . ... .. .. .. .. .

. .. Mj,1:j1 j ...... ..................................... .................................... . .. .. . . ..... . . ...... .... . . . . . .... . .. ... . . . . . . . . . .... . ..... .. . . ... ..... . 1:j1 . .. . .. . . . . ... .... . .. . . . . ....... . . .. .. . . . . .. .. . . . . . . . . . . . .. . . . . . . . . .. . .. . . . . . . . . . . . . . . .. ... .. . . .. . . . . .. . .. . . .. . . . . 1.. . . ....... . . . . . .. . . . . .. . . .. . . . .. . .. . . . . .. . . . . . .. . .. . . . . ... . .. . . . . . . . ... . . . . . . . . . .

e Dveloppons la j `me colonne de A = LDM e

A1:n,j = Lv

avec v = DM ej .

An didentier les inconnues Lj+1:n,j , dj et Mj,1:j1 on rsoud A1:n,j = Lv en deux e tapes en considrant dabord les j premi`res quations e e e e
A1:j,j .....
.. .. .. .. .. .. .. . .. . . . .... .. . . . .. L1:j,1;j .. v1:j1 .. .. .. .. .. .. .. .. .. . ... ... .. ... . ........ . . . . . .. . .... . . ... . . ... . .. . . . . . . . .. . ... ... .. . . . . .. .. . . . .... . . . .. . .. . .. . . . . . . ... ... 1 . . . . ....... .... . ................ ... . . . . . . .. . .... ... . . . . .............. .... .. ............... .... . ... ... . . . ... ... . . . .. v .. j . . . . . . . . . . . . . . . . . . . . . . . . . . . . ... .

Les lments v1:j sont la solution du syst`me triangulaire ee e L1:j,1:j v1:j = A1:j,j . Etant donn que les di , pour i = 1, . . . , j 1 sont connus et mjj = 1 on peut calculer e mj,i = vi /di dj = vj . i = 1, . . . , j 1

Dans une deuxi`me tape on consid`re les n j equations restantes e e e

. .. . . . . . . . . . . . . . .v . . . . ... 1:j1 ... . .. . ... . ..... . .. . ... . . .. . Aj+1:n,j ..... ... ... . . . ........ . .. . . .. . .... ... . . .. ... ... .. .................. . . .. ... ... . .............. .. .. . .. . . . . ..... ... ... . . . .. . . . ... ... . . . . . . . ..... ... v . . . . .. . j . . . . .... . . . . . . . . . . . ..... . . . . . .. . . . . . . . ..... . ... . . . ... .. . . ... . . .... . . . . .. . . .. . ... . . ...... . ................. .. ... . .. .................... ... .. .. .. .. .. .. .. .. .. .. .L . Lj+1:n,1:j1 . j+1:n,j

qui permettent didentier le vecteur Lj+1:n,j . Ce vecteur est dni par le syst`me e e Lj+1:n,1:j v1:j = Aj+1:n,j

que lon rcrit sous la forme ee Lj+1:n,j vj = Aj+1:n,j Lj+1:n,1:j1 v1:j1 do` lon obtient directement Lj+1:n,j en divisant le membre de droite par vj . La u factorisation se rsume alors par lalgorithme suivant :1 e
d1 = A11 L2:n,1 = A2:n,1 /A11 for j = 2 : n 1 do solve L1:j,1:j v1:j = A1:j,j (Syst`me triangulaire) e Mj,1:j1 = v1:j1 ./d1:j1 dj = vj Lj+1:n,j = (Aj+1:n,j Lj+1:n,1:j1 v1:j1 )/vj end for solve L1:n,1:n v1:n = A1:n,n (Syst`me triangulaire) e Mn,1:n1 = v1:n1 ./d1:n1 dn = vn

Comme pour la factorisation LU, il est possible de formuler lalgorithme de sorte que A soit remplac par L, D et M. e Soit A Rnn et admettant une factorisation LU, alors lalgorithme 10 produit la factorisation A = LDM et remplace A par L, D et M. On a aij = lij pour i > j, aij = di pour i = j et aij = mji pour i < j. Algorithme 10 (ldm) Factorisation A = LDM .
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: A2:n,1 = A2:n,1 /A11 for j = 2 : n 1 do v1:j = fsub1(A1:j,1:j , A1:j,j ) A1:j1,j = v1:j1 ./diag(A1:j1,1:j1 ) Ajj = vj Aj+1:n,j = (Aj+1:n,j Aj+1:n,1:j1 v1:j1 )/vj end for v1:n = fsub1(A1:n,1:n, A1:n,n ) A1:n1,n = v1:n1 ./diag(A1:n1,1:n1 ) An,n = vn

Les instructions 1 et 6 remplacent A par L, en 4 et 9 on remplace A par M et en 5 et 10 cest v qui remplace A. En 3 et 8 on rsout un syst`me triangulaire unitaire. e e Lalgorithme 10 ncessite 2 n3 + 2n2 + e 3 lalgorithme de factorisation LU.
L, D et M .
1

10 n 3

9 ops. Il est de mme complexit que e e

Exemple 7.1 Factorisation dune matrice A en LDM . Lalgorithme remplace A avec

3 3 6 1 0 0 3 0 0 1 1 2 3 1 2 A = 9 15 42 = 3 1 0 0 6 0 0 1 4 A = 3 6 4 12 54 194 4 7 1 0 0 2 0 0 1 4 7 2

Cette version spare les tapes pour j = 1 et j = n de ltape gnrique. Voir les fonctions e e e e e ldm1, ldm2, ldm3 et ldm4.

On pourra vrier que lalgorithme eg conduit a la mme matrice L et que U = DM . e ` e

7.2

Factorisation LDL

Lorsque la matrice A est symtrique la factorisation LDM comporte des oprations e e rdondantes. e Thor`me 7.2 Si A = LDM est la factorisation dune matrice non-singuli`re et e e e symtrique A, alors L = M. e Dmonstration 2 La matrice M 1 A(M )1 = M 1 LD est a la fois symtrique et e ` e

triangulaire infrieure, donc diagonale. Comme D nest pas singuli`re M 1 L est aussi e e 1 L est aussi une matrice triangulaire unitaire et donc M 1 L = I. diagonale. Mais M

Les premi`res j 1 composantes du vecteur v peuvent donc tre obtenues par des e e simples produits terme ` terme et vj est obtenue ` partir de lquation L1:j,1:j v = a a e A1:j,j , cest-`-dire a vj = Ajj Lj,1:j1 v1:j1 ce qui donne lalgorithme :
for j = 1 : n do for i = 1 : j 1 do vi = Lj,i di end for vj = Ajj Lj,1:j1 v1:j1 dj = vj Lj+1:n,j = (Aj+1:n,j Lj+1:n,1:j1 v1:j1 )/vj end for

Ainsi il sera possible dconomiser la moiti des oprations dans lalgorithme 10 si la e e e matrice A est symtrique. En eet, ` ltape j on conna dj` Mj,1:j1 tant donn e a e t e a e e que M = L. Prcdemment le vecteur v tait dni comme v = DM ej et comme e e e e M = L le vecteur v devient d1 Lj,1 . . . v= . dj1 Lj,j1 dj

A nouveau il est possible de remplacer la matrice A avec le rsultat. Soit A e nn R , symtrique et admettant la factorisation LU alors lalgorithme 11 produit la e factorisation A = LDL et remplace A. On a aij = lij pour i > j et aij = di pour i = j. En 6 on remplace A avec D et en 7 A est remplace par L. Lalgorithme 11 ncessite e n3 ops, cest-`-dire la moiti de ce que ncessite la factorisation LU. a e e 3

Algorithme 11 (ldl) Factorisation A = LDL.

1: for j = 1 : n do 2: for i = 1 : j 1 do 3: vi = Aj,i Ai,i 4: end for 5: Ajj = Ajj Aj,1:j1 v1:j1 6: Aj+1:n,j = (Aj+1:n,j Aj+1:n,1:j1 v1:j1 )/vj 7: end for

Exemple 7.2 Factorisation dune matrice symtrique A en LDL . Lalgorithme reme

place A par 10 A = 20 30 L et D. 20 30 1 0 0 10 0 0 1 2 3 10 20 30 45 80 = 2 1 0 0 5 0 0 1 4 A = 2 5 80 80 171 3 4 1 0 0 1 0 0 1 3 4 1

7.3

Matrices symtriques dnies positives e e

Les syst`mes linaires Ax = b avec une matrice A dnie positive constituent une e e e classe, parmi les plus importantes, des syst`mes linaires particuliers. Rappelons e e nn quune matrice A R est dnie positive si quelque soit x Rn , x = 0 on a e x Ax > 0. a a Considrons le cas symtrique avec A= 11 12 . Si A est dnie positive on a e e e a12 a22 x Ax = a11 x2 + 2a12 x1 x2 + a22 x2 > 0 ce qui implique pour 1 2
x = [1 x = [0 x = [1 x = [1 0] 1] 1] -1] a11 > 0 a22 > 0 a11 + 2a12 + a22 > 0 a11 2a12 + a22 > 0

Les deux derniers rsultats impliquent que |a12 | (a11 + a22 )/2 et que le plus e grand lment de A se trouve sur la diagonale et quil est positifs. Il en rsulte ee e quune matrice symtrique et dnie positive a des lments diagonaux qui sont e e ee relativement grands ce qui rendra le pivotage inutile.
e e e Proprit 7.1 Soit A Rnn et dnie positive, alors A est non-singuli`re. Dmonstration : e e Car sinon il existerait x = 0 tel que x Ax = 0 et donc Ax = 0 ce qui nest pas possible pour A non-singuli`re. e

Proprit 7.2 Si A Rnn et dnie positive et si X Rnk est de rang k alors e e e B = X AX Rkk est aussi dnie positive. e Dmonstration 3 Montrons que le contraire nest pas possible. Soit un vecteur z = 0 e
et z Rk qui satisfait 0 z Bz = (Xz) A(Xz) alors Xz = 0. Mais, le fait que les colonnes de X sont linairement indpendantes implique z = 0. e e

Corollaire 1 Si A est dnie positive, alors toutes les sous-matrices diagonales e sont dnies positives. En particulier tous les lments diagonaux sont positifs. e ee Dmonstration 4 Soit v Rk un vecteur entier avec 1 v1 < . . . < vk n, alors e

X = In (:, v) est une matrice de rang k compose des colonnes v1 , . . . , vk de la matrice e didentit. A partir de la proprit 7.2 il sensuit que Av,v = X AX est dnie positive. e ee e

Corollaire 2 Si A est dnie positive, alors la factorisation A = LDM existe et e tous les lments de D sont positifs. ee Dmonstration 5 A partir du premier corollaire il dcoule que les sous-matrices A1:k,1:k e e sont non-singuli`res pour k = 1, . . . , n. A partir du thor`me 7.1 il dcoule que la face e e e torisation A = LDM existe. Si lon applique la proprit 7.2 avec X = (L1 ) alors ee B = L1 LDM (L1 ) = DM (L1 ) = L1 A(L1 ) est dnie positive. Comme M (L1 ) e est triangulaire suprieure unitaire, B et D ont la mme diagonale, qui doit tre positive. e e e Une situation pratique courante ou lon rencontre des matrices dnies positives est e donne lorsque A = X X avec toutes les colonnes de X linairement indpendantes. e e e (c.f. proprit 7.2 avec A = In ). ee

7.4

Factorisation de Cholesky

Nous savons que pour des matrices symtriques il existe la factorisation A = LDL . e Montrons quil existe encore une autre factorisation. Thor`me 7.3 (Factorisation de Cholesky) Si A Rnn est symtrique et dnie e e e e positive, alors il existe une matrice triangulaire infrieure unique G Rnn avec les e lments diagonaux positifs et telle que A = GG . ee
e e Dmonstration 6 A partir du thor`me 7.2 on sait que la factorisation A = LDL e

existe. Comme les dk sont positifs, la matrice G = L diag( d1 , . . . , dn ) est une matrice relle triangulaire infrieure avec diagonale positive. On a A = GG . Lunicit dcoule de e e e e lunicit de la factorisation LDL . e

La matrice G est appelle le triangle de Cholesky. La dmonstration du thor`me e e e fournit une mthode pour obtenir G. Il existe cependant une mthode plus ecace e e pour obtenir G. Considrons ltape de la factorisation A = GG ou les premi`res j 1 colonnes de e e e G sont connus :

1:j1

j+1:n

.. . . . . . . . . . . . .. . . . . . .. . . . . . . .. . . . . . .. . .. . . . . . . . ..... . . .. . .. .. . . . . . . . ..... . . . . . . .. . . . . . . . .. . . . . . .. . . . . . .. . . . . . . . . . . . .. . . . . . . .

1:j1

=
j+1:n

... . . . . . . . . . . . . . . . . . . . .. . . .. ..... . . . . .. . . . ... . . . .. . . . . . .... .. . . . . . . ... . . . . .. . . .. . ................... .... Gjj . . . . . . . . . . . ...... .. .. . . ....... ............ . . . .. . . . .... . ....... .......... .................... ... . . .. . . . . ... . . . . . . . . . . . ... . ... . . .. . . .. . . . . . . . . ... . ... ... . . .. . . . .... .. . . . . . . . . . ... . ... . .. . . .. . . .. . .. ... . . . . . . . . . . ... . . . .. . .. . .... .. . . . . . . . .... ... .. . . .. . .... . . .. . . . . . . . . . . ... ... ... . . . .............. . . .................. . .. . . ..................................... . . .. . .. . . j ....... Gj:n,j

.. .. Gj,1:j1 j ...... . .. . . . .. ...................................... .......... ....... .................. . . . .. .... . . . . . . . ... . .. . . . . . .. . . . . ... . . . . ... . . . .. .. . . . . . ... . . . . 1:j1 . .. . . . . . .. . . . . . . .. . . . . . .. . . . .. . . . . .... . . . . . . . . . . . . .. . . . . .. . . .. . . . . .. ... . . . . . . . . . .. . . .. . .. . . .. . .. . . . .. .. . . . .. . .. . . . .. .. . . .. . .. . . .. . . .. . . .. . .. . . .. . ...................

e La j`me colonne du produit sexprime alors sous la forme dune combinaison des j colonnes de la matrice G1:n,1:j , soit :
j

A1:n,j =
k=1

Gjk G1:n,k

dou lon peut expliciter

Gjj G1:n,j = A1:n,j

k=1

Gjk G1:n,k v

et comme G1:n,1:j1 est connu, on peut calculer v. Considrons alors lexpression e Gjj G1:n,j = v. Les composantes du vecteur G1:j1,j sont nulles car G est triangulaire e infrieure. La j`me composante scrit Gjj Gjj = vj dou lon tire Gjj = vj . La e e solution du vecteur Gj:n,j scrit alors : e Gj:n,j = vj:n / vj ce qui conduit ` formuler lalgorithme : a
for j = 1 : n do vj:n = Aj:n,j for k = 1 : j 1 do vj:n = vj:n Gjk Gj:n,k end for Gj:n,j = vj:n / vj end for

Il est possible de formuler lalgorithme de sorte ` remplacer le triangle infrieur a e nn de la matrice A par la matrice G. Soit A R , symtrique et dnie positive, e e alors lalgorithme 12 calcule une matrice triangulaire infrieure G Rnn telle que e A = GG . Pour i j, lalgorithme remplace Aij par Gij . Algorithme 12 (cholesky) Factorisation de Cholesky A = GG .
1: for j = 1 : n do 2: if j > 1 then 3: Aj:n,j = Aj:n,j Aj:n,1:j1 A j,1:j1 4: end if 5: Aj:n,j = Aj:n,j / Ajj 6: end for

Lalgorithme ncessite n + n2 + 8 n 2 ops. Ceci correspond ` la moiti de ce que e a e 3 3 ncessite la factorisation LU. Il nest pas ncessaire de pivoter. e e Dtail pour le calcul des ops e
Instruction 5: 3: 3: 3: 3: 5: 3+5 : Produit Addition Indices Total Opration e Nombre doprations e j=1 n+1 j>1 2(n j + 1)(j 1) nj +1 2 2(n j + 1)(j 1) + n j + 1 + 2 = (n j + 1) 2(j 1) + 1 + 2 (n j + 1) + 1 (n j + 1) 2(j 1) + 2 + 3

La somme des oprations lmentaires est n + 1 + n e ee j=2 (n j + 1) 2(j 1) + 2 + 3 . Avec Maple on peut valuer cette expression a laide des commandes : e ` > fl := n +1 + sum ( (n - j + 1 ) * ( 2 * (j - 1 ) + 2 ) + 3 , j =2.. n ); > s i m p l i f y( fl );

Remarques
Avec Matlab on obtient la factorisation de Cholesky avec la commande R=chol(A) ou R est une matrice triangulaire suprieure dou R R = A. e La factorisation de Cholesky sutilise ecacement pour vrier numriquement si e e une matrice est dnie positive. Pour ce faire il sut de vrier qua chaque tape e e e de lalgorithme ajj est positif. Dans le cas contraire la matrice nest pas dnie e positive. Avec Matlab on peut appeler la factorisation de Cholesky avec un deuxi`me are gument [R,p] = chol(A) qui indique llment diagonal qui est nul. Si p est nul ee la matrice est dnie positive sinon seul la sous-matrice form des p 1 premieres e e lignes et colonnes est dnie positive. e La dcomposition de Cholesky est aussi utilise pour la gnration de variables e e e e alatoires multidimensionelles x vriant une variance et covariance V donne. On e e e proc`de comme : e x = Ge o` G est la factorisation de Cholesky de la matrice V et e N(0, I). La matrice de u variance et covariance de x vrie alors e E(xx ) = E(G ee G ) = GG = V .
I

La factorisation de Cholesky est aussi frquemment utilise pour la solution des e e quations normales. e

Chapitre 8 Matrices creuses

8.1 Matrices par bande

Une matrice dont tout lment aij pour |i j| > b est nul est appele une matrice ee e par bande avec b la largeur de bande. Pour de telles matrices les techniques de factorisation vues avant peuvent tre modis de faon a ne pas faire intervenir les e e c ` lments nuls. ee Un cas particulier est constitu par une matrice tridiagonale pour laquelle on a e b = 1. Dans la situation o` il nest pas ncessaire de pivoter pour assurer la stabilit u e e numrique, la factorisation LU = A scrit : e e
1 l1 1 l2 1 l3 1 l4 1 u1 r1 d1 q1 p1 d2 q2 u2 r2 = u3 r3 p2 d3 q3 u4 r4 p3 d4 q4 u5 p4 d5

(8.1)

Lalgorithme 13 eectue une factorisation A = LU dune matrice tridiagonale A avec pi , i = 1, . . . , n 1 la diagonale infrieure, di , i = 1, . . . , n la diagonale et qi , e i = 1, . . . , n 1 la diagonale suprieure. On vrie que ri = qi , i = 1, . . . , n 1. e e Algorithme 13 (lu3diag) Factorisation LU dune matrice tridiagonale A.
1: u1 = d1 2: for i = 2 : n do 3: li1 = pi1 /ui1 4: ui = di li1 qi1 5: end for

Lalgorithme 13 ncessite 4 n 4 oprations lmentaires. e e ee 55

8.2

Matrices creuses irrguli`res e e

Dnition dune matrice creuse (J. H. Wilkinson) : Toute matrice avec susamment e de zros pour quil y ait un gain ` les considrer. Le gain consiste ici ` viter e a e a e des oprations redondantes avec des lments nuls et surtout dviter ` devoir les e ee e a conserver en mmoire. e Dans la pratique on rencontre souvent des syst`mes creux de grande taille. Sans e traitement particulier ces probl`mes resteraient impraticables. e Exemple 8.1 La discrtisation du Laplacien a 5 points avec une grille de 64 64 noeuds e `
conduit a une matrice dordre 4096 4096 avec 20224 lments non nuls. Le tableau qui ` ee suit indique la mmoire requise pour conserver la matrice et le temps ncessaire pour e e eectuer un produit et la solution dun syst`me linaire lorsquon recourt a la technique e e ` traditionnelle et lorsquon exploite la structure creuse.

Mmoire e Dx Dx = b

Plein 128 Mbytes 30 sec > 12 h

Creux 0.25 Mbytes 0.2 sec 10.0 sec

Lorsquon exploite ecacement la structure creuse dune matrice comportant nnz lments non nuls la complexit dalgorithmes faisant intervenir des produits et des ee e resolutions de syst`mes linaires est O(nnz ) et ne fait pas intervenir lordre n des e e matrices.

8.2.1

Reprsentation de matrices creuses e

Dans la suite nous mettrons laccent sur lapproche ralise dans Matlab pour le traie e tement de structures creuses. Pour lutilisateur de Matlab le maniement de structures creuses est presque transparent. On peut envisager de nombreuses faons pour reprsenter une matrice creuse et il c e nexiste pas de schma qui soit le meilleur dans labsolu. Le choix optimal dpend e e du type doprations eectues sur les lments de la matrice. e e ee Matlab stocke les matrices quelles soient pleines ou creuses colonne par colonne. Une matrice creuse est reprsente par la concatenation des lments non nuls des e e ee vecteurs colonnes. Soit la matrice 0 1 0 A= 2 0 3 4 0 0

La commande nnz renseigne sur la nombre dlments non nuls et nonzeros empile ee les lments non-nuls colonne apr`s colonne. ee e
>> nnz ( A ) ans = 4 >> nonzeros ( A ) ans = 2 4 1 3

Linformation est stocke colonne par colonne e 0 2 4 1 0 0 0 3 0

Avec la commande find on obtient le vecteur dindices des lments non nuls dans ee le vecteur des colonnes empiles. e
>> find ( A ) ans = 2 3 4 8

Pour obtenir lindice de la ligne et de la colonne des lments non nuls et leur valeur ee on excute la commande e
>> [i ,j , e ] = find ( A ) i = 2 3 1 2 j = 1 1 2 3

e = 2 4 1 3

Dans Matlab une matrice creuse est reprsente par la liste des lments non nuls e e ee dans chaque colonne. Cette reprsentation ncessite un vecteur de pointeurs avec e e autant dlments que la matrice a de colonnes, une liste qui contient les indices de ee lignes correspondant aux nnz lments non nuls et un vecteur de mme longueur ee e pour conserver la valeur des lments. ee Ci-apr`s cette reprsentation pour la matrice creuse A de lexemple prcdent o` e e e e u n = 3 et nnz = 4.
Colonnes
1

Pointeurs :

1
. . . . . . . . . . . . . .. .. . .

. . . . . . . . . . . .

. . . . . . . . . . .. . . .. . ... . ... . ... ... ........ ..... ... ... ... .. ... ... .. . . . . . . . . . . . . . . . .. .. . .. . . .. . . . . . . .

. . . . . . . . . . . .

n(+1)

Indices des lignes : e Elments :

2 2

. . . . . . . . . . . . . . . . . . . . . .

3 4
. . . . . . . . . . .

1 1
. . . . . . . . . . .

2 3

nnz nnz

Cette reprsentation ncessite n(+1) mots entiers pour le vecteur de pointeurs, nnz mots e e entiers pour les listes des indices de lignes et nnz mots rels pour les lments. Comme e ee un entier ncessite 4 bytes et un rel 8 bytes, le total des bytes ncessaires a cette e e e ` reprsentation est 12 nnz + 4 n bytes. e

8.2.2

Conversion entre reprsentation pleine et creuse e

Contrairement a ce qui est le cas pour la conversion de matrices relles en matrices e complexes, la conversion en matrice creuse ne se fait pas automatiquement. Lutilisateur doit donner une commande explicite. Cependant une fois initi, la reprsentation e e creuse se propage, cest-`-dire que le rsultat dune opration avec une matrice creuse a e e produit une matrice creuse. (Ceci nest pas le cas pour une addition ou soustraction.) Les fonction full et sparse eectuent la conversion. Conversion de la matrice A dni avant : e
>> B = sparse ( A ) B = (2 ,1) 2 (3 ,1) 4 (1 ,2) 1 (2 ,3) 3

>> C = full ( B ) C = 0 1 0 2 0 3 4 0 0 >> whos Name Size A B C 3 x3 3 x3 3 x3

Bytes 72 64 72

Class double array double array ( sparse ) double array

8.2.3

Initialisation de matrices creuses

Il est possible et indiqu de crer directement une matrice creuse (sans recourir ` e e a la conversion dune matrice pleine) en donnant une liste dlments et les indices ee correspondants. La commande
S = sparse (i ,j ,s ,m ,n , nzmax );

cre une matrice creuse avec Si(k),j(k) = s(k). Les lments dans les vecteurs dine ee dices i et j peuvent tre donnes dans un ordre quelconque. Remarque importante : e e Si une mme paire dindices existe plus dune fois les lments correspondants sont e ee additionns. Si s ou lun des arguments i, j est un scalaire ils sajustent automatie quement.
S S S S = = = = sparse (i ,j ,s ,m , n ); % nzmax = length ( s ) sparse (i ,j , s ); % m = max ( i ) n = max ( j ) sparse (m , n ); % S = sparse ([] ,[] ,[] ,n , m ) sparse ([2 3 1 2] ,[1 1 2 3] ,[2 4 1 3]) % matrice A

8.2.4

Principales fonctions pour matrices creuses

help sparfun speye (n , m ) sprandn (n ,m , d ) spdiags spalloc (n ,m , nzmax ) issparse ( S ) spfun % C = spfun ( exp , B ) spy ( S ) [p ,q ,r , s ] = dmperm ( S ) condest ( S ) eigs sprank

8.3

Proprits structurelles des matrices creuses e e

Pour des matrices creuses on peut mettre en vidence des proprits intressantes e ee e qui dpendent uniquement de la structure qualitative de la matrice, cest-`-dire e a de linformation indiquant quels sont les lments non nuls. Lors de la rsolution de ee e syst`mes linaires ou non-linaires les proprits structurelles des matrices dnissant e e e ee e ces syst`mes sont tr`s importantes. e e Pour tudier de telles proprits il convient alors dassocier ` une matrice sa matrice e ee a dincidence. La matrice dincidence M dune matrice A est dnie comme e mij = 1 si aij = 0 . 0 sinon

Dans la suite les proprits structurelles dune matrice seront tudis en analysant ee e e sa matrice dincidence1 .

Rang structurel dune matrice Rappelons que le dterminant dune matrice M Rnn peut sexprimer comme e det M =
pP

sign(p) m1p1 m2p2 mnpn

(8.2)

o` P est lensemble des n! permutations des lments de lensemble {1, 2, . . . , n} et u ee sign(p) est une fonction qui prend des valeurs +1 et 1. De la relation (8.2) on dduit quune condition ncessaire pour la non-singularit e e e de la matrice M est lexistence dau moins un produit non nul dans la somme de lexpression dnissant le dterminant. Ceci correspond ` lexistence de n lments e e a ee non nuls, chaque ligne et chaque colonne de M contenant exactement un de ces lments. ee Une formulation quivalente de cette condition consiste a dire quil doit exister une e ` permutation des colonnes de la matrice M de sorte que la diagonale de la matrice ne comporte pas dlments nuls. On parle alors dexistence dune normalisation. ee Si la matrice correspond ` la matrice jacobienne dun syst`me dquations, lexisa e e tence dune normalisation implique quil est possible dcrire les quations de sorte e e que chaque quation explicite une variable endog`ne dirente. e e e Soit la matrice dincidence dune matrice creuse pour laquelle on donne galement e
La fonction spones produit la matrice dincidence et la fonction spy produit le graphique de la matrice dincidence.
1

le graphe biparti associ. e h1 h2 M = h3 h4 h5 1 1 . 1 1 1 1 . . . . 1 1 . 1 1 1 . 1 . 1 . 1 . . H Y h1 c................................................................................................................................................................................. c y1

........ . .. ... ....... . . .. . .. . ... ... ....... h2 c............................................................................................................................................................................................................ c y2 . . ...... ... ...... ...... ... ...... .. . . .. . . h3 c..................................................................................................................................................................................................... c y3 . ..... . ..... .. .... ... . .. .... ... .... .. ... .... .. h4 c....................................................................................................................................................... c y4 ...... ...... h5 c.. . ... . .. .... .. ... .... .. .... .. .... .. .... . .... ... . . .. .. .. ... .. ... ...... . ... .. .. .... ... .. ........ ... . ... .. . .. .. . .... .. ...... .. . .... .. ....... . .. . ... ....... .. . ... ... ... ..... . . ... ..... . .. . ..... .. . . . ....... ....... ..... . ...... ....... .. ..

y1 y2 y3 y4 y5

On vrie facilement que les lments de lensemble {m51 , m22 , m33 , m44 , m15 } dnissent e ee e une normalisation et que de ce fait, la matrice M ci-dessus vrie la condition e ncessaire pour la non-singularit. e e Dans le graphe biparti associ ` la matrice une normalisation correspond ` un enea a semble dartes non adjacentes qui relient tous les sommets de lensemble H avec e tous les sommets de lensemble Y . Un ensemble dartes est appel un couplage et e e on le note W . Si le couplage couvre tous les sommets du graphe on parle de de couplage parfait. La fonction dmperm2 de Matlab permet la recherche dun couplage W de cardinalit e maximale. La commande
p = dmperm ( M )

produit un vecteur dont les lments sont dnis comme ee e pi = j si mji W . 0 sinon

Si pi = 0, i, i = 1, . . . , n la permutation M(p, :) produit une matrice avec diagonale non nulle. On appelle rang structurel la cardinalit du couplage maximal, ce qui correspond e au nombre dlments non nuls du vecteur p. La fonction Matlab sprank donne ce ee nombre3 . Matrices structurellement singuli`res e Le rang structurel dune matrice M vrie toujours e sprank(M) rank(M). Soit une matrice carre dordre n de rang structurel k avec k < n. Dans une telle e situation on peut vouloir savoir ou introduire des lments nouveaux dans la maee trice an daugmenter son rang structurel. An dillustrer le probl`me considrons e e
2 3

Cette fonction ralise lalgorithme de Dulmage and Mendelsohn (1963). e La fonction sprank calcule ce nombre comme : r = sum(dmperm(M) > 0).

un syst`me de 8 quations et la matrice dincidence de la matrice jacobienne e e associe : e h1 h2 h3 h4 h5 h6 h7 h8 : f1 (y1 , y2 , y3 , y4, y5 ) = 0 : y6 = f2 (y3) : y3 = f3 (y7) : f4 (y1 , y2 , y4 , y6, y7 , y8 ) = 0 : f5 (y5 , y6 ) = 0 : y6 = f6 (y7) : y7 = f7 (y3) : f8 (y3 , y5 ) = 0
h1 h2 h3 h4 h5 h6 h7 h8 1 . . 1 . . . . 1 . . 1 . . . . 1 1 1 . . . 1 1 1 . . 1 . . . . 1 . . . 1 . . 1 . 1 . 1 1 1 . . . . 1 1 . 1 1 . . . . 1 . . . .

h y

y1 y2 y3 y4 y5 y6 y7 y8

Avec la fonction dmperm on vrie que le rang structurel de la matrice jacobienne e est 6. En eet en excutant p = dmperm(M) on obtient e p= 14205630 . On peut montrer (Gilli, 1995; Gilli and Garbely, 1996) que les lignes et les colonnes dune telle matrice peuvent tre rordonnes de sorte ` faire appara une souse e e a tre matrice nulle de dimension v w avec v + w = 2 n k. Pour augmenter le rang structurel de la matrice, de nouveaux lments non nuls devront alors tre introduits ee e exclusivement dans cette matrice nulle mise en vidence. e On peut obtenir les vecteurs de permutation p et q qui mettent en vidence la matrice e nulle avec la fonction dmperm de la faon suivante : c
[p ,q ,r , s ] = dmperm ( M )

Pour notre exemple on obtient les rsultats : e p= 41785632 q= 84215673 r= 139 s= 159 . Les vecteurs r et s dnissent les lignes et les colonnes de la dcomposition. La e e matrice jacobienne rordonne est e e
h4 h1 h7 h8 h5 h6 h3 h2 1 . . . . . . . 1 1 . . . . . . 1 1 . . . . . . 1 1 . . . . . . . 1 . 1 1 . . . 1 . . . 1 1 . 1 1 . 1 . . 1 1 . . 1 1 1 . . 1 1

y8 y4 y2 y1 y5 y6 y7 y3

Dcomposition triangulaire par blocs dune matrice e Les lignes et les colonnes dune matrice carre et creuse qui admet un couplage pare fait peuvent tre permutes de sorte ` faire appara une structure triangulaire par e e a tre blocs avec des matrices diagonales qui sont carres et indcomposables. Les situae e tions extrmes dune telle dcomposition sont une matrice triangulaire, ce qui correse e pond ` un syst`me dquations rcursives et une matrice indcomposable, ce qui cora e e e e respond ` un syst`me dquations interdpendant. Dans les situations intermdiaires a e e e e on peut mettre en vidence une succession de sous-syst`mes interdpendants. e e e Ci-apr`s une matrice dcomposable M, les commandes Matlab ncessaires pour la e e e mise en vidence de la dcomposition, ainsi que la matrice rordonne. e e e e
1 2 3 4 5 6 7 8 9 10 11 12 13 . . . . . . . . . . 1 . . . 1 . . . . . . 1 . . . . . . . . . . . . . . . 1 1 . . . . . . . . . 1 . . 1 . . 1 1 . . . 1 . . 1 . . . . 1 . . . 1 . 1 . . . . . . . . . 1 1 . . . 1 . . . . . . 1 . . . . . . . . . . 1 . . 1 1 . . . . . . . . . 1 . . . 1 . . . . . 1 . . . . . . . . 1 . 1 . 1 . . . 1 . . . . 1 . . . . . . . . . . . 1 . . . . 11 9 2 5 1 10 13 12 7 6 3 8 4 1 . . . . . . . . . . . . 1 . . . . . . . . . . . 1 1 . . . . . . . . . . . . 1 . . . . . . . . . . . 1 1 1 . . . . . . . . . . . 1 1 . 1 . . . . . . . . . . 1 1 . . . . . . . . . . . 1 1 . . . . . 1 . . . . . . 1 . 1 . 1 . . . . . . . 1 1 . . . . . . . . . . 1 1 1 . . . . . . . . . . . . 1 1 . . . . . . . . . 1 1 5

. .

. . . . .

1 1

1 2 3 4 5 6 7 8 9 10 11 12 13

1 13 2 8 12 11 4 3 6 7 9 10

La commande Matlab
[p ,q ,r , s ] = dmperm ( M )

produit les rsultats suivants e p = 11 9 2 5 1 10 13 12 7 6 3 8 4 q = 1 13 2 8 12 11 4 3 6 7 9 10 5 r = 1 2 3 4 5 9 12 14 s = 1 2 3 4 5 9 12 14 et la permutation de la matrice M ` droite sobtient avec la commande M(p,q). a Lintrt dune telle dcomposition appara lors de la rsolution dun syst`me ee e t e e dquation tant donn que lon pourra rsoudre de mani`re recursive les souse e e e e syst`mes interdpendants correspondant aux blocs diagonaux. e e

8.3.1

Exemples de matrices creuses

Les exemples qui suivent illustrent les gains decacit considrables obtenus en e e considrant la structure creuse de matrices. e

close all , clear n = 400; d = .0043; randn ( seed ,270785650); figure (1) A = sprandn (n ,n , d ) + speye ( n ); B = sprandn (n ,100 ,0.5* d );

subplot (221) , spy ( A ) subplot (222) , spy ( B )

100

200

300

400

100

200 300 nz = 1087

400

0 50 100 nz = 86

Fig. 8.1 Syst`me linaire creux. e e

AF = full ( A ); BF = full ( B ); whos Name A AF B BF d n Size 400 x400 400 x400 400 x100 400 x100 1 x1 1 x1 Bytes 14660 1280000 1436 320000 8 8 Class sparse double sparse double double double

array array array array array array

Grand total is 201176 elements using 1616112 bytes

Lexemple suivant illustre le produit de matrices creuses et la rsolution de syst`mes e e linaires creux. e
flops (0) , tic , Y = AF * BF ; fprintf ( \ n * Full : %5.2 f sec flops (0) , tic , X = A * B ; fprintf ( \ n * Sparse : %5.2 f sec flops (0) , tic , Y = AF \ BF ; %7 i flops , toc , flops );

%7 i flops \ n , toc , flops );

fprintf ( \ n \ n \\

Full : %5.2 f sec

%7 i flops ,toc , flops );

flops (0) , tic , X = A \ B ; fprintf ( \ n \\ Sparse : %5.2 f sec * Full : * Sparse : \ Full : \ Sparse : 0.38 sec 0.00 sec 0.83 sec 0.11 sec

%7 i flops \ n ,toc , flops );

32000018 flops 496 flops 9100260 flops 231664 flops

Dcomposition dune syst`me dquations linaires en syst`mes indcomposables. e e e e e e

subplot (223) , mgspypar ( A )
0

100

200

300

400

100

200 300 nz = 1087

400

Fig. 8.2 Syst`me linaire dcompos par blocs. e e e e

function mgspypar ( A ) % Display the b l o c k t r i a n g u l a r pattern of a matrix n = size (A ,1); [p ,q ,r , s ] = dmperm ( A ); % A (p , q ) block - order bottom up spy ( A ( p ( n : -1:1) , q ( n : -1:1))); hold on ncfc = length ( r ) - 1; if ncfc > 1 for k = ncfc : -1:1 nf = r ( k +1) -1; nd = r ( k ); if nf > nd i = n - nd + 1.5; j = n - nf + 0.5; plot ([ i j ] ,[ j j ] , b - , L i n e W i d t h ,.05); plot ([ i j ] ,[ i i ] , b - ); plot ([ i i ] ,[ j i ] , b - ); plot ([ j j ] ,[ j i ] , b - ); end

end end

Chapitre 9 Mthodes itratives stationnaires e e

Les mthodes de rsolution directes proc`dent dabord a une transformation du e e e ` syst`me original, puis ` la rsolution du syst`me transform. La transformation e a e e e consiste ` factoriser la matrice, puis la rsolution sapplique ` des syst`mes tria e a e angulaires. De ce fait, une mthode directe conduit, abstraction faite des erreurs e darrondi, ` la solution exacte, apr`s un nombre doprations lmentaires connu et a e e ee ni. Les mthodes itratives gn`rent une squence dapproximations {xk }, k = 0, 1, . . . e e e e e de la solution et de ce fait fournissent une rponse approximative en un nombre inni e doprations. Cependant dans la pratique lapproximation est souvent satisfaisante e au bout dun nombre ditrations relativement petit. e Etant donn que la complexit des mthodes directes est O(n3 ), la taille des probl`mes e e e e que lon peut rsoudre avec la technologie actuelle est de lordre de quelques mile liers. Souvent les grands syst`mes dequations sont creux, cest-`-dire le nombre e a 2 dlments non nuls ne constitue quune tr`s petite fraction de n . Dans une telle ee e situation une mthode directe eectue un grand nombre doprations redondantes e e du fait quelles impliquent des coecients dont la valeur est nulle. En adaptant les algorithmes ` la structure creuse dun probl`me on peut rduire la complexit en a e e e vitant notamment ces oprations redondantes. Ainsi il existe des mthodes directes e e e dites creuses qui adaptant les mthodes de factorisation ` la structure creuse de la e a matrice A. La mise en oeuvre de telles mthodes nest cependant pas triviale. e Les mthodes itratives proc`dent simplement ` des produits matrice-vecteur et e e e a il devient alors facile de programmer ce produit de sorte ` viter des oprations a e e redondantes impliquant des lments nuls de la matrice. Les mthodes itratives ee e e peuvent dans certaines circonstances constituer une alternative aux mthodes die rectes creuses. Le fait quelles soient aussi simples ` mettre en oeuvre, tout en exploitant la struca ture creuse du probl`me, a beaucoup contribu ` leur popularit (surtout parmi les e ea e conomistes). Leur ecacit depend cependant de la vitesse avec laquelle la squence e e e

{xk } converge vers la solution.

9.1

Mthodes de Jacobi et de Gauss-Seidel e

Ces mthodes remontent ` Gauss,1 Liouville (1837) et Jacobi (1845). e a Soit une matrice A dont tous les lments diagonaux sont non-nuls. Pour une matrice ee creuse on doit pouvoir mettre en vidence une diagonale non nulle par des permue 2 tation des lignes et des colonnes . On appelle normalisation des quations la mise e en vidence dune telle diagonale et son existence constitue une condition ncessaire e e pour la non-singularit de la matrice. Considrons alors le syst`me linaire Ax = b e e e e dordre 3 ci-apr`s : e a11 x1 + a12 x2 + a13 x3 = b1 a21 x1 + a22 x2 + a23 x3 = b2 a31 x1 + a32 x2 + a33 x3 = b3 Explicitons llment diagonal de chaque ligne et rcrivons le syst`me de la faon ee ee e c suivante : x1 = (b1 a12 x2 a13 x3 )/a11 x2 = (b2 a21 x1 a23 x3 )/a22 x3 = (b3 a31 x1 a32 x2 )/a33 Les xi dans le membre droit des quations ne sont pas connus, mais supposons que e x(k) est une approximation pour x = A1 b, alors on peut imaginer dobtenir une nouvelle approximation en calculant x1 = (b1 a12 x2 a13 x3 )/a11 (k+1) (k) (k) x2 = (b2 a21 x1 a23 x3 )/a22 (k+1) (k) (k) x3 = (b3 a31 x1 a32 x2 )/a33 Ceci dnit litration de Jacobi qui est formalise avec lalgorithme 14. e e e Litration de Jacobi est particuli`re dans le sens quon nutilise pas les rsultats les e e e plus rcents et on proc`de ` n rsolutions unidimensionnelles indpendantes. Ainsi e e a e e (k+1) (k) (k+1) lorsquon calcule, par exemple, x2 , on utilise x1 et non x1 qui est dj` connu. ea Si lon modie litration de Jacobi de sorte ` tenir compte des rsultats les plus e a e rcents on obtient litration de Gauss-Seidel : e e
for j = 1 : n do
1 En 1823 Gauss crivait : Fast jeden Abend mache ich eine neue Auage des Tableau, wo e immer leicht nachzuhelfen ist. Bei der Einfrmigkeit des Messungsgeschftes gibt dies immer eine o a angenehme Unterhaltung ; man sieht daran auch immer gleich, ob etwas Zweifelhaftes eingeschlichen ist, was noch w nschenswert bleibt usw. Ich empfehle Ihnen diesen Modus zur Nachahmung. u Schwerlich werden Sie je wieder direct elimieren, wenigstens nicht, wenn Sie mehr als zwei Unbekannte haben. Das indirecte Verfahren lt sich halb im Schlafe ausf hren oder man kann whrend a u a desselben an andere Dinge denken. 2 Ceci a t formalis avec la relation (8.2). ee e

(k+1)

(k)

Algorithme 14 Algorithme de Jacobi

1: Donner un point de dpart x(0) Rn e 2: for k = 0, 1, 2, . . . until convergence do 3: for i = 1 : n do (k+1) (k) 4: xi = bi j=i aij xj /aii 5: end for 6: end for
(k+1) (k+1) (k)

xi = bi end for

j<i

aij xj

j>i

aij xj

/aii

Remarquons quil est dans ce cas possible de remplacer au fur et ` mesure llment a ee (k) (k+1) xi par llment xi ee dans le vecteur x.
x:
1. . . . . . . . .
. . . . . . . . . . . . . . . . . . .

i1 i i+1 . . .
......... .. . . .. . . . .. . .. . . .. . . . .. . . . . .. . . .. . . . . . . . . .

. . . . . . . . .

x(k+1)

(k+1)

x(k)

Les lments 1 ` i 1 du vecteur x contiennent dj` les valeurs de litration k + 1 ee a ea e et les lments i + 1 ` n encore les valeurs de litration prcdente. On obtient alors ee a e e e lalgorithme 15. Algorithme 15 Algorithme de Gauss-Seidel
1: Donner un point de dpart x(0) Rn e 2: for k = 0, 1, 2, . . . until convergence do 3: for i = 1 : n do 4: xi = bi j=i aij xj /aii 5: end for 6: end for

9.2

Interprtation gomtrique e e e

La prsentation gomtrique du fonctionnement des algorithmes de Jacobi et Gausse e e Seidel permet de montrer pourquoi la normalisation et lordre des quations (pour e Gauss-Seidel seulement) inuencent la convergence de ces mthodes. e Pour cette illustration nous rsolvons le syst`me donn par les deux quations suie e e e vantes : eq1 : 2x1 1x2 = 0 eq2 : 5x1 + 9x2 = 0 Ce syst`me dquations admet deux normalisations, soit lquation 1 explique x1 et e e e lquation 2 explique x2 , soit lquation 1 explique x2 et lquation 2 explique x1 . e e e Pour la rsolution on peut considrer les quations soit dans lordre quation 1 puis e e e e

quation 2, soit dans lordre inverse. Ci-apr`s les quatre syst`mes quivalents qui en e e e e rsultent : e N1 O1 : N1 O2 : x1 x2 x2 x1 = 0
5 9 1 2

0
5 9

x1 x2 x2 x1

N2 O1 : N2 O2 :

x2 x1 x1 x2

0 2 9 0 5 0 9 5 2 0

x2 x1 x1 x2

0
1 2

Considrons les itrations de Jacobi appliques aux syst`mes dquations N1 O1 et e e e e e N2 O1 en choisissant x(0) = [6 5] comme point de dpart : e N1 O1 : x1 (1) x2
(2) x1 (2) x2 (1)

= = = = = =

1 (0) x 2 2 5 (0) x 9 1 1 (1) x 2 2 5 (1) x 9 1 1 (2) x 2 2 5 (2) x 9 1

= = = = = =

5 2 10 3 5 3 25 18 25 36 25 27

N2 O1 : x2 = 2x1 = 12 (1) 9 (0) x1 = 5 x2 = 9 x2 = 2x1 = 18 (2) (1) x1 = 9 x2 = 108 5 5

(2) (1)

(1)

(0)

x1 (3) x2

(3)

A chaque itration lalgorithme de Jacobi op`re n rsolutions unidimensionnelles e e e indpendantes. Chaque quation est rsolue par rapport ` la variable spcie par e e e a e e la normalisation, les autres variables restent xes ` la valeur donne au dbut de e a e e litration. e
. . .. . .. . ... .. . .. . . . .. . . . 2 1 .. . . . .. . .. .. . . .. .. ... .. ... . .... .. ... ...... ...... .. . . .. ... . ......................................................... ....... . .. .......................................................... . . . . . . . . . . . . . . ........... . . .. .. .. .. .. . ... .... . . . 6 . . . .. . .... . .... .. .. .. . . x(0) = . .... .. .. . . . . .... . .. . .. .. .... . . .... . . . . . . .. . .. . . . 5 .............. . . .. .. . . . .. . .. . . . . . . . . . . .. . .. . .. .... . .... . . . . . .. . . . .... . .... . .. . . . .... 1 . .... .... . 1 ... . . .. ... ... . . ... .. . . . . .. . .... . .. . . .. .. .... . . .. . . . . ...... . .. .... .. .... ..... . . .. . ..... . ........................................................... ...... .. .. . ......................................................... .. . . . .. ... ... . . . .... . .... .. . .. 5 . . . .. .... .... . (1) . .... .... .. . = 2 . .... ... . . x .. 2 ... ... 10 . . .. .... . .... .. . . 9x2 5x1 =0 3 . .... .... . .. . .... . .... .. . ... ... . . . .. .... .... . . .... .... . .. . . .... .... . .. .. ... ... .... . ... .. .... .. .. .... . .. .. . (2) ....... ... . .. . .. ... x ....... .. .. . .. (3) .. .. ... . x ..... .. .... .... .. . .. . ... .... ... . ... .. .. ..... ... ... . .... .. .... . . .... .... .. .... . .... .. ... ... . .... .. .... . . ..... .. ....... . .. . .. . ..... ... .... .... .... ... . ... . .... .. .... .. . .. . . .... .... . .... .... .. .... ... . .. . .

eq1 :

2x1 x2 =0

N O

5 4 3 2 1

N O

eq2 :

Fig. 9.1 Algorithme de Jacobi. La gure 9.1 montre comment lalgorithme converge pour la normalisation 1 et diverge pour la normalisation 2. Etant donn que les composantes de la solution e

sont calcules indpendamment il appara clairement que lordre des quations est e e t e sans incidence sur le fonctionnement de lalgorithme de Jacobi. Lors de litration de Gauss-Seidel la rsolution unidimensionnelle ne se fait pas de e e mani`re indpendante mais on tient successivement compte de la mise ` jour des e e a composantes du vecteur de solution. Lalgorithme de Gauss-Seidel est donc sensible ` la normalisation et ` lordre des quations comme le montre son application aux a a e quatre syst`mes spcis plus haut : e e e N1 O1 : x1 (1) x2 x1 (2) x2 N2 O1 : x2 (1) x1
(1) (2) (1)

= = = =

1 (0) x 2 2 5 (1) x 9 1 1 (1) x 2 2 5 (2) x 9 1

= 2.5 = 1.38 = 0.69 = 0.38

5 N1 O2 : x2 = 9 x1 = 3.33 (1) (1) x1 = 1 x2 = 1.66 2

(1)

(0)

x2 = 5 x1 = 0.92 9 (2) (2) x1 = 1 x2 = 0.40 2 N2 O2 : x1 = 9 x2 = 9 5 (1) (1) x2 = 2x1 = 18

(1) (0)

(2)

(1)

= 2x1 = 12 9 (1) = 5 x2 = 21.6

(0)

Dans la gure 9.2 on peut suivre le cheminement de la solution en fonction de la normalisation et de lordre des quations. e
.. .. . . . . .. .. . .. . ... .. ... . 2 1 .. . . .. . . 1 1 2 2 .. .... .... . . .. . .. .. .. .. .. .. .. .. .. .. .. .. .. .. . . . . . . . . . . . . . . . ...... ...... .... ... .. ......... .. .. .. .. .. .. .. .. .. .. .. .. .. ... . ......... .......... .. .. . .. . . .. .... .. . .... 6 ... . (0) . .... .. . ... ... = . . . x .... .. . .... . . . .. .. . . 5 ............... . .. . .. . . . . . . .... .... .. . . .... . .... . .. 1 2 ... .... .... . . .. ... .. ... . . .. . .. .. .... .... .. .. . ..... . . ...... . . . ... . .. . .. ... . .. . .... . ... . .. ... ... . . . .... . .. .... . . . . .... .. .... . . . . . .... .... .. . . . .... .... 2 . .. .... .... . .. . ... ... . 9x2 5x1 =0 . .... .... .. . . .... . .... . .. .... .... . . .. . ... ... . . .... .... .. . .. ... . .... .. ... .... .. . .. .... . .... .. . . ...... . . . .... .. .. . .. . ...... .. .. .. .. .. .. ...... .. . .. ........ .. .. .. .. ........ . .... 5 ... ... .. .... .... .. . ... ... . . x(1) = 2 .. . .... .... . .. .... 25 .... .. ... . .. . .... . . .. .. .. .. ...... 18 . ...... . ... . .. ... . .... .. ... . . ..... x(2) .. ....... .. .. ....... . .... .... . ..... ..... .... . .... . ... .. ... .. .... . .... ... ... . .. .... .... . .. .

eq1 :

2x1 x2 =0

N O

5 4 3 2 1

N O

eq2 :

d d

Fig. 9.2 Algorithme de Gauss-Seidel.

9.3

Convergence de Jacobi et Gauss-Seidel

An de pouvoir tudier la convergence des mthodes itratives il convient de formalie e e ser ces algorithmes de mani`re dirente. Considrons la matrice A et sa dcomposition e e e e A=L+D+U :

a11 0 a12 a1n 0 0 . . . . . . .. .. . . a22 a21 . . . . . . + +. . . .. . . .. .. . . . .. . . . an1,n . . . an1 an,n1 0 0 0 ann

Litration gnrique pour la mthode de Jacobi scrit alors e e e e e Dx(k+1) = (L + U)x(k) + b et pour la mthode de Gauss-Seidel on a e (D + L)x(k+1) = Ux(k) + b. Ainsi les mthodes itratives peuvent se gnraliser sous la forme e e e e Mx(k+1) = Nx(k) + b (9.1)

avec A = M N une dcomposition de la matrice A. Pour la mthode de Jacobi e e on a M = D et N = (L + U) et pour la mthode de Gauss-Seidel M = D + L et e N = U. Pour quil ait un intrt pratique, le syst`me linaire (9.1), cest-`-dire ee e e a Mx(k+1) = c doit tre facile ` rsoudre. Ceci est le cas dans nos exemples, tant donn que M e a e e e est diagonale pour la mthode de Jacobi et triangulaire pour la mthode de Gausse e Seidel. Montrons que la convergence de (9.1) vers x = A1 b dpend des valeurs propres de e 1 M N. Dnition 9.1 Soit une matrice B Rnn , alors son rayon spectral est dni e e comme (B) = max{ || tel que est une valeur propre de B}. Thor`me 9.1 Soit b Rn et A = M N Rnn non-singuli`re. Si M est e e e 1 (k) non-singuli`re et si (M N) < 1, alors les itrations x dnies par Mx(k+1) = e e e Nx(k) + b convergent vers x = A1 b, quelque soit la valeur initiale x(o) . Dmonstration 7 Soit e(k) = x(k) x lerreur a litration k. Comme Ax = b et e ` e

M x = N x + b et en remplacent x par x(k) e(k) on dduit que lerreur pour x(k+1) e

est donn par e(k+1) = M 1 N e(k) = (M 1 N )k e(o) . Lon sait que (M 1 N )k 0 si et e seulement si (M 1 N ) < 1. Ax = b M (x(k+1) x) = N e(k)
e(k+1) (k+1)

M x = N x + b = N (x(k) e(k) ) + b = N x(k) + b N e(k)

M x(k+1)

= M 1 N e(k) = (M 1 N )k e(o)

Ainsi la convergence de (9.1) dpendra de (M 1 N). e

Conditions susantes pour la convergence Une condition susante mais pas ncessaire qui garantit la convergence pour litration e e de Jacobi est la dominance diagonale. Une matrice A Rnn est dite ` diagonale a dominante stricte si
n

|aii | > Dans ce cas on vrie e

j=1

|aij | i = 1, . . . , n.

j=i

(M 1 N)

D 1 (L + U)

= max

1in
j=1

aij < 1. aii

j=i

Pour la mthode de Gauss-Seidel on peut monter que litration converge toujours e e si A est symtrique et dnie positive (condition susante mais pas ncessaire). e e e

9.4

Mthodes de relaxation (SOR) e

Lorsque le rayon spectral de M 1 N est proche de lunit, la mthode de Gauss-Seidel e e converge tr`s lentement. On peut y remdier en modiant literation de Gauss-Seidel e e de la faon suivante : c
xi
(k+1)

= xGS,i + (1 )xi
(k+1)

(k+1)

(k)

avec xGS,i la valeur de xi calcul avec un pas de Gauss-Seidel et le param`tre e e de relaxation. Ceci dnit la mthode de la sur-relaxation successive (SOR). On voit e e quil sagit dune combinaison linaire dun pas de Gauss-Seidel et le pas prcdent e e e de SOR. On a lalgorithme :

(k+1)

Algorithme 16 Sur-relaxation successive (SOR)

1: Donner un point de dpart x(0) Rn e 2: for k = 0, 1, 2, . . . until convergence do 3: for i = 1 : n do (k+1) (k+1) (k) 4: xi = xGS,i + (1 )xi 5: end for 6: end for

Choix du param`tre de relaxation e An dtudier le choix des valeurs de formalisons litration matriciellement. On e e peut crire e M x(k+1) = N x(k) + b avec M = D + L et N = (1 )D U. La valeur du param`tre de relaxation e qui minimise le rayon spectral nest connue que pour quelques probl`mes particuliers. e En gnral on dnit par ttonnement ou par balayage. e e e a Montrons encore que le param`tre de relaxation doit vrier 0 < < 2. En eet e e on a
1 det(M N ) = det(D + L)1 det (1 )D U

= (1 )n < 1

= det (1 )D / det(D)

1 1 et comme dautre part det(M N ) = n i on voit que la condition (M N ) < i=1 1 implique 0 < < 2.

Rappelons que si > 1 videmment aucune des mthodes ne converge. e e

9.5

Structure des algorithmes itratifs e

Les mthodes itratives ncessitent un nombre inni ditrations pour converger vers e e e e la solution. Ainsi, pour rendre ces mthodes oprationnelles, il convient dintroduire e e un crit`re darrt pour le procd itratif. On arrtera les itrations lorsque les e e e e e e e changements dans le vecteur de solution, cest-`-dire lerreur devient susamment a petite. Comme il a t suggr ` la section 1.3 on choisira une combinaison entre erreur ee ee a absolue et erreur relative ce qui conduira ` stopper les itrations lorsque la condition a e suivante est satisfaite |xi
(k+1) (k)

|xi | + 1

xi |

(k)

i = 1, 2, . . . , n

avec une tolrance donne. e e

Structure des algorithmes itratifs e La mise en oeuvre dalgorithmes itratifs peut alors se schmatiser avec les instruce e tions suivantes :
Initialiser x(1) , x(0) , et le nombre maximum ditrations e while converged(x(0), x(1) , ) do x(0) = x(1) (Conserver litration prcdente dans x(0) ) e e e (1) avec Jacobi, Gauss-Seidel ou SOR Evaluer x Test sur le nombre ditrations e end while

La fonction Matlab converged se prsente comme : e

function res = converged(y0,y1,tol) res = all(abs(y1-y0)./(abs(y0)+1) < tol);

Ltape qui value x ncessite n 2 p oprations lmentaires pour Jacobi et Gausse e e e ee Seidel et n (2 p + 3) oprations lmentaires pour SOR, ou p est le nombre moyen e ee dlments non nuls dans une ligne de la matrice A dordre n. D`s lors on peut ee e envisager de prfrer la mthode de Jacobi ou Gauss-Seidel ` une mthode directe ee e a e (qui nexploite pas la structure creuse de A) si n2 3p

2knp <

2 3

ou k est le nombre ditrations. e

9.6

Mthodes itratives par bloc e e

Certains probl`mes peuvent tre naturellement dcomposs en sous-probl`mes plus e e e e e au moins fortement lis entre eux. Ceci est notamment le cas pour des mod`les e e dcrivant linterdpendance de lconomie de plusieurs pays (multi-country moe e e dels). Les conomies nationales tant lies entre elles par les relations de commerce e e e extrieur seulement. Une autre situation constituent des mod`les multi-sectoriels e e dsagrgs ou lon observe que peu de liens entre les secteurs compar aux liens e e e e existant ` lintrieur dun secteur. a e

eq. pays 1

eq. pays 2

. . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. ............................................................................. ........................................................... ............... . . . . . . . . . . .. . . . . . . .. . . . .. . . . .. . . . . .. . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . .. . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . . . . .. . . .. . . .. . . .. . . .. . . . . . .

Une mthode par bloc est alors une technique ou lon it`re sur les blocs. La technique e e utilise pour rsoudre le bloc peut tre quelconque et nest pas importante pour la e e e dnition de la mthode. e e Les mthodes itratives par bloc peuvent tre appliques indiremment aux syst`mes e e e e e e linaires et aux syst`mes dquations non-linaires. Pour simplier la prsentation e e e e e on consid`re dans la suite la solution dun syst`me linaire Ay = b. Dans le cas none e e linaire la matrice A constitue la matrice Jacobienne et du fait de la non-linarit la e e e valeur des lments de la matrice A et du vecteur b changent ditration en itration. ee e e Supposons la partition suivante de notre matrice A A11 A12 A1N A21 A22 A2N A= . . . . . . . . . AN 1 AN 2 AN N

avec Aii , i = 1, . . . , N des matrices carres. En partitionnant le syst`me Ay = b de e e faon correspondante on obtient c y1 b1 A11 A1N . . . = . . . . . . . . . . . . AN 1 AN N
N

que lon peut aussi crire e

Aij yj = bi
j=1

i = 1, . . . , N.

Si les matrices Aii , i = 1, . . . , N sont non-singuli`res alors on peut mettre en oeuvre e lalgorithme suivant : Il existe une plthore de variantes pour le choix de la partition et la faon dont le e c bloc est rsolu. Une variante consiste ` choisir N = 2 avec A11 triangulaire infrieure, e a e donc facile ` rsoudre, et A22 quelconque mais de taille nettement infrieure que lon a e e rsoudra avec une mthode directe dans le cas dun syst`me linaire et la mthode e e e e e de Newton pour un syst`me non-linaire. e e

Algorithme 17 Mthode itrative par bloc e e

1: Donner un point de dpart y (0) Rn e 2: for k = 0, 1, 2, . . . until convergence do 3: for i = 1 : N do
N

solve Aii yi

(k+1)

= bi

Aij yj
j=1 j=i

(k)

5: end for 6: end for

A11

A22

Une autre alternative consiste ` rsoudre les blocs de mani`re approximative (ina e e complete inner loops). La convergence des mthodes itratives par bloc dpend comme pour les mthodes e e e e 1 itratives stationnaires des valeurs propres de M N. e

Chapitre 10 Mthodes itratives non e e stationnaires

Chapitre 11 Equations non-linaires e

Contrairement ` ce qui a t le cas pour les syst`mes dquations linaires, la a ee e e e rsolution des syst`mes non-linaires sav`re beaucoup plus dlicate. En gnral il e e e e e e e nest pas possible de garantir la convergence vers la solution correcte et la complexit e des calculs cro tr`s vite en fonction de la dimension du syst`me. Ci apr`s, ` titre t e e e a dexemple un syst`me de deux quations : e e
3 1 y1 y2 2 = 0 2 4 y1 y2 + c = 0

(11.1) (11.2)

Dans la Figure 11.1 on montre quen fonction de la valeur que prend le param`tre c e ce syst`me admet entre zro et quatre solutions. e e
1 c=0 0.5 y
1

1 c~0.7 0.5 y
1

0.5 1 1 y 1 c=1 0.5 y

1 2

0.5 0 1 1 1 y 1 c=1.5 0.5 y

1 2

0.5 1 1 y
2

0.5 0 1 1 1 y
2

Fig. 11.1 Solutions du syst`mes dquations (11.111.2) en fonction de c. e e

On consid`re un syst`me dquations e e e h(y, z) = 0 F (y) = 0

o` h : Rn Rn Rn est drivable dans un voisinage de la solution y Rn et u e m z R pour les variables exog`nes. Si au moins une des quations est non-linaire, e e e le syst`me sera dit non-linaire. e e Dans la pratique, et surtout dans le domaine de la modlisation macroconomtrique, e e e la matrice Jacobienne h/y peut souvent tre permute de sorte ` faire appara e e a tre une structure rcursive par blocs comme indique dans la gure ci-apr`s1 : e e e
.. .. .. .. .. . .... . . .... . . . ... . .. . . ... . . . . ... . ... . . ... . ... . ... . ... . . . . ............................................... . . . ... . ... . . ... . ... . ... . ... . . . .............................................. ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. ..................................... . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. .................................... . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. ..................................... . . . ... . ... . . ... . ... . ... . ... . . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . . ............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. . . . ... . ... . . ... . ... . ... . ... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . ............................................. .................................... . . . .............................................. . . . ... . ... . . ... . ... . ... . ... . . . ...................................... ... . . . . . . . . . .. .. . . . . . . . . . . . . . ... . ... . ... . . . . . . . . . . . . . . ........... ... . ... . . . . . . . . . . . . . . ........... ........... ... . ... ... . ... . . . . . . . . . . . . . ........... ........... ... . ... ... . ... . . . . . . . . . . . . . . ........... ........... ... . ... ... . ... . . . . . . . . . . . . . ........... ......... ... . ... . . . . . . . . . . . . . . ............. . . . .. . . . . . . . . . . . . . . . . ... .. . . . . . . . . . . . . . . . . . ... .. . . . . . . . . . . . . . . . . . ... .. . . . . . . . . . . . . . . . . . . ... .. .

Dans ce cas la solution du syst`me dquations sobtient en rsolvant une succession e e e de syst`mes rcursifs et de syst`mes interdpendants. Dans la suite on supposera e e e e toujours que le syst`me ` rsoudre est interdpendant, cest-`-dire que sa matrice e a e e a Jacobienne ne peut se mettre sous forme triangulaire par blocs. Les trois approches suivantes constituent les techniques les plus souvent utiliss pour e la rsolution de syst`mes dquations non-linaires : e e e e Mthode du point xe e Mthode de Newton e Rsolution par minimisation e

11.1

Equation ` une variable a

Toute valeur de x qui satisfait f (x) = 0 est appele zro de la fonction. Souvent e e dans la pratique, soit il nexiste pas de solution analytique ` ce probl`me, soit son a e obtention est tr`s laborieuse. Remarquons que toute galit peut tre mise sous cette e e e e forme. Le probl`me de la recherche des zros dune fonction est frquemment rencontr. e e e e Mentionnons ici, par exemple, le calcul de la volatilit implicite des options (donner e
1

La recherche de cette dcomposition a t illustr dans la section 8.3.1. e ee e

un exemple en conomie). En ralit, il sagit de rsoudre une quation o` , apr`s e e e e e u e avoir ramen tous les termes ` gauche du signe de lgalit, elle se prsente comme e a e e e

f (x) = 0 .

Ici on va traiter uniquement le cas o` f est univarie, f : R R. Le cas multivari, u e e f : Rm Rn , sera trait dans le cadre de la rsolution de syst`mes (linaires et e e e e non-linaires). e Il existe une grande varit de situations, plus ou moins dlicates ` traiter. Hormis ee e a dans les cas linaires, on proc`de toujours par itrations (voir mthodes itratives). e e e e e Le cas reprsent ` gauche dans la gure 11.2 est un des plus simples et favorables. e ea La fonction varie doucement, sannule et change de signe.

1 0.5

1 0.8 0.6

0 0.5

0.4 0.2 0

1 1

0.2 1

Fig. 11.2 Gauche : f (x) = cos(x + 3/8). Droite : f (x) = (x 1/2)2.

Le cas ` droite de la gure 11.2 reste relativement ais, mais est toutefois moins a e favorable que le prcdent car la fonction varie doucement, sannule mais ne change e e pas de signe. Les cas de la gure 11.3 sont similaires au cas ` la droite de la gure 11.2, sauf que a les erreurs numriques peuvent conduire ` deux situations direntes. Notamment e a e la fonction sannule presque mais reste numriquement strictement positive. Dans e ce cas, doit-on considrer que x = 1/2 est la solution dune fonction mal value, ou e e e nest pas la solution car la fonction ne sannule eectivement pas ?

1 0.8 0.6 0.4 0.2 0 0.2 1 0 1

Fig. 11.3 Gauche : f (x) = (x 1/2)2 + . Droite : f (x) = (x 1/2)2 . Dans lautre cas la fonction sannule deux fois car elle prend des valeurs tr`s faiblee ment ngatives : dans ce cas, peut-on considrer que les deux solutions distinctes e e sont acceptables ou sont seulement les consquences dune mauvaise valuation de e e la fonction ? Dans la gure 11.4 ` nouveau deux situations dlicates. Dans le premier cas, la a e fonction prsente plusieurs zros ; dans ce cas, comment peut-on obtenir toutes les e e solutions ? Un algorithme donn converge vers une solution particuli`re, mais lae e quelle ? Comment passer dune solution ` la suivante ? a
1 0.5 0 0.5 1 0.2 100 50 0 50 100 0

0.4

0.6

0.8

0.5

Fig. 11.4 Gauche : f (x) = cos(1/x2 ). Droite : f (x) = 1/(x 1/2). Le second cas correspond ` une fonction prsentant une singularit et changeant a e e de signe de part et dautre de la singularit ; dans un tel cas, lalgorithme risque de e croire que le point de singularit est la solution cherche ; dautre part, lvaluation e e e de la fonction en ce point risque de poser probl`me (division par zro, grande erreur e e de prcision, . . . ) e Les fonctions de la gure 11.5 prsentent une discontinuit. La fonction H tant e e e dnie comme e 1 x < 0 H(x) = 1 x>0

Une telle discontinuit ne peut pas tre prise en compte numriquement, elle ncessiterait e e e e une prcision innie ; dans lexemple illustr ici, la valeur x = 0 serait considre e e ee comme solution de lquation f (x) = 0. e
f(x) = 2*H(x)1 1.5 1 0.5 0 0.5 1 1.5 1 0 1 1.5 1 0.5 0 0.5 1 1.5 1 0 1

Fig. 11.5 Gauche : f (x) = 2H(x) 1. Droite : f (x) = 1/(x 1/2). Il est en gnral assez dlicat de donner une rponse universelle ` ces questions ; e e e e a seule une tude particuli`re adapte ` la situation permet de rsoudre ce probl`me. e e e a e e Ces quelques exemples montrent quil est dicile de dnir un algorithme numrique e e universel, et quil est indispensable de conna a priori certaines informations relatre tives ` la fonction et ` son/ses zro(s). Cest alors en fonction de cette connaissance a a e a priori que lon cherchera la mthode numrique la plus adapte. e e e

Crit`res pour le choix de la mthode e e

Avant de choisir une mthode numrique pour la recherche des zros dune fonce e e tion, o` dexplorer lexistence dune solution analytique il convient de se poser les u questions suivantes : Doit-on procder ` des recherches rptes ? Si oui, peut-tre convient-il dexplorer e a e ee e lexistence dune solution analytique. Quelle est la prcision dsire ? Si quelques digits susent, alors une mthode e e e e quelconque peut-tre utilise. e e La mthode doit-elle tre rapide et robuste (insensible aux choix des points de e e dpart) ? e Sagit-il dun polynme ? Dans ce cas il existe des procdures spciales. o e e La fonction comporte-elle des singularits ? e

Principe de la recherche des zros e

Considrons la gure ?? qui illustre une fonction avec deux zros, x1 et x2 dans le e e domaine des x considr. Toute recherche de zro dmarre avec une solution initiale ee e e

qui sera amliore dans des itrations successives, et dans des conditions normales, e e e sapproche de la solution. Quelques remarques : Le choix du point de dpart determine vers quelle solution on se dirigera. e Il faut dnir un crit`re darrt pour les itrations. e e e e Il convient de substituer la solution dans la fonction et de vrier que lon a e f (xsol ) 0. La gure ?? illustre quune fonction qui, pour un intervalle donn, vrie des e e signes opposs peut soit contenir un zro, soit une singularit. e e e

11.1.1

Recoupement par intervalles (bracketing)

Il sagit de construire des intervalles susceptibles de contenir un zro. Ultrieurement e e on ranera la recherche du zro ` lintrieur de lintervalle. e a e On subdivise un domaine donn en intervalles rguliers et on examine si la fonction e e traverse laxe des x en analysant le signe de la fonctions value aux bords de e e lintervalle (voir gure ??).

Algorithme 18 Etant donn f (x), xmin , xmax et n on identie les intervalles suscepe tibles de contenir des zros. e
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: = (xmax xmin )/n a = xmin i=0 while i < n do i=i+1 b=a+ if sign f (a) = sign f (b) then [a, b] peut contenir un zro, imprimer a et b e end if a=b end while

La fonction bracketing ralise lalgorithme 18 et la gure 11.6 illustre la recherche e de ces intervalles pour la fonction f (x) = cos(1/x2 ) dans le domaine x [0.3, 0.9] et n = 25.

1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 11.6 f (x) = cos(1/x2 ), n = 25.

Code

/Teaching/MNE/Ex/bracketing.m:

function ab = bracketin g (f , xmin , xmax , n ) k = 0; x = linspace ( xmin , xmax , n ); a = xmin ; fa = feval (f , a ); for i = 2: n b = x ( i ); fb = feval (f , b ); if sign ( fa )~= sign ( fb ) k = k + 1; ab ( k ,:) = [ a b ]; end a = b; fa = fb ; end

Pour tester si les fonctions sont de signe oppos en a et b on pourrait crire e e f (a) f (b) < 0 . Dun point de vue numrique il est cependant prfrable davoir recours ` la fonction e ee a sign. Lexemple qui suit montre comment pour des valeurs tr`s petites de f le e rsultat peut correspondre ` un underow. e a
>> fa = 1e -140; >> fa * fb ans = 0 fb = 1e -180;

11.1.2

Itrations de point xe e

Soit une fonction f (x) = 0, on isole un terme contenant x de la sorte ` pouvoir a crire e xnew = g(xold ) . On appelle la fonction g(x) la fonction ditration et une itration du point xe est e e dnie par lalgorithme qui suit. e Algorithme 19
1: Initialiser x(0) 2: for k = 1, 2, . . . jusqu` convergence do a 3: x(k) = g(x(k1) ) 4: end for

A propos de lalgorithme du point xe il convient de faire les remarques suivantes : La notion de convergence doit tre approfondie. e Les valeurs successives de x ne doivent pas tre conserves dans un vecteur. Il e e (k) (k+1) sut de conserver deux valeurs successives x et x dans des variables x0 et x1 et de remplacer x0 par x1.
while ~ c o n v e r g e d x1 = g ( x0 ) x0 = x1 end

Le fait que la solution xsol vrie xsol = g(xsol ) fait quon appelle xsol un point xe. e Le choix de la fonction g(x) est dterminant pour la convergence de la mthode. e e

Code /Teaching/MNE/Ex/FPI.m:
function x1 = FPI (f , x1 , tol ) % FPI (f , x0 ) I t e r a t i o n s du point fixe % if nargin == 2 , tol = 1e -8; end it = 0; itmax = 100; x0 = realmax ; while ~ converged1 ( x0 , x1 , tol ) x0 = x1 ; x1 = feval (f , x0 ); it = it + 1; if it > itmax error ( Maxit in FPI ); end end

Exemple 1 f (x) = x x4/5 2 = 0 g1 (x) = x4/5 + 2

f = inline(x.^(4/5)+2);

g2 (x) = (x 2)5/4

FPI(f,8) ans = 6.4338 f = inline((x-2).^(5/4)); FPI(f,8) ??? Error using ==> fpi Maxit in FPI

Exemple 2 f (x) = x 2x4/5 + 2 = 0 g1 (x) = 2x4/5 2

bracketing(f,0,50,30) ans = 3.4483 5.1724 18.9655 20.6897 f = inline(2*x.^(4/5)-2); FPI(f,1) ??? Error using ==> fpi Maxit in FPI FPI(f,9) ??? Error using ==> fpi Maxit in FPI f = inline(((x+2)/2).^(5/4)); FPI(f,1) ans = 3.7161 FPI(f,20) ??? Error using ==> fpi Maxit in FPI

g2 (x) =

x+2 2

5/4

Convergence de la mthode de point xe e Les itrations de point xe convergent pour x [a, b] si lintervalle contient un zro e e et si |g (x)| < 1 pour x [a, b] .

Si 1 < g (x) < 0 les itrations oscillent autour de la solution et si 0 < g (x) < 1 e les itrations convergent de faon monotone. Voir les gures qui suivent. e c

11.1.3

Bisection

On consid`re un intervalle qui contient le zro, on le divise en deux et on cherche e e le demi-intervalle qui contient le zro. La procdure est ritre sur le nouveau e e e ee

intervalle contenant le zro. A chaque itration lintervalle est divis par deux. Les e e e centres a+b c= 2 des intervalles successifs convergent vers le zro. En valuant le centre comme e e c=a+ ba 2

on obtient des rsultats qui sont numriquement plus stables. e e

f (a)

f (c)

f (b)

Lalgorithme 20 dcrit la mthode. e e Algorithme 20 Bisection ( est la tolrance derreur). e

1: Vrier que f (a) f (b) < 0 e 2: while |a b| > do 3: c = a + (b a)/2 4: if f (c) f (a) < 0 then 5: b=c (z est a gauche de c) ` 6: else 7: a=c (z est a droite de c) ` 8: end if 9: end while

Le code qui suit donne la ralisation avec Matlab de lalgorithme de la bisection. e

function c = bisection (f ,a ,b , tol ) % R e c h e r c h e zero d une f o n c t i o n avec m e t h o d e de la b i s e c t i o n % On c h e r c h e zero dans i n t e r v a l l e [a , b ] avec t o l e r a n c e tol % La f o n c t i o n doit etre de signe oppose en a et b % if nargin == 3 , tol = 1e -8; end fa = feval (f , a ); fb = feval (f , b ); if sign ( fa ) == sign ( fb ) error ( fonction pas de signe oppose en a et b ); end fait = 0; while abs ( b - a ) > 2* tol & ~ fait c = a + ( b - a ) / 2; % C h e r c h e r centre i n t e r v a l l e fc = feval (f , c ); % E v a l u e r f au centre

if sign ( fa ) ~= sign ( fc ) % zero a gauche de c b = c; fb = fc ; elseif sign ( fc ) ~= sign ( fb ) % zero a droite de c a = c; fa = fc ; else % on tombe e x a c t e m e n t sur zero fait = 1; end end

Exemple 3
f = inline(x-2*x.^(4/5)+2); bisection(f,0,5,1e-8) ans = 3.7161 bisection(f,5,20,1e-8) ans = 19.7603

Convergence Soit 0 = b a lintervalle initial, alors lintervalle ` litration n se dduit comme a e e 1 0 2 1 1 2 = 1 = 0 2 4 . . . 1 n n = 0 2 1 = La rduction relative apr`s n itrations est n / 0 = 2n dou on tire le nombre e e e ditrations ncessaires pour atteindre une rduction donne de lintervalle, soit e e e e n = log2 Voir crit`re de convergence. e n . 0

11.1.4

Mthode de Newton e

Parfois appele Newton-Raphson. e

11.2
11.2.1

Syst`mes dquations e e
Mthodes itratives e e

Gauss-Seidel et Jacobi Les mthodes de Gauss-Seidel et Jacobi vus dans la section 9.1 pour la solution de e syst`mes linaires peuvent tre tendus pour la rsolution de syst`mes non-linaires. e e e e e e e En gnral on normalise les quations, cest-`-dire on les rcrit sous la forme e e e a ee yi = gi (y1, . . . , yi1 , yi+1 , . . . , yn , z) i = 1, . . . , n.

Rappelons que les quations gi peuvent maintenant tre non-linaires. e e e Pour la mthode de Jacobi litration gnrique scrit alors : e e e e e yi
(k+1) (k) = gi (y1 , . . . , yi1, yi+1 , . . . , yn , z) (k) (k) (k)

i = 1, . . . , n.

Litration gnrique de la mthode de Gauss-Seidel utilise les i 1 composantes e e e e mises ` jour de y (k+1) aussitt quelles sont disponibles : a o yi
(k+1)

= gi (y1

(k+1)

(k) , . . . , yi1 , yi+1, . . . , yn , z)

(k+1)

(k)

i = 1, . . . , n.

On utilise le mme crit`re darrt que pour le cas linaire, cest-`-dire on stoppera e e e e a les iterations lorsque la condition suivante |yi
(k+1) (k)

est satisfaite ou est une tolrance donne. e e La convergence ne peut tre tablie au pralable car la matrice M 1 N du thor`me 9.1 e e e e e change ` chaque itration. La structure de lalgorithme itratif est identique ` celle a e e a prsente pour le cas linaire : e e e
Initialiser y (1) , y (0) , et le nombre maximum ditrations e (0) , y (1) , ) do while converged(y y (0) = y (1) (Conserver litration prcdente dans y (0) ) e e e Evaluer y (1) avec Jacobi, Gauss-Seidel ou SOR Test sur le nombre ditrations e end while

|yi | + 1

yi |

(k)

i = 1, 2, . . . , n

11.2.2

Formalisation de la mthode du point xe e

Il sagit dune formulation gnrale des mthodes itratives ou le syst`me dquations e e e e e e est formalis comme e y = g(y)

et litration du point xe scrit e e y (k+1) = g(y (k)) k = 0, 1, 2, . . .

tant donn un vecteur de dpart y (0) . La condition de convergence est e e e g(y sol ) < 1 o` u g
1 sol y1 =y1

est la matrice Jacobienne value ` la solution y sol . e e a

g(y sol ) = . . .

g1 yn

. . .

sol yn =yn

gn y1 y =y sol 1 1

Exemple 11.1 Illustrons la mthode du point xe en rsolvant le syst`me dquations e e e e

non-linaires suivant : e F (y) = 0 f1 (y1 , y2 ) : y1 + y2 3 = 0 2 2 f2 (y1 , y2 ) : y1 + y2 9 = 0 0 3

La Figure 11.7 montre les deux solutions y = syst`me dquations. e e

et y =

quadmet ce

2 0 y
2

Fig. 11.7 Solutions du syst`me dquations de lexercise 11.1. e e

En choisissant comme point de dpart y (0) = 1 5 e et une tolrance de 105 , lalgoe rithme du point xe pour ce probl`me peut scrire comme : e e

y1 = [1 5]; tol = 1e-5; y0 = 1+y1; % Pour premier test dans while while ~converged(y0,y1,tol) y0 = y1; y1(1) = -y0(2) + 3; y1(2) = sqrt(-y0(1)^2 + 9); end La fonction Matlab converged a t prsente a la page 75. La Figure 11.8 montre comee e e ` ment lalgorithme chemine vers la solution.
Iter k 0 1 2 3 4 5 6 7 8
2.2 2.8 y
2

0.76 0.17 y1

2 5

Solution (k) y2 1.0000 5.0000 2.0000 2.8284 0.1716 2.2361 0.7639 2.9951 0.0049 2.9011 0.0989 3.0000 0.0000 2.9984 0.0016 3.0000 0.0000 3.0000
(k) y1

(k) y1

Erreur (k) sol y1 y2 y2 1.0000 2.0000 2.0000 0.1716 0.1716 0.7639 0.7639 0.0049 0.0049 0.0989 0.0989 0.0000 0.0000 0.0016 0.0016 0.0000 0.0000 0.0000
sol

Fig. 11.8 Fonctionnement de lalgorithme du point xe avec y (0) = [ 1 5 ].

En choisissant comme valeur initiale y = 0 2 lalgorithme oscille entre y = et y = 3 3 sans converger. Ceci est montr dans la Figure 11.9. e

0 0

2 y1

0 y2

Fig. 11.9 Fonctionnement de lalgorithme du point xe avec y (0) = [ 0 2 ].

11.2.3

Mthode de Newton e

La mthode de Newton pour la rsolution dun syst`me dquations est une gnralisation e e e e e e de lalgorithme de Newton pour la recherche du zro dune fonction unidimensione nelle. Le fait quun algorithme pour un probl`me unidimensionnel puisse tre ee e cacement gnralis au cas multidimensionnel est une situation exceptionnelle. La e e e mthode de Newton est souvent aussi appele mthode de Newton-Raphson. e e e Dans la formulation classique on crit le syst`me dquations comme : e e e f1 (y) = 0 . . F (y) = 0 . f (y) = 0 n

Cette criture est quivalente ` la notation h(y, z) = 0, ` la dirence que les e e a a e variables z sont directement dans F . On approche la solution y par la squence {y (k)}k=0,1,2,.... Etant donn y (k) Rn et e e une valuation de la matrice Jacobienne e f1 f1 yn (k) (k) yn =yn y1 y1 =y1 . . (k) . F (y ) = . . . fn fn yn (k) (k) y1
y1 =y1 yn =yn

on construit une approximation meilleure y (k+1) de y . Pour ce faire on approche F (y) dans le voisinage de y (k) par une fonction ane : F (y) F (y (k)) + F (y (k) )(y y (k) ). On peut rsoudre ce mod`le local pour obtenir une valeur de y qui satisfasse e e F (y (k)) + F (y (k))(y y (k)) = 0 cest-`-dire a y = y (k) F (y (k))
1

F (y (k)).

On construira alors un nouveau mod`le local autour de la valeur obtenue pour y. A e litration k on a e 1 y (k+1) = y (k) F (y (k)) F (y (k)) et y (k+1) est solution du syst`me linaire e e F (y (k) ) (y (k+1) y (k)) = F (y (k) ) .
J s b

Algorithme 21 Mthode de Newton. e

1: Donner y (0) 2: for k = 0, 1, 2, . . . until convergence do 3: Evaluer b = F (y (k) ) et J = F (y (k) ) 4: Vrier la condition de J e 5: solve J s = b 6: y (k+1) = y (k) + s 7: end for

Exemple 11.2 Illustrons la mthode de Newton en rsolvant le syst`me dquations de e e e e

lexemple 11.1 soit F (y) = 0 f1 (y1 , y2 ) : y1 + y2 3 = 0 2 2 f2 (y1 , y2 ) : y1 + y2 9 = 0

et dont la matrice Jacobienne scrit e F (y (k) ) = f

f1 f1 y1 y =y (k) y2 y =y (k) 1 2 1 2 f2 2 y1 y =y (k) y2 y =y (k) 1 2 1 2

1 1 = . (k) (k) 2y1 2y2

Si lon choisit y (0) = 1 5

alors 3 17 et F (y (0) ) = 1 1 2 10 .

F (y (0) ) =

En rsolvant le syst`me e e 1 2 lon obtient s(0) = 13/8 11/8 y (1) = y (0) + s(0) = .625 , 3.625

1 10 do` u

s(0) =

3 17

F (y (1) ) =

0
145 32

F (y (1) ) =

1 1 . 5 29 4 4

En rsolvant a nouveau le syst`me e ` e 1 1 5 29 4 4 on obtient s(1) = 145/272 145/272

s(1) = do` u

0 145 32

y (2) = y (1) + s(1) =

.092 . 3.092

En se servant du code Matlab qui suit (en initialisant y1) :

tol = 1e-5; y0 = 1+y1; % Pour premier test dans while while ~converged(y0,y1,tol) y0 = y1; F(1) = y0(1) + y0(2) - 3; F(2) = y0(1)^2 + y0(2)^2 - 9; b = -F; J = [1 1; 2*y0(1) 2*y0(2)]; s = J \ b; y1 = y0 + s; end on peut calculer les itrations jusqu` la convergence et la Figure 11.10 montre comme ces e a itrations convergent vers la solution y = 0 3 . Rappelons que la fonction Matlab e converged a t introduite a la page 75. ee `

y1 0.09 0.62

Iter k 0 1 2 3 4

Solution (k) (k) y1 y2 1.0000 5.0000 0.6250 3.6250 0.0919 3.0919 0.0027 3.0027 0.0000 3.0000

Erreur (k) (k) sol sol y1 y1 y2 y2 1.0000 2.0000 0.6250 0.6250 0.0919 0.0919 0.0027 0.0027 0.0000 0.0000

3.09 3.62 y2

Fig. 11.10 Fonctionnement de lalgorithme de Newton avec y (0) = [ 1 5 ].

Remarquons que la mthode de Newton converge aussi lorsquon choisit comme valeur e ce qui est illustr dans la Figure 11.11. e initiale y (0) = 0 2

Convergence Le taux de convergence r dune mthode itrative est dni comme e e e lim e(k+1) =c e(k) r

ou e(k) est lerreur ` litration k et c une constante nie. On distingue alors les cas a e particuliers suivants : r = 1 et c < 1, convergence linaire e r > 1, convergence super-linaire e r = 2, convergence quadratique

3.25 3

0.25 0 y2

Fig. 11.11 Fonctionnement de lalgorithme de Newton avec y (0) = [ 0 2 ].

Dans la pratique ceci correspond ` un gain de prcision par itration dun nombre a e e de digits constant pour une convergence linaire ; laccroissement de la prcision par e e itration va en augmentant pour une convergence super-linaire, et dans le cas dune e e convergence quadratique la prcision double ` chaque itration. e a e La mthode de Newton converge quadratiquement sous certaines conditions. On e vrie e y (k+1) y y (k) y 2 k = 0, 1, 2, . . . (11.3) o` mesure la non-linarit relative u e e Lipschitz. F (y )1 < 0 et est la constante de

La convergence est seulement garantie si y (0) se trouve dans un voisinage de y , voisinage quil faut dnir. e Pour la solution de mod`les macroconomtriques, y (0) est dni par yt1 ce qui e e e e constitue en gnral un bon voisinage. e e La convergence quadratique des itrations de lexemple 11.2 peut tre vrie dans e e e e (k) le tableau de la Figure 11.10. En eet pour chaque composante yi , on vrie e |yi
(k+1)
sol sol yi | < |yi yi |2 .

(k)

La complexit de la mthode de Newton est dtermine par la rsolution du syst`me e e e e e e linaire. Si lon compare les mthodes itratives avec la mthode de Newton pour e e e e une itration on remarque que la dirence de travail rside dans lvaluation de e e e e la matrice Jacobienne et la solution dun syst`me linaire. Ceci peut para e e tre un dsavantage de la mthode de Newton. Dans la pratique cependant on peut montrer e e que lorsque on rsoud le syst`me linaire avec une mthode directe creuse et lorsquon e e e e utilise des drives analytiques les deux mthodes sont de complexit comparable e e e e tant donn que le nombre ditrations pour Newton est nettement infrieur aux e e e e nombre ditrations pour les mthodes itratives. e e e

11.2.4

Quasi-Newton

La mthode de Newton ncessite ` chaque itration le calcul de la matrice Jacoe e a e 2 bienne, ce qui requiert lvaluation de n drives, et la rsolution dun syst`me e e e e e linaire (O(n3)). Dans la situation ou lvaluation des drives est co teuse on peut e e e e u remplacer le calcul exact de la Jacobienne par une simple mise ` jour. Ceci donne a lieu ` des nombreuses variantes de la mthode de Newton, appellees mthodes quasia e e Newton.

Mthode de Broyden e Dans ce cas la mise ` jour de la matrice Jacobienne se fait avec des matrices de rang a un. Etant donn une approximation B (k) de la matrice Jacobienne ` litration k on e a e calcule lapproximation de litration k + 1 comme e dF (k) B (k) s(k) s(k) s(k) s(k)

(k+1)

(k)

ou dF (k) = F (y (k+1)) F (y (k) ) et s(k) est la solution de B (k) s(k) = F (y (k) ). Lalgorithme de Broyden est formalis ci-apr`s. e e Algorithme 22 Algorithme de Broyden.
1: Donner y (0) et B (0) (approximation de F (y (0) )) 2: for k = 0, 1, 2, . . . until convergence do 3: solve B (k) s(k) = F (y (k) ) 4: y (k+1) = y (k) + s(k) 5: = F (y (k+1) ) F (y (k) ) 6: B (k+1) = B (k) + B (k) s(k) s(k) / (s(k) s(k) ) 7: end for

Le code Matlab qui suit ralise la mthode de Broyden. e e y0 = - y1; tol = 1e-5; B = eye(2); F1 = feval(ExSNL,y1); while ~converged(y0,y1,tol) y0 = y1; F0 = F1; s = B \ -F0; y1 = y0 + s; F1 = feval(ExSNL,y1); dF = F1 - F0; B = B + ((dF - B*s)*s)/(s*s); end

Lvaluation de la fonction F (y) a t dans ce cas transfre dans une fonction Mate ee ee lab ExSNL ce qui rend le code indpendant de la rsolution dun syst`me dquations e e e e particulier. function F = ExSNL(y) F = repmat(NaN,2,1); F(1) = y(1) + y(2) - 3; F(2) = y(1)^2 + y(2)^2 - 9; La Figure 11.12 donne les rsultats de mthode de Broyden en partant avec B (0) = I e e (0) et pour une valeur initiale y = [ 2 0 ].
Iter k 0 1 2 3 4 5 6 7 8 9 10 Solution (k) y2 2.0000 0 3.0000 5.0000 2.1667 0.8333 1.6585 1.4634 4.4582 1.4984 2.3304 0.6762 2.7384 0.2614 3.0847 0.0852 2.9921 0.0079 2.9998 0.0002 3.0000 0.0000
(k) y1 (k) y1

4.45

3 y

1.49

0 y
2

Erreur (k) sol y1 y2 y2 1.0000 0 0 5.0000 0.8333 0.8333 1.3415 1.4634 1.4582 1.4984 0.6696 0.6762 0.2616 0.2614 0.0847 0.0852 0.0079 0.0079 0.0002 0.0002 0.0000 0.0000
sol

Fig. 11.12 Rsultats de lalgorithme de Broyden avec B (0) = I et y (0) = [ 2 0 ]. e

Lorsque lalgorithme dmarre avec B (0) = I aucune information sur la matrice e Jacobienne nest donn et il nest pas surprenant que la performance soit mdiocre e e (comparable ` la mthode du point xe). Ci-apr`s on fait dmarrer lalgorithme a e e e avec B (0) = F (y (0) ) et on peut observer que dans ce cas sa performance est bien suprieure. Ces rsultats sont reproduits dans la Figure 11.13. e e

11.2.5

Newton amorti (damped )

Si la valeur initiale se trouve loin de la solution, la mthode de Newton et ses e variantes convergent en gnral tr`s mal. En fait la direction et surtout la longueur e e e du pas sont peu ables. An de dnir un pas plus prudent on peut calculer y (k+1) e ` litration k comme a e y (k+1) = y (k) + (k) s(k) ou (k) est un scalaire ` dterminer. Ainsi on choisira 0 < (k) < 1 lorsquon est a e (k) loin de la solution et = 1 lorsque y (k) est proche de la solution. Une faon de c contrler le param`tre (k) est de le lier ` la valeur de F (y (k)) . o e a

2 y1

Iter k 0 1 2 3 4 5 6
0 y2 3

Solution (k) (k) y1 y2 2.0000 0 2.6250 0.3750 2.4044 0.5956 3.1100 0.1100 2.9739 0.0261 2.9991 0.0009 3.0000 0.0000

Erreur (k) (k) sol sol y1 y1 y2 y2 1.0000 0 0.3750 0.3750 0.5956 0.5956 0.1100 0.1100 0.0261 0.0261 0.0009 0.0009 0.0000 0.0000

Fig. 11.13 Algorithme de Broyden avec B (0) = F (y (0) ) et y (0) = [ 2 0 ].

11.2.6

Solution par minimisation

Pour rsoudre F (y) = 0 on peut minimiser la fonction objectif suivante e g(y) = F (y)
p

o` p est une norme quelconque dans Rn . Une raison qui motive cette alternative est u quelle introduit un crit`re pour dcider si y (k+1) est une approximation meilleure e e pour y que y (k). Comme ` la solution F (y ) = 0 on peut tre tent de comparer les a e e (k+1) (k) vecteurs F (y ) et F (y ) en calculant leur norme respective. Ce que lon dsire e est que F (y (k+1)) p < F (y (k)) p ce qui conduit ` minimiser la fonction objectif a 1 min g(y) = F (y)F (y) y 2 si lon choisit p = 2 pour la norme. On utilisera lalgorithme de Gauss-Newton ou Levenberg-Marquardt.

11.3

Tableau synoptique des mthodes e

La gure 11.15 rsume les combinaisons possibles des mthodes pour la resolution e e des syst`mes dquations linaires et non-linaires. e e e e

Fig. 11.14 Graphe de

F (y)

du syst`me dquations de lexemple 11.1. e e

... ...... ....... .. .... ....... ....... ....... ....... ....... ......................................................... ... . ... .......................................................... ..... . .. ....... .. ....... . ....... ....... .. ... ....... ....... .. ..... .... ... .. ... .. .. . .. .. . ................................... . .................................... .. .. .. .. .. .. .. .. .. .. .. . .. .. .. .. ....... ....... .. .. .. .. . ..... ....... .. .. .. ....... .. ....... .. ..................................... .. ....... ............................... .. ....... ....... .. ....... .. ....... ....... .. .. ....... ....... ..... .. .... .. .. .. .. .. .. .. . . .................................... . ..................................... .. . .. . .. ... ... .. ... ... . .. . . .. .. . . .. .... . .. .. ....... . .. ....... . ....... .. .. ....... .. ....... . ....... .. .. . ..... .. ....... . ................................ .................................. .. .. . ..... .. ....... .. . .. ....... .. .. ....... .. .. ....... .. .. ....... .. .. .. .. ....... .. .... . .. .. .. .. .. . .. .. .. .. .. . . . .. .. .. .. .. . .................................... .. .. .................................. .. . .. .. .. .. .. . .. .. .. . .. .. .. .. .... . .. .. .. ....... ....... .. .. ....... .. ....... .. .. ....... .. ....... . .. ................................... .... .............................. ... . ....... .. ....... ....... ....... . .. ....... ....... . ....... .. ....... .. . .. .. .. . .. . .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . .. .. .. . .. .. .. .. .. ................................... . .. .................................... .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. . .. .. .. .. .. .. .. .. .. .......................................................................................... .......................................................................................... .. . .. . .. .. . . ... ... .. ... ... . . .. .. . . .. . .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. ................................. .................................

pleine

directe

LU Cholesky QR LU QR Jacobi Gauss-Seidel SOR Block L BiCGStab GMRES QMR

creuse

linaire e L

stationaire

itrative e

non-stat

1er ordre

Jacobi Gauss-Seidel SOR Block N

non-lin

Newton

Minimisation

Fig. 11.15 Tableau synoptique des mthodes pour la solution de syst`mes e e dquations. e

Chapitre 12 Dcompositions orthogonales e

12.1 La factorisation QR

La factorisation QR dune matrice A Rmn est donne par e

................ .. .............. .. . . . . . .. . . . . .. . . . . .. . . . .. . . . .. . . .. . . .. . .. .. . .. .. .. .

avec Q Rmm une matrice orthogonale et R Rmn une matrice triangulaire suprieure. e

Exemple 12.1 La factorisation QR dune matrice A sobtient avec la commande [Q,R]=qr(A)

1 2 A = 0 1 1 4
3 3 0 3 3 2 2 33

2 2

6 6 26 6 6 6

2 3 3 R= 0 3 0 0

et on vrie que Q Q = I. Le calcul de R ncessite 3n2 (m n et le calcul de Q ncessite e e e 3 4m2 n.

La transformation de Gauss, utilise pour la factorisation LU, est une technique e pour mettre ` zero certaines composantes dun vecteur. a Dautres techniques consistent soit ` faire subir au vecteur une rotation de sorte quil a se trouve dans la direction dune des axes du repaire orthonorm, soit doprer une e e reexion du vecteur a travers sa bisectrice ce qui lui donne a nouveau la direction dune des axes du rep`re. e 105

12.1.1

Matrices de Givens

Une rotation du vecteur x de (radians) dans le sens contraire aux aiguilles dune montre est obtenu avec la transformation Gx avec G= cos() sin() sin() cos()

et une reexion du vecteur x par rapport ` sa bisectrice [cos(/2) sin(/2)] , cest-`a a dire la droite dangle /2 entre x et le rep`re, avec langle du vecteur x, est obtenu e au moyen de la transformation F x avec F = cos() sin() sin() cos() .

La matrice G est appell matrice de rotation de Givens et la matrice F est appell e e matrice de reexion de Givens. Les matrices G et F sont des matrices orthogonales et on rapelle qune matrice F Rnn est dite orthogonale si F F = In . Les matrices orthogonales jouent un rle important dans les considrations qui suivent. o e Exemple 12.2 Soit le vecteur x = [ 3 1] et son angle (avec MATLAB = atan2(x(2), x(1))).
Alors la rotation r = Gx avec = /2 donne r = [0 2] et une rotation avec = donne s = [2 0] . Une reexion z = F x donne z = [2 0] .
2
. . . . . . . . . . .. . . . . . . . .. ... . .. .. ... . .. . . .. r . . . .. . . .. . . .. . . . . .. . 3 . . .. . . . .. . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . ..... . ..... . . .. . . . ..... . .. . . .. . .... . . . . ... ... . . . cos(/2) . . ....... .... .......... .... . . .. . .. .. .... . ... . . ........ .. . .... . ...... . .... ... . ... . . sin(/2) .... .... . . . .. .... ..... . ...... . . .. ..... . . ... . . . ... . . . ........... .. . . . . ... . .... . . . 2 . . . .. ... ....... . . . . .... ...... .. .... ....... . . . . . . .z . . ....... . ....... . . . . .. ... ................................................................................. . ..................................................................................... . .. . .. ...

s 3 2

On voit dans le graphique que la transformation orthogonale ne modit pas la longeur euclidienne du vecteur transform. On a F x = z et Gx = r . e

12.1.2

Rexion de Householder e

Soit un vecteur v Rn , alors une matrice de la forme 2vv H=I vv est appelle une matrice de Householder. On peut vrier que la matrice H est e e symtrique et orthogonale. e

Une matrice de Householder transforme un vecteur x en le reechissant ` travers a lhyperplan qui est perpendiculaire au vecteur v. Ceci est illustr dans lexemple de e gauche de la gure ??. Maintenant on voudrait choisir v de telle sorte que Hx devienne un multiple de e1 , ce qui a pour consequence dannuler toutes les composantes de x sauf la premi`re. e Pour x Rn le produit Hx scrit e Hx = (I 2vv 2v x )x = x v vv vv .

Si lon veut que Hx engendre lespace dni par e1 il faut que v soit une combinaison e lineaire de x et e1 . Posons donc v = x + e1 , alors on obtient v x = (x + e1 )x = x x + x1 et v v = (x + e1 )(x + e1 ) = x x + 2x1 + 2 et donc x x + x1 2v x Hx = (1 2 )x e1 x x + 2x1 + 2 vv .

Pour pouvoir annuler le coecient de x dans lexpression ci-dessus il faut poser = x 2 dou lon tire que v = x et on verie que Hx = (I 2 x
2

vv )x = v v

Exemple 12.3 Soit le vecteur x = [2 .5] et v = [.7 1]. A gauche est represent la e
reexion de x a travers lhyperplan perpendiculaire a v. Ensuite on rdenit v = x ` ` e x 2 e1 = [.0616 .5] et on peut verier que la reexion de x a travers lhyperplan ` perpendiculaire a v se confond avec e1 . `

12.1.3

Algorithme de Householder - Businger - Golub

Il sagit de factoriser une matrice A Rmn en A = QR o` Q est une matrice u orthogonale et R une matrice triangulaire. On cherche H1 = I
2v1 v1 v1 v1

telle que 1 0 H1 A = . . . . . . . . .

puis on cherche H2 = I

Finalement la decomposition sera donn par e

telle que 1 0 H2 H1 A = 0 0 0

2v2 v2 v2 v2

2 0 . . .

. . . 0 . . . . . .

Hn H2 H1 A = R .
Q

(TP : Golub p.231, prob. 5.3.6)

12.2

Dcomposition singuli`re SVD e e

U = [u1 , . . . , um] Rmm et V = [v1 , . . . , vn ] Rnn

Soit une matrice A Rmn alors il existe deux matrices orthogonales telles que
.. .. .. .. .. .. .. .. .. .. .. .. .. .. .. ..

avec = diag(1 , . . . , p ) Rmn et p = min(m, n). Les i sont appells les valeurs e singuli`res et elles vrient 1 2 . . . p 0. Les vecteurs ui et vi sont e e appells les vecteurs singuliers de gauche respectivement de droite. e Exemple 12.4 Avec MATLAB on obtient la dcomposition singuli`re avec la commande e e
[U,S,V]=svd(A).

12.2.1

Approximation dune matrice par une somme de matrices de rang un

Comme les matrices U et V sont orthogonales on vrie que la matrice A peut e sexprimer comme une somme de matrices de rang un, soit
A = UV = 1 u1 v1 + . . . + p up vp

Si k < r = rang(A) lapproximation de A denie comme

Ak =
i=1

i uivi

verie A Ak
2=

k+1

ce que lon peut interprter comme une msure de la qualit de lapproximation. e e e

12.2.2

Calcul numrique du rang dune matrice e

La dcomposition singuli`re dune matrice constitue un algorithme numriquement e e e tr`s stable. De ce fait cet algorithme se prete bien pour tablir le rang numrique e e e dune matrice dans une situation ou le probl`me est mal conditionn. Si les valeurs e e singuli`res dune matrice A verient e 1 2 . . . r > r+1 = . . . = p = 0 alors on a rang(A) = r. Dans MATLAB la fonction rank correspond ` : a s=svd(X); tol=max(size(X))*s(1)*eps; rang=sum(s > tol); Il existe dautres algorithmes dont la complexit est infrieure pour tablir le rang e e e dune matrice. La dcomposition singuli`re ore cependant les resultats les plus e e ables dans des situations ou le probl`me est mal conditionn. e e

12.2.3

Interprtation des valeurs singuli`res e e

Il y a un lien entre les valeurs propres de A A et AA avec les valeurs singuli`res de e A. Soit = U AV

qui est une matrice diagonale, alors = (U AV ) U AV = V A UU AV = V A AV et on conclut que les valeurs singuli`res au carr correspondent aux valeurs propres e e de A A. De faon analogue on peut crire c e = U AV (U AV ) = U AV V A U = UAA U
2 et on conclut que les i sont aussi les valeurs propres de AA .

Les valeurs propres de A A ont une interpretation gometrique. Il sagit de la demie longeur des axes de lhyperellipse x A Ax 1 = 0. Ils dcrivent comment le nuage e des points (observations) est dispers dans lespace. Le rapport de deux valeurs e singuli`res correspondantes correspond au rapport des axes de lellipse. Ainsi les e valeurs singuli`res caracteristent llongation de lhyperellipse. Ceci est illustr dans e e e la gure 12.1.

. .. .. .. ..

.. .

.. . ...... .. .. . . .. ......... . . .. .. ..

Fig. 12.1 Ellipse denie par x A Ax 1 = 0. Il sensuit que la condition dune matrice A dun syst`me linaire correspond ` e e a 2 (A) = max (A) min (A) .

12.2.4

Complexit e

La dcomposition singuli`re necessite 4m2 n+8mn2 +9n3 > 21n3 oprations lmentaires. e e e ee Il existe des versions plus conomiques de lalgorithme lorsque on ne sinteresse que e soit ` S, soit ` S et U ou soit ` S et V . On remarque que cette dcomposition est a a a e tr`s couteuse. Son intrt est nanmoins de grande importance. e ee e

12.2.5

Calcul du pseudoinverse

Soit une matrice A Rmn et U, V et les matrices qui correspondent ` sa a dcomposition singuli`re. Si lon dnit e e e A+ = V + U avec + = diag( 1 1 , , , 0, . . . , 0) 1 r

alors A+ vrie les 4 conditions de Moore-Penrose (1955). e AA+ A = A A+ AA+ = A+ (AA+ ) = AA+ (A+ A) = A+ A .

e e Exemple 12.5 Soit une matrice A et sa dcomposition siguli`re obtenu avec la commande [U,S,V]=svd(A) 2 1 6 1 0.863 0.505 0.002 A = 1 1 1 2 U = 0.286 0.485 0.827 0 1 2 3 0.416 0.714 0.563 7.304 0 0 S= 0 2.967 0 0 0 0.918 0 0 0

alors son pseudoinverse correspond a V( :,1 :r)*diag(ones(r,1)./s)*U( :,1 :r) ou ` r est le nombre de valeurs singuli`res non nulles. e 0.275 0.177 0.896 0.137 0 0 0.863 0.286 0.416 0.214 0.234 0.285 A+ = 0 0.337 0 0.505 0.485 0.714 0.862 0.376 0.339 0 0 1.089 0.002 0.827 0.563 0.367 0.879 0.041 0.061 0.788 0.576 0.015 0.303 0.106 = 0.167 0.333 0.167 0.106 0.121 0.258

0.275 0.177 0.896 0.302 0.214 0.234 0.285 0.905 V = 0.862 0.376 0.339 0.000 0.367 0.879 0.041 0.302

Le pseudoinverse peut tre obtenu aussi plus simplement avec la commande pinv(A). e

Chapitre 13 Moindres carrs e

Un probl`me fondamental en sciences consiste ` adapter un mod`le ` des obe a e a servations qui sont entaches derreurs. Nous nous intresserons ici au cas o` le e e u nombre dobservations est suprieur au nombre de param`tres. Ceci nous conduira e e ` rsoudre, suivant le choix du mod`le, des syst`mes linaires ou non-linaires a e e e e e surdtermins. e e La meilleure solution dun syst`me surdtermin peut tre dnie de plusieurs e e e e e faons. Etant donn les observations y, les variables indpendantes X et le mod`le c e e e f (X, ) avec le vecteur de param`tres, la solution retenue correspond a la solution e ` du probl`me de minimisation e min

f (X, ) y

2 2

(13.1)

La notation (13.1) est principalement utilis en statistique. En analyse numrique e e on note le probl`me comme suit : e 1 1 min g(x) = r(x) r(x) = xRn 2 2
m

ri (x)2
i=1

(13.2)

o` r(x) est le vecteur des rsidus fonction des param`tres x. u e e Dans le cas linaire le mod`le scrit Ax b avec b les observations, A les variables e e e indpendantes et x le vecteur de param`tres. Le vecteur des rsidus est dans ce cas e e e r = Ax b . Si on a un mod`le non-linaire f (ti , x), la fonction f est non-linaire par rapport aux e e e param`tres x. Les ti sont les variables indpendantes et les yi sont les observations. e e Le vecteur des rsidus scrit alors e e ri (x) = f (ti , x) yi i = 1, . . . , m.

Ce choix particulier de minimiser la norm Euclidienne des rsidus est motiv par e e deux raisons. La premi`re raison est de nature statistique car la solution correspond e 113

au BLUE (best linear unbiased estimator) dans le cas du mod`le linaire classique e e (rsidus i.i.d. et N(0, )). La deuxi`me raison est de nature numrique car la solution e e e peut sobtenir avec des procdures numriques relativement simples. e e Il semble que ce soit Gauss qui en 1795 (` lage de 18 ans) a dcouvert la mthode a e e des moindres carrs. Le dveloppement des mthodes numriques modernes pour la e e e e rsolution des moindres carrs sest fait dans les annes soixante (QR, SVD) et plus e e e rcemment on a dvelopp les mthodes pour des syst`mes qui sont grands et creux. e e e e e On prsentera dabord les mthodes pour rsoudre les moindres carrs linaires e e e e e puis on prsentera le cas non-linaire. La solution des moindres carrs non-linaires e e e e sobtient par des mthodes itratives qui ncessitent la solution dun probl`me de e e e e moindres carrs linaires ` chaque tape. e e a e

13.1

Moindres carrs linaires e e

Souvent, et cest le cas plus particuli`rement pour des probl`mes qui apparaissent e e en conomtrie et statistique, on nest pas seulement intress dobtenir la solution e e e e x du probl`me des moindres carrs Ax b, mais on dsire aussi calculer la matrice e e e des variances et covariances 2 (A A)1 . Ainsi on aura besoin dvaluer dune part b Ax 2 ` la solution, tant donn e e e 2 a 2 2 que cette valeur intervient dans le calcul de = b Ax 2 /(m n), et dautre part dune mthode ecace et numriquement stable pour calculer (A A)1 ou un e e sous-ensemble dlments de cette matrice. ee Dans la prsentation qui suit on donnera des indications comment obtenir ces e lments numriquement suivant lapproche choisi. ee e

13.1.1

Mthode des quations normales e e

La solution de (13.2) peut tre obtenue de plusieurs mani`res. Une faon consiste e e c ` driver (13.2) par rapport ` x et dcrire les conditions du premier ordre pour le a e a e minimum, soit : 2A Ax 2A b = 0 do` on tire le syst`me des quations normales u e e A Ax = A b . (13.3)

Les quations normales sobtiennent galement en considrant que la solution Ax la e e e plus proche ` b sobtient par une projection orthogonale de b dans lespace engendr a e par les colonnes de A.

b................ ........... r = b Ax . ... . ... . . . ... . ...

.. .. ... . ... . .... .. . ... .

De ce fait A r = A (b Ax) = 0 do` u A Ax = A b .

Si la matrice A est de rang complet A A est dnie positive et le syst`me des e e quations normales peut tre rsolu en utilisant la factorisation de Cholesky (cf. Thor`me 7.3). e e e e e En posant C = A A et c = A b les tapes de rsolution des quations normales A Ax = A b proc`dent dabord ` la e e e e a factorisation de Cholesky C = GG puis ` la rsolution du syst`me triangulaire infrieur a e e e Gz = c et nalement ` la rsolution du syst`me triangulaire suprieur a e e e G x = z. Factorisation de la matrice borde e Si lon dsire obtenir r e
2 2

il convient de former la matrice borde e A= Ab

et de considrer la factorisation de Cholesky de la matrice C = A A e C= C c c b b = GG

avec

G 0 z

La solution x et la norme des rsidus sobtiennent alors ` partir de e a G x = z et Ax b

= .

Pour dmontrer ce rsultat il sut de dvelopper lgalit entre les matrices A A et e e e e e GG C c G 0 G z GG Gz = = c bb z 0 z G z z + 2

do` lon voit que G correspond ` la matrice de Cholesky de A A et que u a Gz = c et b b = z z + 2 .

Comme r = b Ax est orthogonal ` Ax, a Ax do` u Ax b Calcul de (A A)1 La matrice S = (A A)1 peut tre obtenu de la faon suivante : e c S = (A A)1 = (GG )1 = G1 G1 ce qui ne ncessite que le calcul de T = G1 et de former le triangle infrieur du e e produit T T . Linverse T dune matrice triangulaire infrieure G est dni par e e tij =
3

2 2

= (r + Ax) Ax = b Ax = b AG1 z = z z
z

2 2

= x A Ax 2 b Ax +b b = b b z z = 2 .
zz zz

1/gii
i1 k=j

i=j gik tkj /gii i=j

5 e ee et ncessite n + 3 n2 6 n oprations lmentaires. Il est aussi possible de calculer e 3 2 la matrice S sans devoir inverser la matrice G (cf. Bjrck (1996, p. 119) et Lawson o and Hanson (1974, p. 6773)).

Si lon na besoin que des lments diagonaux de S, il sut de calculer le carr de ee e la norme Euclidienne des vecteurs de T , soit
n

sii =
j=i

t2 ij

i = 1, 2, . . . , n.

Les tapes de la rsolution sont rsumes dans lalgorithme qui suit : e e e e

C = A A Former C = Calculer G = solve G x = z (Calculer seulement la portion triangulaire infrieure) e C A b b A b b G z 0 (syst`me triangulaire) e

2 = 2 /(m n) S = 2T T

Calculer T = G1

(inversion matrice triangulaire)

(matrice des variances et covariances)

Lalgorithme ncessite oprations lmentaires. e e ee La mthode des quations normales est sduisante car elle repose sur des procdures e e e e tr`s simples (factorisation de Cholesky, produit matriciel, forward and backward e substitution). Dautre part la matrice initiale de dimension m n est comprime e en une matrice n n ce qui peut tre un avantage certain dans des situations avec e m n.
mthode des quations normales avec factorisation de Cholesky est : e e 1 1 2 4 10 5 2.0 0 0 1 2 1 G = 5.0 2.236 A= 0 1 3 b = 1 C = 10 30 11 5 11 7 2.5 .670 0.547 1 4 1 solve : G x = z T = G1 = x= 2.0 0.3 2 = .5472 /2 = 0.15 0.225 .075 .075 0.030

Exemple 13.1 Soit le probl`me Ax b. Alors la solution des moindres carrs par la e e =

0.500 0 1.118 0.447

S = 2 T T =

13.1.2

Solution ` partir de la factorisation QR a

En appliquant une transformation orthogonale au vecteur des rsidus on ne moe die pas la norme Euclidienne du vecteur et la minimisation conduira aux mmes e rsultats. e Pour la discussion qui suit il convient de partitioner les matrices Q et R comme suit : .. . .. . . .. .. . . .. . .. . . .. R1 .. . . .. . .. . .. . .. . . .. . .. .. .. Q1 ........ Q2 .. = . A ............... .
. . . . . . . . . . . . .

On remplace A par sa factorisation QR et on transforme le vecteur Ax b avec Q sans que ceci modie sa norme Q Q Rx Q b
I 2 2

Soit alors la matrice A Rmn et sa factorisation A = QR. Le probl`me est de e minimiser Ax b 2 . 2

On a alors R1 Q1 x b 0 Q2 que lon peut galement crire e e R1 x Q1 b

2 2

Q2 b

2 2

Du syst`me triangulaire R1 x = Q1 b on obtient la solution pour x et la somme des e rsidus au carr g(x) correspond ` Q2 b 2 . e e a 2 Dans la pratique on ne conservera pas la matrice Q mais on calcule ses colonnes au fur et ` mesure que lon en a besoin. a Dautre part on a A A = R Q QR = R R tant donn que Q est orthogonale et pour e e le calcul de linverse de A A on procdera comme on la indiqu prcdemment dans e e e e le cas de la factorisation de Cholesky. Exemple 13.2 Soit le probl`me Ax b avec A et b donns a lexemple 13.1. Alors la e e ` =
solution des moindres carrs par la factorisaton QR est : e .500 .670 .500 .223 Q1 = .500 .223 .500 .670 Q b = c = 1 .023 .547 .439 .712 Q2 = .807 .217 .392 .382 2.500 .670 0.023 .547

R1 =

2.000 5.000 0 2.236 2.0 .3

solve : R1 x = c x = et

Q b = d = 2

2 = d d/2 = 0.15

13.1.3

Dcomposition ` valeurs singuli`res SVD e a e

La dcomposition SVD est un outil particuli`rement appropri pour la rsolution e e e e du probl`me des moindres carrs. Comme prcdemment avec la factorisation QR, e e e e on exploite le fait que les transformations orthogonales ne modient pas la norme Euclidienne dun vecteur. Etant donn Ax b, avec A Rmn et rang de A gal ` r p = min(m, n), la e e a = solution du probl`me des moindres carrs peut scrire e e e x=V 1 0 r U b 0 0
A+ r

o` bien u

x=
i=1

uib vi i

(13.4)

o` les matrices V , et U sont celles de la dcomposition singuli`re de A = UV . u e e

Notons z = V x =

z1 z2

c = U b =

c1 c2

avec z1 et c1 Rr . Alors considrons la norme Euclidienne du vecteur des rsidus e e et appliquons lui une transformation orthogonale b Ax
2

= = =

U (b AV V x) c1 r 0 c2 0 0 c1 r z1 c2
2

z1 z2

Le vecteur x tel que dni par (13.4) minimise la norme Euclidienne des rsidus e e tant donn que c1 r z1 = 0. Donc SSR = m (uib)2 . Remarquons que (13.4) e e i=r+1 fournit aussi la solution lorsque A nest pas de rang plein. A partir de la dcomposition ` valeurs singuli`res de la matrice A on dduit que e a e e 1 2 (A A) = V r V o` r est la sous-matrice carr de . Il est alors possible dexu e pliciter un lment particulier de la matrice S = (A A)1 comme ee
n

sij =
k=1

vik vjk . 2 k

solution des moindres carrs par la factorisaton SVD est : e 0.219 .807 .023 0.547 5.779 0 0.383 .391 .439 .712 0 .773 = V = .322 .946 U = 0.547 .024 .807 .217 0 0 .946 .322 0.711 .441 .392 0.382 0 0 2.080 1.539 c 0.360 z1 = 1 c1 = c = U b = 1 = r 0.023 c2 1.990 0.547 x = V z1 = 2.0 .3 2 = c c2 /(m n) = 0.15. 2

Exemple 13.3 Soit le probl`me Ax b avec A et b donns a lexemple 13.1. Alors la e e ` =

13.1.4

Comparaison des mthodes e

Dans la situation ou Ax b 2 est petit et 2 (A) est grand lapproche QR est 2 suprieure ` lapproche des quations normales. e a e Dans la situation ou Ax b 2 est grand et 2 (A) est grand les deux approches 2 prouvent des dicults comparables. e e 2 (A A) = 2 (A)
2

!!!

Dans une situation avec m n lapproche des quations normales ncessite a peu e e pr`s la moiti moins doprations lmentaires et de place mmoire. e e e ee e Bien quil ne soit pas possible de dsigner un algorithme qui soit prfrable dans e ee toutes les situation, la factorisation QR est a considrer comme la mthode standard e e pour ce type de probl`me (avec Matlab on obtient cette solution avec la commande e x = A \ b). SVD produit les rsultats qui sont les plus stables numriquement. La mthode a e e e cependant un co t qui est plus lev. u e e

13.2

Moindres carrs non-linaires e e

1 1 g(x) = r(x) r(x) = 2 2
m

Dans le cas des moindres carrs non-linaires il sagit de minimiser la fonction e e ri (x)2
i=1

(13.5)

avec o` f (ti , x) est une fonction non-linaire (le mod`le) avec ti les variables indpenu e e e n dantes et x R le vecteur de param`tres ` estimer. e a An dcrire le mod`le quadratique pour la minimisation de (13.5) nous avons besoin e e des drives premi`res et secondes de g(x). La drive premi`re scrit e e e e e e e
m

ri (x) = yi f (ti , x)

i = 1, . . . , m

g(x) = avec

i=1

ri (x) ri (x) = r(x) r(x)

r1 (x) x1

(13.6)

la matrice dite Jacobienne. Le vecteur ri (x) =

r(x) =

. . .

r1 (x) xn

. . .

rm (x) x1

rm (x) xn

ri (x) . . .
ri (x) xn

correspond ` la ligne i de la matrice Jacobienne. a La drive seconde scrit e e e

g(x) = avec S(x) =

i=1

ri (x) ri (x) + ri (x) 2 ri (x) (13.7)

m 2 i=1 ri (x) ri (x).

= r(x) r(x) + S(x)

Exemple 13.4 Soit le mod`le non-linaire f (t, x) = x1 ex2 t et les donnes suivantes : e e e
t y 0.0 2.0 1.0 0.7 2.0 0.3 3.0 . 0.1

On a le vecteur des rsidus e

et la matrice Jacobienne

y1 x1 ex2 t1 y x1 ex2 t2 r(x) = 2 y3 x1 ex2 t3 y4 x1 ex2 t4

r1 (x) x1 r2 (x) x1 r3 (x) x1 r4 (x) x1 r1 (x) x2 r2 (x) x2 r3 (x) x2 r4 (x) x2

r(x) =

ex2 t2 = ex2 t3 ex2 t4

ex2 t1

x1 t1 ex2 t1

x1 t2 ex2 t2 . x2 t3 x1 t3 e x2 t4 x1 t4 e
4 x2 ti i=1 ri (x)e 4 x2 ti i=1 ri (x)x1 ti e

Les drives premi`res sont : e e e g(x) =

g(x) x1 g(x) x2

On a les m matrices de drives secondes e e 2 2 2 ri (x) =

= r(x) r(x) =

ri (x) ri (x) x1 x2 x2 1 2 ri (x) 2 ri (x) x1 x2 x2 2

et la matrice des drives secondes de g(x) est e e (ex2 ti )2 i=1 2 g(x) = 4 x1 (ti ex2 ti )2
i=1

0 ti ex2 ti

ex2 ti

x1 t2 ex2 ti i

x1 ti (ex2 ti )2 i=1 4 (x1 ti ex2 ti )2

i=1

4 i=1

ri (x)

ti ex2 ti

ti ex2 ti x1 t2 ex2 ti i

Considrons mc (x), lapproximation quadratique de g(x) dans un voisinage de xc e 1 mc (x) = g(xc ) + g(xc ) (x xc ) + (x xc ) 2 g(xc )(x xc ) 2 et cherchons le point x+ = xc + sN pour lequel mc (x+ ) = 0 ce qui est la condition du premier ordre ncessaire pour que x+ minimise mc . On a e mc (x+ ) = g(xc ) + 2 g(xc ) (x+ xc ) = 0
sN

o` sN est le pas de Newton que lon determine en rsolvant le syst`me linaire u e e e 2 g(xc ) sN = g(xc ) . D`s lors on pourrait envisager la solution du probl`me (13.5) avec la mthode de e e e Newton dont litration k est dnie comme e e Solve 2 g(x(k) ) sN = g(x(k) ) x(k+1) = x(k) + sN
(k) (k)

Dans la pratique on proc`de cependant diremment tant donn que lvaluation e e e e e de S(x) dans (13.7) peut savrer tr`s dicile, voire impossible. e e Les mthodes qui seront prsentes ci-apr`s se direncient par la faon dont ils e e e e e c 2 approchent la matrice des drives secondes g(x). e e

13.2.1

Mthode de Gauss-Newton e

Cette mthode utilise une approximation de la matrice des drives secondes (13.7) e e e en omettant le terme S(x). Comme S(x) est compos dune somme de termes e 2 ri (x) ri (x) cette simplication se justie dans une situation ou les rsidus ri (x) e sont petits. Cest une situation ou le mod`le pouse bien les donnes. La mthode e e e e de Gauss-Newton se rsume dans lalgorithme qui suit. e Algorithme 23 Mthode de Gauss-Newton pour les moindres carrs non-linaires e e e
1: Choisir x(0) 2: for k = 0, 1, 2, . . . until convergence do 3: Calculer r(x(k) ) 4:

Solve

5: Update 6: end for

x(k+1) = x(k) + sGN

r(x(k) ) r(x(k) ) sGN = r(x(k) ) r(x(k) )

(k)

On remarque que le syst`me linaire qui dnit le pas de Gauss-Newton sGN est un e e e syst`me dquations normales et que le pas est aussi la solution dun probl`me des e e e moindres carrs. En accord ` ce qui a t dit prcdemment au sujet des moindres e a ee e e carrs linaires on prfra la solution du syst`me surdtermin e e ee e e e r(x(k) ) sGN r(x(k) ) obtenu ` partir de la factorisation QR ` la solution des quations normales. a a e Lalgorithme peut ne pas converger si le point de dpart est choisi trop loin de la e solution. Lorsque les rsidus sont grands au point de la solution, lapproximation de la matrice e des drives secondes peut savrer insusante, avec comme consquence soit une e e e e convergence tr`s lente ou pas de convergence du tout. e
(k)

(k)

13.2.2

Mthode de Levenberg-Marquardt e

Lorsque la mthode de Gauss-Newton choue, notamment lorsque la solution du e e sous-probl`me de la solution du syst`me linaire nest pas de rang complet, la e e e mthode de Levenberg-Marquardt constitue une alternative intressante. Dans cette e e mthode on approche la matrice S(x) par une matrice diagonale I. Nous avons alors e lalgorithme suivant : Algorithme 24 Mthode de Levenberg-Marquardt pour les moindres carrs none e linaires e
1: Choisir x(0) 2: for k = 0, 1, 2, . . . until convergence do 3: Calculer r(x(k) ) et k

Solve

5: Update 6: end for

x(k+1) = x(k) + sLM

r(x(k) ) r(x(k) ) + k I sLM = r(x) r(x)

(k)

que lon rsoudra en utilisant la factorisation QR. Ainsi on vite le calcul du produit e e r(x(k) ) r(x(k) ) qui peut introduire des instabilits numriques. e e k est choisi ` partir des considrations sur la trust region. Dans la pratique on a e 2 choisit la valeur 10 . Suivant le choix de k , 0 k < le pas de Levenberg-Marquardt se situe entre le pas de Gauss-Newton pour k = 0 et un pas steepest descent. Si k est choisi judicieusement la mthode de Levenberg-Marquardt sav`re tr`s e e e robuste en pratique. Elle constitue lalgorithme de base dans un grand nombre de logiciels spcialiss. e e Lalgorithme peut ne pas converger si le point de dpart est choisi trop loin de la e solution. Lorsque les rsidus sont grands au point de la solution, lapproximation de la matrice e des drives secondes peut savrer insusante, avec comme consquence soit une e e e e convergence tr`s lente ou pas de convergence du tout. e

Dans ce cas la solution du pas sLM est dnie ` partir du probl`me des moindres e a e carrs linaires e e (k) (k) r(x ) (k) r(x ) sLM 1/2 k I 0

(k)

Bibliographie
Bjrck, A. (1996). Numerical Methods for Least Squares Problems. SIAM. Philadelo phia, PA. Dulmage, A.L. and N.S. Mendelsohn (1963). Two algorithms for bipartite graphs. SIAM 7, 183194. Gilli, M. (1995). Graph-Theory Based Tools in the Practice of Macroeconometric Modeling. In : Methods and Applications of Economic Dynamics (S. K. Kuipers, L. Schoonbeek and E. Sterken, Ed.). Contributions to Economic Analysis. North Holland. Amsterdam. Gilli, M. and M. Garbely (1996). Matching, Covers, and Jacobian Matrices. Journal of Economic Dynamics and Control 20, 15411556. Golub, G. H. and C. F. Van Loan (1989). Matrix Computations. Johns Hopkins. Baltimore. Lawson, C. L. and R. J. Hanson (1974). Solving Least Squares Problems. PrenticeHall. Englewood Clis, NJ. Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vetterling (1986). Numerical Recipes. Cambridge University Press. Cambridge.

125

Index
limination e de Gauss, 31, 32, 35 de Gauss-Jordan, 45 algorithme Broyden, 99 Gauss-Newton, 122 Gauss-Seidel, 69 Jacobi, 69 Levenberg-Marquardt, 123 Newton, 96 SOR, 74 bsub, 28 cholesky, 53 egpp, 42 eg, 37 fsub, 28 ldl, 51 ldm, 49 lu3diag, 55 perm, 40 rsl, 43 tgm1, 34 fsub1, 43 back substitution, 28 base, 2 catastophic cancellation, 10 chopping, 4 classication des algorithmes, 20 complexit e dun algorithme, 18 non-polynmiale, 20 o polynomiale, 20 condition dun probl`me, 11 e dune matrice, 12, 13 convergence linaire, 97 e quadratique, 97 super-linaire, 97 e taux de, 97 digits signicatifs, 8 eps, 6 erreur absolue, 10 darrondi, 5 de troncature, 1 rsiduelle, 14 e relative, 9 exposant, 2 factorisation de Cholesky, 52 LDL, 50 LDM, 47, 48 LU, 31 ops, 20 forward substitution, 27 Gauss limination de, 32 e matrice dlimination de, 33 e multiplicateur de, 33 Gauss-Jordan, limination de, 45 e Gauss-Seidel non-linaire, 92 e Inf, 4 instabilit numrique, 11 e e itration e de Gauss-Seidel, 68 de Jacobi, 68 Jacobi 126

non-linaire, 92 e mantisse, 2 normalisation, 2, 3 Matlab code FPI, 88 bracketing, 87 converged, 75 matrice dchange, 39 e dlimination de Gauss, 33 e dnie positive, 51 e de Hessenberg, 24 de permutation, 38 de Tplitz, 24 o par bande, 55 transformation de Gauss, 33 triangulaire unitaire, 28 tridiagonale, 55 multiplicateur de Gauss, 33 Newton, 95 amorti, 100 convergence, 98 Newton-Raphson, 95 oset, 4 oprations lmentaires, 20 e ee ordre dune fonction, 19 overow, 4 perfect rounding, 4 pivot, 33, 35 pivotage complet, 41 partiel, 41, 42 point xe, 93 prcision machine, 5, 6 e probl`me mal conditionn, 11, 12 e e Quasi-Newton, 99 rang structurel, 61 rayon spectral, 72 roundo error, 9 SOR, 73

param`tre de relaxation, 73 e sur-relaxation successive, 73 syst`me e triangulaire infrieur, 27 e triangulaire suprieur, 28 e underow, 4 virgule ottante, 13, 8 ensemble des nombres en, 3 oprateur, 4 e

Vous aimerez peut-être aussi

Erreurs en calcul numérique
Pas encore d'évaluation
Erreurs en calcul numérique
132 pages
OptiAlgo PDF
Pas encore d'évaluation
OptiAlgo PDF
53 pages
Opti Algo
Pas encore d'évaluation
Opti Algo
50 pages
Analyse numérique : erreurs et méthodes
Pas encore d'évaluation
Analyse numérique : erreurs et méthodes
247 pages
Analyse numérique : cours et méthodes
Pas encore d'évaluation
Analyse numérique : cours et méthodes
261 pages
Calcul Matriciel
100% (1)
Calcul Matriciel
83 pages
AnaNumI - 15 09 2020
Pas encore d'évaluation
AnaNumI - 15 09 2020
229 pages
Cours MN
Pas encore d'évaluation
Cours MN
47 pages
Arithmétique des ordinateurs et erreurs
Pas encore d'évaluation
Arithmétique des ordinateurs et erreurs
170 pages
Anlyse M3
Pas encore d'évaluation
Anlyse M3
228 pages
Cours MN
Pas encore d'évaluation
Cours MN
47 pages
Polycope AN1 2021 2022
Pas encore d'évaluation
Polycope AN1 2021 2022
120 pages
Travaux Pratiques D'analyse Numérique Avec Python
Pas encore d'évaluation
Travaux Pratiques D'analyse Numérique Avec Python
25 pages
Livre
Pas encore d'évaluation
Livre
442 pages
Livre
Pas encore d'évaluation
Livre
442 pages
Cours de Calcul Numérique à Polytech'Lille
Pas encore d'évaluation
Cours de Calcul Numérique à Polytech'Lille
107 pages
Support Cours PHY200
Pas encore d'évaluation
Support Cours PHY200
34 pages
Notes Dec Ours
Pas encore d'évaluation
Notes Dec Ours
95 pages
Cours
100% (1)
Cours
140 pages
Analyse Numérique et Optimisation
Pas encore d'évaluation
Analyse Numérique et Optimisation
150 pages
SUPPORT DE COURS Introduction Au Calcul Scientifique
Pas encore d'évaluation
SUPPORT DE COURS Introduction Au Calcul Scientifique
52 pages
Méthodes Numériques Systèmes Linéaires
100% (1)
Méthodes Numériques Systèmes Linéaires
72 pages
Cours d'Analyse Numérique 2019/2020
Pas encore d'évaluation
Cours d'Analyse Numérique 2019/2020
84 pages
Analyse Matricielle LAKEHALI PDF
Pas encore d'évaluation
Analyse Matricielle LAKEHALI PDF
72 pages
Analyse Numérique - EHTP
Pas encore d'évaluation
Analyse Numérique - EHTP
85 pages
Un Cours D'Analyse Numérique
Pas encore d'évaluation
Un Cours D'Analyse Numérique
168 pages
Cours Analyse Numérique
Pas encore d'évaluation
Cours Analyse Numérique
57 pages
Analyse Numerique
Pas encore d'évaluation
Analyse Numerique
59 pages
Méthodes Numériques pour Systèmes Linéaires
Pas encore d'évaluation
Méthodes Numériques pour Systèmes Linéaires
68 pages
Notes de cours sur les méthodes numériques
Pas encore d'évaluation
Notes de cours sur les méthodes numériques
94 pages
2012 2013 M33 L2
100% (1)
2012 2013 M33 L2
176 pages
Résolution Numérique de Systèmes Linéaires
Pas encore d'évaluation
Résolution Numérique de Systèmes Linéaires
72 pages
COURS Analyse Numérique Nantes PDF
Pas encore d'évaluation
COURS Analyse Numérique Nantes PDF
135 pages
Calcul Formel Et Numerique in French Handbook of The Numerical Analysis Course of ULB Computer Science Department
Pas encore d'évaluation
Calcul Formel Et Numerique in French Handbook of The Numerical Analysis Course of ULB Computer Science Department
166 pages
Analyse Numérique : Méthodes et Applications
Pas encore d'évaluation
Analyse Numérique : Méthodes et Applications
135 pages
Analyse Numerique Saad
100% (1)
Analyse Numerique Saad
135 pages
Traitement numérique des données
100% (3)
Traitement numérique des données
38 pages
Analyse Numérique : Concepts et Méthodes
100% (1)
Analyse Numérique : Concepts et Méthodes
62 pages
Guide Complet sur l'Analyse Numérique
Pas encore d'évaluation
Guide Complet sur l'Analyse Numérique
168 pages
An Numerique
Pas encore d'évaluation
An Numerique
168 pages
Méthodes Numériques et Systèmes Linéaires
Pas encore d'évaluation
Méthodes Numériques et Systèmes Linéaires
73 pages
Analyse numérique et arithmétique flottante
Pas encore d'évaluation
Analyse numérique et arithmétique flottante
87 pages
COURS MethodesNumeriques2015 - SOUK - AHRAS PDF
0% (1)
COURS MethodesNumeriques2015 - SOUK - AHRAS PDF
125 pages
TD Et Cours
Pas encore d'évaluation
TD Et Cours
125 pages
Cours d'Analyse Numérique en Python
Pas encore d'évaluation
Cours d'Analyse Numérique en Python
93 pages
Cours Anal Num31-10-2018
Pas encore d'évaluation
Cours Anal Num31-10-2018
37 pages
Introduction aux langages numériques
Pas encore d'évaluation
Introduction aux langages numériques
131 pages
Algèbre Linéaire et Analyse Matricielle
100% (2)
Algèbre Linéaire et Analyse Matricielle
167 pages
Blog
Pas encore d'évaluation
Blog
1 page
Révolution de l'entraînement au foot avec SenseBall
Pas encore d'évaluation
Révolution de l'entraînement au foot avec SenseBall
2 pages
Introduction au langage Java et concepts clés
Pas encore d'évaluation
Introduction au langage Java et concepts clés
59 pages
Les Principes Fondamentaux Des Echecs
100% (6)
Les Principes Fondamentaux Des Echecs
124 pages
Méthodes de Conception Orientées Objet
Pas encore d'évaluation
Méthodes de Conception Orientées Objet
62 pages
Cours MCSI PDF
Pas encore d'évaluation
Cours MCSI PDF
38 pages
Droit des TIC et Enjeux Sociaux
Pas encore d'évaluation
Droit des TIC et Enjeux Sociaux
22 pages
09 Complex
Pas encore d'évaluation
09 Complex
11 pages