UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Chapitre-III : Méthode de résolution directe des systèmes d’équations linéaires

III.1 Introduction :

L’objectif de ce chapitre est une introduction aux méthodes numériques pour la résolution des systèmes
linéaires A.X=B. La plupart du temps, A sera réelle, mais pas toujours carré.

On distingue deux grandes familles de méthodes :

- Les Méthodes directes : On obtient la solution en un nombre fini d’opérations. Si elle est
numériquement « mauvaise », à cause des erreurs d’arrondis par exemple, on ne peut en général plus
rien faire pour l’améliorer. Par contre, on a l’assurance de « Terminer ».

- Les Méthodes itératives : On construit une suite qui tend vers la solution. Si on n’est pas satisfait du
résultat fourni, on peut « continuer ». Cette fois, la qualité de la convergence est primordiale.

III.2 Description : On se propose de résoudre un système linéaire à matrice carrée n (n équations à n

inconnues).

  +   + … … … . +  = 

  +   + … … … . +  = .

..
   +    + … … … . +  =

Ou encore :

'   =  ) = 1, 2, … , ,
(

Ou sous forme matricielle.


⎡  ⎤ ⎡ ⎤


⎢ . ⎥ ⎢ ⎥
.
. =    .  ⎢ . ⎥ =  ⎢ ⎥
⎢ ⎥ ⎢ . ⎥
,

⎢ ⎥. ⎢ . ⎥
⎣ ⎦ ⎣ ⎦

On suppose que la solution existe et qu’elle unique, d’où : ∆= det(A) ≠ 0

La solution d’un tel système est donnée comme suit :

∆
 = ; ∆= det(A)
∆

22
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

∆i : est le déterminant de A après avoir remplacé la colonne i dans A par le membre constant B.

Cette méthode est appelée méthode de Cramer, c‘est une méthode très connue pour la résolution des
systèmes linéaires. Toutefois, elle est la moins recommandable quand n est relativement grand (n > 8).

III.3 Système linéaire à matrice particulière

III.3.1 Système à matrice A triangulaire supérieur : Soit à résoudre le système suivant :

  … … 
⎡  … …  ⎤
⎢ . ⎥
.  =  ù =⎢ ⎥
.
⎢  ⎥
⎣ ⎦

A est une matrice triangulaire supérieure.

On a :

'   =  ) = 1, 2, … , ,
(

Où :  = 0 .) / )

 = '   = '   +   ) = 1, 2, … , ,

( (0

xi : S’obtient donc aisément par le calcul suivant :

1
 = 1 − '   3 ) = ,, , − 1, … ,1
 
(0

III.3.2 Système à matrice A diagonale : Soit le système A.X=B, où A est une matrice diagonale. La solution
X s’obtient immédiatement par :

 = ) = 1, 2, … , ,



III.4 Système à matrice A quelconque (Méthode directe)

A) Méthode de Gauss.

23
UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Principe : La méthode de Gauss consiste à transformer le système A.X=B à matrice quelconque en un

système équivalent A’ . X=B’ où A’ est une matrice triangulaire supérieure.

(4) Ÿ (4 5 )

  +   + 7 7 = (1)

Soit à résoudre le système suivant :

6  +   + 7 7 =  (2)
7  + 7  + 77 7 = 7 (3)

Etape1 :

Elimination de  de l’équation (2) et (3) :

de (1) nous aurons :  = 9 ( −  .  − 7 . 7 )



::

  +   + 7 7 = 

⎧
⎪ 0 + ? −  .  @  + ? −  . 7 @  =  . 
−
 
  7
 7 

⎨  .  . 7 . 
⎪ 0 + ?7 − 7  @  + ?77 − 7 7 @ 7 = 7−
⎩   
 .  7 . 
⎧A −  B =  (1) A7 − (1) ⎫
 B = 7
⎪  . C  . ⎪
A7 −  7 B = 7 (1) A77 −  7 B = 77 (1)
 C 
⎨ ⎬
⎪ A −  .  B =  .
(1) A 7 − 7  B = 7 (1) ⎪
⎩     ⎭

  +   + 7 7 =  (1′)

6 0.  +  ()  + 7 () 7 = 
()
(2′)
0.  + 7 ()  + 77 () 7 = 7
()
(3′)

Etape2 :

Elimination de  de la 3ème équation :

de (2) nous aurons :  =  ( − 23 (1) . 7 )

 (1)
(1) 2
22

  +   + 7 7 =

⎧ 
⎪ 0.  +   + 7 7 =
() ()

()

⎨ 32 . 23 32 (1) .

(1) (1) ()
⎪0.  + 0.  + H33 − I 3 =

−
(1) ()
⎩ 22 (1) 22 (1)
7

7 () . 7 () 7 () . 2

(1)
J77 () − = 77 (2) (2)K − = (2)L
(1) (2)
 () 3
 () 7

24
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

  +   + 7 7 = 

6 0.  +  ()  + 7 () 7 = 
()

0.  + 0.  + 77 () 3 = 7 ()

Un système (4 5 ) équivalent au système initial et à matrice triangulaire supérieur.

Résumé :

Etape1 :

 (M) .  (M)

)/ (1) = H (M) − I ; ) = 2,3
 (M)
 (M) . 1
(0)
= − ; / = 2,3
(1) (0)
 )
 (M)

Etape2 :

 () .  ()

)/ (2) = H () − I ; ) = 3,3
 ()
 () . 2
(1)
= − ; / = 3,3
(2) (1)
 )
 ()

La formule générale pour la transformation des  NO  d'une étape à une autre par la méthode de GAUSS
dans le cas général est:

RSU (UVW) . RUT (UVW)

RST (U) = HRST (UVW) − I; U = W, . . . , X − W
RUU (UVW)
RSU (UVW) . YU
(UVW)
YS (U) = YS − ; S = U + W, . . . , X ; T = U + W, . . . , X
(UVW)
RUU (UVW)

Notion de Pivot: Dans la pratique, il se peut que l'un des pivots de la méthode de Gauss soit nul. Dans ce cas,
on cherche dans la même colonne s'il y a une ligne où le pivot n'est pas nul, et on fait une permutation entre
les deux lignes. On ne choisit pas n'importe quelle ligne. En effet, dans un ordinateur, plus on divise par des
membres petits, plus on risque de faire des erreurs d'arrondis importantes. On choisit donc le nouveau pivot
de manière à minimiser ces erreurs ; on prend celui dont la valeur absolue est la plus grande.

Tabulation de la méthode de GAUSS.

Exemple: Soit à résoudre le système suivant:

2 +  − 27 = 1
P  +  + 37 = 4
 + 2 − 7 = −1

Etape (0) :

25
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

2 1 −2 1
Z1 1 3 4[
1 2 −1 −1

Etape (1) :

2 1 −2 1
Z 0 1/2 4 7/2 [
0 3/ 2 0 − 3/2

Etape (2) :

2 1 −2 1
Z0 1/2 4 7/2[
0 0 − 12 − 12

Après la résolution on trouve : 7 = 1;  = −1;  = 2

Remarque: Si A a subit des permutations de lignes (cas d'un pivot nul) et leurs nombres est P.

fNO( ) = fNO( 5 ) = (−1)g . h 5

(

Exercice : Soit à résoudre le système suivant par la méthode de Gauss.

 + 3 + 37 = 0
2
P  + 2 = 4
3 + 2 + 67 = 11

Algorithme :

I) Soit à résoudre _. ` = a (_X,X , aX )

II) Varier k de 1 jusqu'à n-1 (n-1 étapes)

II.1) Vérification du pivot

Si bb = 0

Alors

Chercher cb ≠ 0 (d = e + 1, . . . , ,)

Si cb ≠ 0 existe

Alors

On permute entre la ligne L et la ligne k.

Sinon

Ecrire (" Déterminant nul, pas de solution unique")

Finsi

Finsi

26
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

II.2) Varier i de k+1 jusqu'à n (n lignes)

II.2.1) Calculer :b = 9 ij

9
jj

II.2.2) Varier le j de k+1 jusqu'à n (colonnes)

II.2.2.1)  =  − Ωb . b

II.2.2.2)  =  − Ωb . b

III) Calcul des solutions Xi

III.1) Varier le i = n jusqu'à 1 (Remonter le système).

III.1.1) S ←
YS
RSS

III.1.2) varier le j de i+1 à n (colonnes)

III.1.2.1) S ← −
(RST . mn )
 RSS

IV) Afficher les solutions.

Remarque :

- Pour un système d'ordre n=10

Cramer ≈ 3.109 opérations

Gauss ≈ 805 opérations

- Pour un système d'ordre n=50

Cramer ≈ 7.1067 opérations

Gauss ≈ 9.104 opérations

B) Méthode de Jordan (Gauss-Jordan).

Principe : La méthode de Jordan consiste à transformer le système Cramérien A.X=B en un système

A’ . X=B’ où A’ est une matrice identité d'ordre n.

Cette diagonalisation est opérée en n étapes qui se composent d'une opération de normalisation suivie
d'une opération de réduction.

Etape1 : Elimination de  de toutes les équations sauf de la première ligne (même chose que Gauss).

Etape2 : Elimination de  de toutes les équations sauf de la deuxième ligne (y compris la première
équation).

Etape (n-1) : Elimination de  V de toutes les équations sauf de l'équation (n-1).

Etape (n) : Elimination de  de toutes les équations sauf de la dernière.

27
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Le développement de l'algorithme de Jordan pour un système d'ordre n d'une étape (k-1) à une étape
(k) est comme suivant:

RST (U) = ARST (UVW) − ΩSU

(UVW)
. RST (UVW) B ; U = W, . . . , X
YS (U) = YS − ΩSU . YU ; S = W, . . . , X (S ≠ U); T = U + W, . . . , X
(UVW) (UVW)

Tabulation de la méthode de Jordan:

Exercice : Soit à résoudre le système suivant par la méthode de Jordan.

1 3 3  0
Z2 2 0 [ x Zp [ = Z2[
3 3 6 q 9

 = 3 ;  = −2 ; 7 = 1

Etape (0) :

1 3 3 0
Z2 2 0 2[
3 3 6 9

Etape (1) :

1 3 3 0
Z0 −4 −6 2[
0 −6 −3 9

Etape (2) :

28
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Normalisation Réduction

⎡1 0 −
7 7
1 3 3 0 ⎤
⎢ ⎥
10 1 − 3 ⎢0 1 − ⎥
7  7
  
0 −6 −3 9 ⎢ ⎥
Ÿ

⎣0 0

 ⎦

Etape (3) :

Normalisation Réduction

1 0 −
7 7
  1 0 0 3
s0 1
7
− t
 Z0 1 0 − 2[
  0 0 1 1
Ÿ
0 0 1 1

1) Résolution simultanée de m systèmes à matrice A identique.

Supposons que l’on parte avec la matrice u =  , ()

, ()
,…, (v)
 ù ()
, ()
,…, (v)
sont les m
vecteurs à n dimensions.

La méthode de Jordan donne :

w= V
.u = V
 , ()
, ()
,…, (v)


w= V
. u = x, V
. ()
, V
. ()
,…, V
. (v)


Les m vecteurs  V
. ()
, V
. ()
,…, V
. (v)
 sont les solutions respectives des m systèmes linéaires
suivants :

.  () = ()
;
⎧
⎪ . = ;
() ()

.
⎨ .
⎪ .
⎩ .  (v) = (v)

2) Calcul de la matrice inverse A-1 par la méthode de Jordan.

Soit à calculer les solutions des systèmes suivants :

1 11 12 … … 1, 1
⎡0⎤ ⎡ ⎤
⎢ ⎥ ⎡ 21 22 … … 2, ⎤ ⎢0⎥
. ⎢ . ⎥ ⎢.⎥
.  =  ù  = ⎢ ⎥ V
=⎢ ⎥  ⎢.⎥
⎢.⎥ ⎢ . ⎥
⎢.⎥ ⎦ ⎢.⎥
;

⎣ ,1 ,2 ,,

⎣0⎦ ⎣0⎦

.  =  Ÿ  = V
. 

29
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique


⎡  ⎤
⎢ . ⎥
= V
.  = ⎢ . ⎥ La solution du système .  =  donne la première colonne de la matrice A-1.
⎢ ⎥
⎢ . ⎥
⎣  ⎦

0 
⎡1⎤ ⎡ ⎤
⎢ ⎥ ⎢ . ⎥
.
.  =  ù  = ⎢ ⎥ ; .  =  Ÿ  = V
.  = ⎢ . ⎥
⎢.⎥ ⎢ ⎥
⎢.⎥ ⎢ . ⎥
⎣0⎦ ⎣  ⎦

La solution du système .  =  donne la deuxième colonne de la matrice A-1.

0 
⎡0⎤ ⎡  ⎤
⎢ ⎥ ⎢ . ⎥
.
. =  ù  = ⎢ ⎥ ; . =  Ÿ  = V
. = ⎢ . ⎥
⎢.⎥ ⎢ ⎥
⎢.⎥ ⎢ . ⎥
⎣1⎦ ⎣ ⎦

La solution du système .  =  donne la nème colonne de la matrice A-1.

.  =  donne la 1er colonne de A-1 ;

.  =  donne la 1er colonne de A-1 ;

.  =  donne la nème colonne de A-1 .

En résolvant les n systèmes on trouvera la matrice A-1.

identité x , .
Comme les n systèmes ont la même matrice A on utilisera la matrice A augmentée en lui adjoignant la matrice

Où : x , = z | | … … | ]

On aura : u = z , x].

Après les opérations de l’algorithme de Jordan, on obtiendra :

w= V
.u = V
z , x] ; w = zx, V
]

Résumé :

} = z_, ~] Ÿ  = z~, _VW ]

30
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

par ce dernier, c’est-à-dire diviser la ligne k à l’étape k par l’élément b,b .

Remarque : L’étape de normalisation dans la méthode de Jordan consiste à diviser, à chaque ligne du pivot

Application : Calculons l’inverse de A:

=? 7 V@
  V€
7
€  V

Etape (0) :

2 1 −4 1 0 0
Z3 3 −5 0 1 0[
4 5 −2 0 0 1

Etape (1) :

Normalisation Réduction

1 −2 0 0
 
1 −2 0 0
 
 
13 03 s0 0t
 
3 −5 0 1 1 − 1
7 7
 
4 5 −2 0 0 1
Ÿ
0 3 6 −2 0 1

Etape (2) :

Normalisation Réduction

1 −2 0 0 1 0 −7 1 −7 0
  ƒ 
 
s0 1 −1
 
0t s0 1

−1 7

0t
7 7 7
Ÿ
0 3 6 −2 0 1 0 0 4 1 −2 1

Etape (3) :

Normalisation Réduction

⎡1 0 −7 1 −7 0⎤ ⎡1 0 0  −  ⎤
ƒ  „ 7 ƒ

⎢ ⎥ ⎢ ⎥
⎢ 0 1 −1 0⎥ ⎢1 1 0 − 1 − ⎥
  ƒ
7 7
⎢  ⎥ ⎢  ⎥
Ÿ

⎣ 0 0 1 − ⎣ 0 0 1 € − € ⎦
   
€  € ⎦

31
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

1 19 − 18 7
_VW = . Z−14 12 −2[
12
3 −6 3

Exercice : Soit à résoudre le système suivant A.X=B et trouver la matrice inverse de A.

3 + 5 + 27 = 1

P3 + 4 − 7 = 5
4 − 4 + 27 = 0

 = 0,9149 ;  = 0,2128 ; 7 = −1,4043

−0,0426 0,1915 0,1383

_VW = Z0,1064 0,0213 − 0,0957 [
0,2979 − 0,3404 0,0319

III.5 Méthode de factorisation LU à matrice A quelconque

En générale, les systèmes linéaires sont issus de problèmes physiques complexes (équations aux
dérivées partielles ou équations différentielles) et la dimension de la matrice est souvent très grande. Par
conséquent, la résolution peut avoir un coup de calcul plus ou moins important suivant la nature/structure de
la matrice. Pour cette raison, il existe différentes méthodes numériques dans la littérature. Nous allons nous
intéresser à la factorisation LU d’une matrice A afin de résoudre un système Ax = b.

III.5.1 Principe de la méthode

triangulaires L et U. Comment cela nous permet-il de résoudre le système ⃗ = ˆ⃗? Il suffit de remarquer que :
Supposons un instant que nous ayons réussi à exprimer la matrice A en un produit de deux matrices

⃗ = dŠ⃗ = ˆ⃗

Et de poser Š⃗ = ‹⃗. La résolution du système linéaire se fait alors en deux étapes :

d‹⃗ = ˆ⃗

Š⃗ = ‹⃗

‹⃗ et par la suite une remontée triangulaire sur la matrice U pour obtenir la solution recherchée ⃗.
Qui sont deux systèmes triangulaires. On utilise d’abord une descente triangulaire sur la matrice L pour obtenir

Il faut tout de suite souligner que la décomposition LU n’est pas unique. On peut en effet écrire un
nombre réel comme le produit de deux autres nombres d’une infinité de façons. Il en est de même pour les
matrices.

 = Š V ‹ = Š V dV = V

32
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

a) Description des méthodes de factorisation

En développant l’équation A=LU, on obtient :

 = ∑b( b . b

 = 0 .) e > )
Où, par définition : ŽŠb = 0 .) e > /
b

D’où :

 = ∑b( b . b

‘’(,)

La partie triangulaire supérieure de A a pour termes :

“ = ∑“b( “b . b j = r, r+1,…, n

La partie triangulaire inférieure de A a pour éléments :

“ = ∑“b( b . b“ i = r, r+1,…, n

D’où :

“ =
9”n V∑”–:
j—: •”j .˜jn
•””
j = r, ……..,n

“ =
9i” V∑”–:
j—: •ij .˜j”
˜””
i = r, ……..,n

Ces équations étant vraie, ∀r =1,2,…,n.

Supposons que l’on connaisse les n2 éléments (aij) et que l’on cherche les (n2+n) éléments non nuls des
matrices L et U.

Les systèmes (aij) ou (urj + lir) sont donc de n2 équation à (n2+n) inconnues. Il y a donc n paramètres à fixer
arbitrairement, dans le système (urj + lir), on pourra par exemple fixer les éléments diagonaux urr ou lrr des
matrices U ou L. On débouche sur les deux algorithmes les plus connus.

- algorithme de Doolittle : lii = 1, ∀i ;

- algorithme de Crout : uii = 1, ∀i ;

Les algorithmes s’écrivent donc :

• Algorithme de Factorisation : (A → LU ) de Doolittle.

33
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

lii =1 i = 1,………..,n

“ = “ − ∑“V

b( “b . b j = r,………...,n

r =1,.........,n

“ = “ − ∑“V

b( b . b“ /““ i = r+1,…….,n

Application : Soit à résoudre le système d’équations linéaires suivant, avec la méthode de Doolittle:

.  =  … … (1)
+‹+™ =5 ⎧
⎪ = dŠ … … (2)
P + 2‹ + 2™ = 6 dŠ.  = 
 + 2‹ + 3™ = 8 ⎨ d. ‹ =  … … (3)
⎪
⎩ Š.  = ‹ … … (4)

1ère méthode (Factorisation de la matrice A selon l’algorithme de Doolittle):

1 1 1
= Z1 2 2[
1 2 3

Š =  −  . Š Ÿ Š =  Ÿ Š = 1 , avec :  =  = 77 = 1
r=1

Š =  −  . Š Ÿ Š =  Ÿ Š = 1

Š7 = 7 −  . Š7 Ÿ Š7 = 7 Ÿ Š7 = 1

 = ( −  . Š )/Š Ÿ  =  /Š Ÿ  = 1/1 = 1
----------------------------------------------------------------------------------

7 = (7 − 7 . Š )/Š Ÿ 7 = 7 /Š Ÿ 7 = 1/1 = 1

Š =  −  . Š Ÿ Š = 2 − 1.1 Ÿ Š = 1

r=2

Š7 = 7 −  . Š7 Ÿ Š7 = 2 − 1.1 Ÿ Š7 = 1

7 = (7 − 7 . Š )/Š Ÿ 7 = (2 − 1.1)/1 Ÿ 7 = 1/1 = 1

----------------------------------------------------------------------------------

Š77 = 77 − 7 . Š7 − 7 . Š7

r=3

Š77 = 3 − (1.1) − (1.1) Ÿ Š7 = 3 − 2 = 1

323

Avec la factorisation, notre matrice A de départ est transformée en deux matrices comme suit :

1 0 0 1 1 1
d = Z1 1 0[ Š = Z0 1 1[
1 1 1 0 0 1
,

34
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

2ème méthode :

Etape (1) :

1 1 1  5
Z 1 2 2 [ . š‹› = Z6[
1 2 3 ™ 8

A . x = B

Etape (2) :

1 0 0 f N 
= dŠ d = Z 1 0[ Š = Z0 ž ℎ[
œ 1 0 0 )
,

1 0 0 f N 
= Z 1 0 [ . 1 0 ž ℎ3
œ 1 0 0 )

f N  1 1 1
= 1f N + ž  + ℎ 3 = Z 1 2 2 [
f N + œž  + œℎ + ) 1 2 3

f=1 N=1 =1

= 1f = 1 ( = 1) N + ž = 2 (ž = 1)  + ℎ = 2 (ℎ = 1) 3
f = 1 ( = 1) N + œž = 2 (œ = 1)  + œℎ + ) = 3 () = 1)

Avec la factorisation, notre matrice A de départ est transformée en deux matrices comme suit :

1 0 0 1 1 1
d = Z1 1 0[ Š = Z0 1 1[
1 1 1 0 0 1
,

Etape (3) :

Déterminant d’abord y par la résolution de :

/
d. ‹ =  = Ze [

et

1 0 0 / 5
Z1 1 0[ . Ze[ = Z6[
1 1 1  8

On trouve :

35
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

/=5 / 5
6e = 6 − 1. (5) = 1 = Ze[ = Z1[
 = 8 − 1. (5) − 1. (1) = 2  2
Ÿ

Etape (4) :

Š.  = ‹ ¡ = š‹›
™
et

Maintenant, avec le vecteur y les inconnues du système Ux = y se déterminent comme suit :

1 1 1  5
Z0 1 1[ . š‹ › = Z 1[
0 0 1 ™ 2

Par un remplacement par arrière, le vecteur x est :

 4
¡ = š‹› = Z−1[
™ 2
On peut vérifier facilement en remplaçant les valeurs du vecteur x dans la première équation du système
d’origine :

+‹+™ =5
P4 − 1 + 2 = 5
5=5

• Algorithme de Factorisation : (A → LU) de Crout.

uii =1 i = 1,………..,n

“ = “ − ∑“V

b( b . b“ i = r,………...,n

r =1,.........,n

“ = “ − ∑“V

b( “b . b /““ j = r+1,…….,n

Application : Soit à résoudre le système d’équations linéaires suivant, avec la méthode de Crout:

.  =  … … (1)
 + 2‹ + 3™ = 14 ⎧
⎪ = dŠ … … (2)
P 2 + 5‹ + 2™ = 18 dŠ.  = 
3 + 2‹ + 5™ = 22 ⎨ d. ‹ =  … … (3)
⎪
⎩ Š.  = ‹ … … (4)

36
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

1ère méthode (Factorisation de la matrice A selon l’algorithme de Crout):

1 2 3
= Z2 5 2[
3 2 5

 =  −  . Š Ÿ  =  Ÿ  = 1 , avec : Š = Š = Š77 = 1
r=1

 =  −  . Š Ÿ  =  Ÿ  = 2

7 = 7 − 7 . Š Ÿ 7 = 7 Ÿ 7 = 3

Š = ( −  . Š )/ Ÿ Š =  / Ÿ Š = 2/1 = 2
----------------------------------------------------------------------------------

Š7 = (7 −  . Š7 )/ Ÿ Š7 = 7 / Ÿ Š7 = 3/1 = 3

 =  −  . Š Ÿ  = 5 − 2.2 Ÿ  = 1

r=2

7 = 7 − 7 . Š Ÿ 7 = 2 − 3.2 Ÿ 7 = −4

Š7 = (7 −  . Š7 )/ Ÿ Š7 = (2 − 2.3)/1

----------------------------------------------------------------------------------

Š7 = −4/1 = −4

77 = 77 − 7 . Š7 − 7 . Š7 Ÿ 77 = 5 − (3.3) − (−4. −4)
r=3

7 = −20
323

Avec la factorisation, notre matrice A de départ est transformée en deux matrices comme suit :

1 0 0 1 2 3
d = Z2 1 0 [ , Š = Z0 1 −4[
3 −4 −20 0 0 1

2ème méthode :

Etape (1) :

1 2 3  14
Z 2 5 2 [ . š‹› = Z18[
3 2 5 ™ 22

A . x = B

Etape (2) :

 0 0 1 ž ℎ
= dŠ d=Z œ 0[ Š = Z0 1 )[
f N  0 0 1
,

37
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

 0 0 1 ž ℎ
=Z œ 0[ . Z0 1 ) [
f N  0 0 1

 ž ℎ 1 2 3
=1 ž+œ ℎ + œ) 3 = Z 2 5 2 [
f fž + N fℎ + )N +  3 2 5

=1 ž = 2 (ž = 2) ℎ = 3 (ℎ = 3)
=1 =2 ž + œ = 5 (œ = 1) ℎ + œ) = 2 () = −4) 3
f=3 fž + N = 2 (N = −4) fℎ + )N +  = 5 ( = −20)

Avec la factorisation, notre matrice A de départ est transformée en deux matrices comme suit :

1 0 0 1 2 3
d = Z2 1 0 [ , Š = Z0 1 −4[
3 −4 −20 0 0 1

Etape (3) :

Déterminant d’abord y par la résolution de :

/
d. ‹ =  = Ze [

et

1 0 0 / 14
Z2 1 0 [ . Ze [ = Z 18[
3 −4 −20  22

On trouve :

/ = 14 / 14
6e = 18 − 2. (14) = −10 = Ze[ = Z−10[
 = ¢22 − 3. (14) − (−4). (−10)£/−20 = 2  3
Ÿ

Etape (4) :

Š.  = ‹ ¡ = š‹›
™
et

Maintenant, avec le vecteur y les inconnues du système Ux = y se déterminent comme suit :

1 2 3  14
Z0 1 −4[ . š‹› = Z−10[
0 0 1 ™ 3

38
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Par un remplacement par arrière, le vecteur x est :

 1
¡ = š‹› = Z2[
™ 3
On peut vérifier facilement en remplaçant les valeurs du vecteur x dans la première équation du système
d’origine :

 + 2‹ + 3™ = 14
P 1 + 2. (2) + 3. (3) = 14
14 = 14

a) Algorithme de Factorisation (A → LDU)

matrice U telle que LU=A. L’algorithme de Crout génère une matrice Š 5 unitaire et une matrice d5 telle
L’algorithme de Doolittle génère une matrice L unitaire (à éléments diagonaux égaux à l’unité) et une

que d5 Š 5 = .

En fait, on remarque que  = 5 , ∀i = 1,n.

En notant dii = uii on a : A = LU = (LD) (D-1U) = d5 Š 5 = dwŠ 5 = d5 wV Š

Ces algorithmes sont donc des cas particuliers des décompositions de A en :

A = LDU

Où L est une matrice triangulaire inférieure, U une matrice triangulaire supérieure et D une matrice
diagonale. La méthode de Cholevski sans racine carré est également une méthode de ce type.

Après décomposition de A → LDU, la résolution de Ax = b peut se décomposer en la résolution

successive des trois systèmes :

d‹ =
Pw™ = ‹
Š = ™

Ces systèmes linéaires sont respectivement résolus par les formules vues au début du chapitre (descente
triangulaire, remontée triangulaire, système à matrice diagonale). La résolution de Dz = y est inutile dans le
cas (D = I) des algorithmes de Crout et Doolittle.

III.5.2 Principe de la méthode de Cholevsky

Si A est une matrice symétrique définie positive, alors elle peut être décomposée en :

= dd¤

Où L est une matrice réelle triangulaire inférieure.

39
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

a) Décomposition de la matrice A

Développons l’équation : = dd¤

 = ' b . b

b(

b = 0 .) e > )
Où, par définition : Ž
b = 0 .) e > /

On a donc :

 = ∑b( b . b ,

‘’(,)
avec : i = 1,n et j = 1,n

Pour la partie triangulaire supérieure de A, la rème ligne s’écrit :

“ = ∑“b( “b . b , avec : j = r,n

Soit :

“ = ∑“V
b( “b . b + ““ . “ , avec : j = r,n

Ce qui donne l’algorithme suivant :

• Algorithme de Factorisation : (A → dd¤ ) de Cholevsky.

““ = ““ − ∑“V

b( “b 
 /

r =1,....,n

“ = “ − ∑“V

b( “b . b /““ j = r+1,….,n

Exemples:

4 1
=¥ ¦ ,
1 9
a)

r=1  = √ = 2 ,  =  / = 1/2

r=2  = ¨ − 


= ¨9 − 1/4

2 0
d’où : d = © ª
1/2 ¨9 − 1/4

On peut vérifier facilement que : = dd¤

40
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

1 1 3
= Z1 5 5[
3 5 19
b)

 = √ = 1 ,  =  / = 1/1 = 1 , 7 = 7 / = 3/1 = 3
 = ¨ −  = √5 − 1 = √4 = 2 , 7 = 7 −  . 7 / = 5 − 1.3/2 = 2/2 = 1
r=1


77 = ¨77 − 7 − 7 = √19 − 3 − 1 = √9 = 3

r=2
 
r=3

1 0 0
d’où : d = Z1 2 0[
3 1 3

On peut vérifier facilement que : = dd¤

2ème méthode :

1 1 3
= Z1 5 5[
3 5 19

= d. d¤ , avec A : une matrice symétrique

 0 0  f
d=Z œ 0[ d¤ = Z0 œ N[
f N  0 0 
,

 0 0  f   f
= d. d = Z
¤ œ 0[ . Z0 œ N [ = 1 
+ œ f + œN 3
f N  0 0  f f + œN f + N  +  

  f 1 1 3
= 1 
+ œ f + œN 3 = Z1 5 5[
f f f + œN f + N  +   3 5 19

=1  = 1 ( = 1) f = 3 (f = 3)
=1  
+ œ = 5 (œ = 2)

f + œN = 5 (N = 1) 3
f f + œN f + N  +   = 19 ( = 3)


 0 0 1 0 0
d’où : d = Z œ 0[ = Z1 2 0[
f N  3 1 3

On peut voir clairement que c’est le même résultat que la méthode précédente.

b) Résolution des systèmes linéaires : Ax = b par la méthode de Cholevsky

Après décomposition de A, la résolution de dd¤  = s’écrit :

d‹ =
Ž ¤
d=‹

41
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Et la résolution peut se faire de la même manière pour résoudre un système triangulaire (inférieur et
supérieur).

Remarques :

1. Si le terme (““ − ∑“V

b( “b ) < 0 , alors la matrice n’est pas définie positive. Ce test constitue un


moyen simple pour vérifier qu’une matrice A est bien définie positive.
2. On peut résoudre un système à matrice symétrique par une des méthodes précédentes (Gauss et Gauss-
Jordan). Toutefois, celles-ci ne tenant pas compte de la propriété de symétrie, elles impliqueront des
calculs redondants et donnant généralement des résultats moins précis. Cependant, on peut modifier la
méthode de Gauss, par exemple, afin qu’elle tienne compte de la symétrie et elle peut s’appliquer
même si A n’est pas définie positive.

c) Résolution d’un système [A,b] par factorisation de Cholevsky (technique de la matrice augmentée).

Après application de l’algorithme de décomposition, on a :

dd¤ =
Ž
d‹ =

C’est donc la même transformation dV qui fait passer de A en d¤ et de b en y. On peut donc simplifier
l’algorithme en formant la matrice augmentée [A,b] et en appliquant la décomposition de Cholevsky, la
matrice [A,b] devient :

dV z , ] = zd¤ , ‹]

Ceci fait, il suffit de résoudre le système :

d¤  = ‹

Avec la technique de la matrice augmentée les transformations sont faites sur place. La factorisation peut se
noter :

z , ] o zd¤ , ‹]

La résolution du système triangulaire résultant :

zd¤ , ‹] o zd¤ , ]

Permettra d’obtenir la solution x qui prendra la place du vecteur y qui lui-même a pris la place de b.

La résolution d’un système linéaire à l’aide de la factorisation = dd¤ de Cholevsky s’écrit donc :

• Résolution du système [A, b] par l’algorithme de Cholevsky.

42
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

1. Décomposition : dV z , ]
““ = ““ − ∑“V
b( b“ 
 /

“ = “ − ∑“V

b( b“ b /““
r = 1,…,n
j = r+1,…,n+1

2. Résolution de d¤  = ‹
, 0 = , 0 − ∑(0  . , 0 / i = n, n-1,….,1

Application : Soit à résoudre le système d’équations linéaires suivant, avec la méthode de Cholevsky:

.  =  … … (1)
+‹+™ =3 ⎧
= d. d¤ … … (2)
P − ‹ + ™ = 4
+‹−™ =5 ⎨ d. ‹ =  … … (3)
⎩d¤ .  = ‹ … … (4)

Etape (1) :

1 1 1  3
Z 1 − 1 1 [ . š‹› = Z4[
1 1 −1 ™ 5

A . x = B

Etape (2) :

 0 0  f
= d. d¤
d=Z œ 0[ d = Z0
¤ œ N[
f N  0 0 
,

 0 0  f
=Z œ 0[ . Z0 œ N [
f N  0 0 

  f 1 1 1
= 1 
+ œ f + œN 3 = Z 1 − 1 1 [
f f f + œN f + N  +   1 1 −1

=1  = 1 ( = 1) f = 1 (f = 1)
=1  
+ œ  = −1 (œ = )√2) f + œN = 1 (N = 0)
3
f f + œN f + N +  = )√2 ( = )√2)
  

Avec la factorisation, notre matrice A de départ est transformée en deux matrices comme suit :

1 0 0 1 1 1
d = Z1 )√2 0 [ , d = Z0 )√2
¤ 0 [
1 0 )√2 0 0 )√2

43
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Etape (3) :

Déterminant d’abord y par la résolution de :

/
d. ‹ =  = Ze [

et

1 0 0 / 3
Z1 )√2 0 [ . Z e [ = Z 4[
1 0 )√2  5

On trouve :

/=3 3
⎧ /
3 + )√2. e = 4 Ÿ e = 1/)√2 = () 2). () 2) = −)√2 /2
 )√2
= Ze[ = s −)  t
√

⎨
√ √

⎩  = )√2 = −)√2 −)√2

Ÿ

Etape (4) :

d . = ‹
¤
¡ = š‹›
™
et

Maintenant, avec le vecteur y les inconnues du système Lt.x = y se déterminent comme suit :

1 1 1  3
Z0 )√2 0 [ . š‹› = s −)  t
√

0 0 )√2 ™ −)√2

Par un remplacement par arrière, le vecteur x est :

 9/2
¡ = š‹› = Z−1/2[
™ −1
On peut vérifier facilement en remplaçant les valeurs du vecteur x dans la première équation du système
d’origine :

+‹+™=3
9 1

− −1=3
2 2
3=3

III.5.3 Résolution d’un système à matrice tridiagonale (Méthode de Thomas/TDMA Method).

Quand un système linéaire :

=‹

44
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

particulière de A. Elle est applicable aux matrices diagonalement dominantes c.à.d. |  | ≥ |œ | + | |,
est à matrice tridiagonale, alors l’algorithme de Gauss peut être simplifié, en tenant compte de la structure

(i=1,…,n). Ce cas particulier est important car on obtient très souvent un tel système lors de la résolution
numérique de certaines équations différentielles partielles, ou en opérant une interpolation à l’aide de fonctions
Splines.

Adoptons pour les éléments non nuls de A, la notation suivante :

œ
⎡  ⎤
⎢   œ 0 ⎥
⎢ 7 7 œ7 ⎥
=⎢ ………..… ⎥
⎢  œ ⎥
(*)
⎢  ⎥
⎢ 0 ………. … … . œ V ⎥
⎣  ⎦

Ce qui signifie que l’on peut ne mémoriser que les (3n-2) termes non nuls de A dans trois vecteurs a, b, c ayant
respectivement (n-1), (n) et (n-1) éléments, plutôt que de mémoriser (et traiter) tous les termes de A dont (n2-
3n+2) sont nuls. L’algorithme de Thomas est le suivant :

Normalisons la première ligne du système :

 ‹
œ ⎡ ⎤ ⎡ ⎤
⎡  ⎤ ⎢  ⎥ ⎢‹ ⎥


⎢ 
0 ⎥ ⎢ ⎥ ⎢ ⎥
⎢   œ ⎥ ⎢ 7 ⎥ ⎢‹7 ⎥
⎢ 7 7 œ7 ⎥ ⎢ ⎥ ⎢ ⎥
⎢ ………..… ⎥ .. = ⎢..⎥
⎢  ⎥ ⎢ ⎥ ⎢ ⎥
⎢  œ ⎥ ⎢  ⎥ ⎢ ‹ ⎥
⎢ 0 ………. ……. œ V ⎥ ⎢ . . ⎥ ⎢ ⎥
..
⎣  ⎦ ⎢ ⎥ ⎢ ⎥
⎣ ⎦ ⎣‹ ⎦

Puis, réduisons la seconde en lui soustrayant  fois la première :

‹
œ  ⎡ ⎤
⎡  ⎤ ⎡ ⎤ ⎢
‹ ⎥


⎢ 
œ ⎥ ⎢  ⎥ ⎢‹ −  .  ⎥
⎢0 ? −  . @ œ 0 ⎥ ⎢ 7 ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ‹7

 ⎥
⎢ 7 œ7 ⎥⎢ .. ⎥ = ⎢
.. ⎥
⎢ ⎥ ⎢ ⎥
7
………..… ⎢ ⎥
⎢ ⎥ ⎢  ⎥
⎢   œ ⎥ ⎢ ⎥ ⎢ ‹  ⎥
⎢ 0 ………. ……. œ V ⎥ ⎢ . . ⎥ ⎢ .. ⎥
⎣  ⎦ ⎣ ⎦ ⎢ ⎥
⎣ ‹ ⎦

45
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

Normalisons la seconde ligne :

œ ‹
⎡  ⎤  ⎡ ⎤
⎡ ⎤ ⎢ ⎥
⎢ ⎥ 
 
⎢ ⎥ ⎢ ‹ −  . ‹ / ⎥
⎢ ⎥
 
œ 0 ⎢ ⎥ ⎢  −  . œ /
⎢0 1 ­ ⎥
œ ® ⎥ 7
⎢ ⎥ ⎢ ‹7 ⎥
⎢ −  . ⎥
⎢ .. ⎥ = ⎢ ⎥

⎢

⎥ 
7 œ7 ⎢ ⎥ ⎢ .. ⎥
⎢ 7
⎥ 
………..… ⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ‹ ⎥
⎢  œ ⎥ ⎢ .. ⎥ ⎢

.. ⎥
⎢ 0 ………. ……. œ V ⎥
⎣ ⎦ ⎢ ⎥
⎣  ⎦ ⎣ ‹ ⎦

Puis, réduisons la troisième :

œ
⎡  ⎤ 
⎢ 
⎥ ⎡ ⎤
⎢ œ ⎥ ⎢  ⎥
⎢0 1 ­ œ ® ⎥ ⎢ ⎥
⎢  −  . 0 ⎥ ⎢ 7 ⎥
⎢ ⎥ ⎢ € ⎥


⎢ œ ⎥⎢ ⎥
⎢0 0 1 ­ 7 − 7 œ ® œ7 ⎥ ⎢ ..⎥
⎢  −  . ⎥ ⎢ ⎥

⎥ ⎢ ⎥

⎢  œ
⎢ …€… … . . … ……
€
. … €… ⎥ ⎢ ..⎥
⎢   œ ⎥ ⎢ ⎥
⎢ 0 ………. ……. ⎥ ⎣ ⎦
⎣  ⎦
‹
⎡ ⎤
⎢ 
⎥
⎢ ‹ −   . ‹ /  ⎥
⎢  −  . œ /  ⎥
⎢ ‹ −  . ‹ / ⎥
= ⎢‹7 − 7 . −  . œ / ⎥
⎢    ⎥
‹
⎢ € ⎥
⎢ .. ⎥
⎢ ‹  ⎥
⎢ .. ⎥
⎣ ‹ ⎦

En continuant ainsi pour les n lignes du système, on obtient un système à matrice bidiagonale unitaire,
(éléments diagonaux égaux à 1)

46
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

1 ¯  °
⎡ ⎤ ⎡ ⎤ ⎡  ⎤
1 ¯ 0 °
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
7
⎢ 1 ¯7 ⎥ ⎢ ⎥ ⎢ °7 ⎥
⎢ …… ⎥⎢ .. ⎥ = ⎢ .. ⎥
⎢ 1 ¯ ⎥ ⎢  ⎥ ⎢ ° ⎥
⎢ 0 ……. ⎥ ⎢ .. ⎥ ⎢ .. ⎥
⎣ 1 ⎦ ⎣ ⎦ ⎣° ⎦

L’algorithme de Thomas, dans un format approprié pour la programmation informatique est résumé comme
suit :

• Algorithme de Thomas pour la résolution d’un système linéaire Ax = y à matrice tridiagonale.

Soit la notation A= [a, b, c] définie en (*)

• Triangularisation (matrice tridiagonale devient bidiagonale).

œ
¯ =


¯ = (² V9 i.³
±
i i i–: )
i = 2,…,n-1

‹
° =


° = (²i V9i .³ i–:)

´ V9 .µ
i i i–:
i = 2,…,n

• Résolution du système à matrice bidiagonale (cas particulier d’un

système à matrice triangulaire supérieure).

 =°

 = ° − ¯ . 0 i = n-1, n-2,…, 1

Le code Matlab qui permet de résoudre un système tridiagonal selon l’algorithme de Thomas :

function x = thomas(a,b,c,y,N)

% La fonction thomas:
% Inverse un système tridiagonal A.x=y dont les diagonales inférieure,
% principale et supérieure sont respectivement donnés par les vecteurs
% a, b et c.
% y est le membre droit du système, et N est la taille du système.
% Le résultat est placé dans le vecteur x.
% >> a=[0 2 1 1];b=[2 3 4 3];c=[1 1 2];y=[1 2 3 4];N=4;
% x = thomas(a,b,c,y,N)

if (b(1)==0)

47
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

fprintf(1,'Réorganiser les équations pour le solveur tridiagonal...')

pause
end

Gamma(1)=c(1)/b(1); Betta(1) = y(1)/b(1);

% Lancer la décomposition et la substitution directe

for j = 2:N-1
D = b(j)-a(j-1)*Gamma(j-1);
Gamma(j) = c(j)/D;

if (D==0)
fprintf(1,'Le solveur tridiagonal a échoué...')
pause
end
end

for j = 2:N
D = b(j)-a(j-1)*Gamma(j-1);
Betta(j) = (y(j)-a(j-1)*Betta(j-1))/D;

end

% Effectuer la substitution arrière

x(N)=Betta(N);
for j = 1:(N-1)
k = N-j ;
x(k) = Betta(k) - Gamma(k)*x(k+1) ;
end

Application : Soit à résoudre le système d’équations linéaires suivant (A.x = y), avec la méthode de Thomas:

2 +  = 1
 + 3 + 7 = 2


 + 47 + 2€ = 3
7 + 3€ = 4

Etape (1) :

2 1 0 0 1 1
12 3 1 03 . s2 t = s2t
0 1 4 2 3 3
0 0 1 3 4 4

A . x = y

48
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

En utilisant l’algorithme de Thomas on obtiendra un système à matrice bidiagonale unitaire, (éléments

diagonaux égaux à 1) à la place d’un système tridiagonal.

1 ‹1 ¯1 °
œ1 0 0 1 0 1
0 ⎡ 1⎤
‹2 ⎢° ⎥

  œ2 0 2 0 1 ¯2 0 2
s t
œ3 3. s t = s ‹3 t s t . s t = ⎢ 2 ⎥
0 7 0 0 1 ¯3 3 °
7
4 ‹4  ⎢ 3⎥
Ÿ
0 0 € € 0 0 0 1 4
⎣°4 ⎦
œ 1
¯ = ¯ =



2
‹ 1
° = ° =
Ÿ
 2

Etape (2) :

1 1
œ ⎧¯ = =
1
⎧¯ = ( −  . ¯ ) ⎪
⎪ (3 − 2. 2) 2
  V
1
⎨ ° = ‹ −  . °V ⎨ 2 − 2 .2 1
⎪° = =
Ÿ
⎩  (  −  . ¯V ) ⎪ 1
⎩ A3 − 2. 2B 2

Etape (3) :

2 4
œ7 ⎧¯7 = =
1
⎧¯7 = ( −  . ¯ ) ⎪ (4 − 1. 2) 7
7 7 7V
‹ −  . ° 1
⎨° = 7 7 7V ⎨ 3 − 1 .2 5
⎪°7 = =
Ÿ
⎩ 7
( 7 − 7 . ¯7V ) 1
⎩ (4 − 1. 2) 7

Etape (4) :

5
‹€ − € . °€V 4 − 1. 7 23/7 23
6°€ = °€ = = =
( € − € . ¯€V ) 4
A3 − 1. 7B 17/7 17
Ÿ

Etape (5) :

Par un remplacement par arrière, le vecteur x est :

 =° ;

 = ° − ¯ . 0 , avec : i = n-1, n-2,…, 1

⎧ € = ƒ
7
4 = °4  4/17
⎧ ⎪
⎪ 7 =  − € . 7 = − 

= °3 − ¯3 . 4 9/17
¡ = s t = s t
3 ƒ ƒ ƒ ƒ

⎨2 = °2 − ¯2 . 3 ⎨ = −  . A− ƒB = ƒ −1/17

   „
7
€
Ÿ Ÿ

⎩1 = °1 − ¯1 . 2 ⎪
⎪

23/17
⎩ = −  . AƒB = ƒ
  „ €


49
 UHBC, Département de mécanique 2ème LMD (S4) Cours d’analyse numérique

III.5.4 Conclusions sur les méthodes directes

Les méthodes directes sont surtout employées pour la résolution des petits (n < 20) et moyens systèmes
(n < 100) d’équations linéaires.

Comme le temps de calcul croit aussi avec le nombre d’opérations, la meilleure méthode directe au sens du
temps de calcul est celle qui demande le moins d’opérations (la moins complexe).

Pour cette raison, la méthode de Cramer n’est jamais employée pour résoudre des systèmes d’ordre n supérieur
à 5.

La méthode de Gauss (ou sa variante de Crout) est souvent préférée pour les systèmes à matrice A quelconque.
Encore que manuellement, celle de Jordan est plus pratique d’emploi.

Quand la matrice A est symétrique, l’algorithme de Gauss peut se simplifier en conséquence. Si A est
tridiagonale, l’algorithme de Gauss simplifié s’appelle algorithme de Thomas.

La méthode de Cholevsky est recommandée pour la résolution des systèmes linéaires à matrice symétrique
définie positive. De tels systèmes sont courants dans la description des phénomènes physiques, ou dans les
problèmes d’optimisation.

50