Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 PDF

, Chapitre 2 Méthodes directes de résolution des systèmes linéaires Pr HADDADOU Hamid HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Plan 1 Méthode directe de résoluti...

, Chapitre 2 Méthodes directes de résolution des systèmes linéaires Pr HADDADOU Hamid HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Plan 1 Méthode directe de résolution : Position du problème et idée de résolution 2 Méthode d’élimination de Gauss 3 Factorisation et méthode LU 4 Factorisation et méthode de Cholesky 5 Factorisation et méthode QR 6 Conditionnement du problème et stabilité des méthodes 7 Exemple motivant 8 Annexe 1 : La factorisation QR par Householder HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Méthode directe de résolution : Position du problème et idée de résolution Méthode directe de résolution : Position du problème et idée de résolution HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Méthode directe de résolution : Position du problème et idée de résolution Méthode directe de résolution Définition Considérons le système AX = b avec A ∈ Mnn (C) inversible et b ∈ Cn tq A = (a1 , · · · , an ) et b = t (b1 , · · · , bn ) où ak est le k me vecteur colonne de A. On dit qu’une méthode est directe si elle donne la solution exacte X en un nombre fini d’opérations élémentaires. Exemple La formule de Cramer est une méthode directe qui donne la solution exacte X = t (xk , · · · , xn ) avec det(a1 , · · · , abk · · · , an ) xk = , avec abk = b det A Position du problème Sachant que le calcul directe d’un déterminant d’une matrice pleine d’ordre n demande ' n(n!) opérération élémentaires. Si on dispose d’un ordinateur qui effectue 109 d’opérations par seconde, le calcul d’un n(n! déterminant s’effectue en un temps équivalent à 109 s. 50(5!) Par exemple si n = 50, le calcul s’effectue en un temps équivalnt à 109 s. 50(50!) C’est à dire 109 ×3600×24×365 années ' 4, 82 1048 années HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Méthode directe de résolution : Position du problème et idée de résolution Idée de réolution Cas de matrice triangulaire et l’idée de résolution Si A est triangulaire supérieure inversible, le système à résoudre se présente sous la forme suivante : a11 x1 + · · · · · · · · · · · · · · · + a1n xn = b1  .   .... .   a(n−1)(n−1) xn−1 + an(n−1) xn = bn−1    ann xn = bn On calcule un puis on remonte xn = abn    nn 1 
 xn−1 = a bn−1 − an(n−1) xn   (n−1)(n−1).. ⇒ Nombre d 0 opérations = n2    .  1  x = 1 a11 (b1 − (a12 x2 + · · · + a1n xn )). Idée de résolution Dans le cas où A n’est pas triangulaire supérieure ou inférieur, l’idéé est de transformer AX = b à un système équivalent dont la matrice est triangulaire et de tel sorte que la transformation nécessite un nombre d’opérations élémentaires négligeable devant n!. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss La méthode d’élimination de Gauss HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss La méthode d’élimination de Gauss AX = b ⇔ MAX = Mb La méthode Gauss est constituée de 3 étapes 1) Procédure d’élimination : déterminer M telle que MA soit triang. sup., 2) calculer Mb, 3) résolution de MAX = Mb. Remarque Dans la pratique, on peut faire les combinaisons linéaires entre les équations et ceci donne directement MA et Mb en même temps. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss Procédure d’élimination (Exemple) On explique cette méthode par un exemple. Considérons le système suivant :   2x1 + 4x2 − 4x3 = 1 3x1 + 6x2 + x3 = 0 −x1 + x2 + 3x3 = 0  Étape 1 : Notons A la matrice du système et posons A1 = A et b1 = b. Alors,  1  l1 ! 2 4 −4 1 2 4 −4 A|b= 3 6 1 0 ∼ l2 − 32 l1  0 0 7 − 32  −1 1 3 0 l3 + 12 l1 0 3 1 1 | {z }| 2 {z } A2 b2   l1 2 4 −4 2 1 ∼ l3  0 3 1 2  3 l2 0 0 7 −2 | {z }| {z } A3 b3 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss Procédure d’élimination (suite de l’exemple 1) Ainsi,   2x1 + 4x2 − 4x3 = 1, AX = b ⇔ 3x1 + 6x2 + x3 = 0, −x1 + x2 + 3x3 = 0,    2x1 + 4x2 − 4x3 = 1, ⇔ 3x2 + x3 = 21 , 7x3 = − 23 ,  3   x3 = − 14 , méthode de 5 ⇒ x2 = 21 , relontée  x1 = 17 42. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss Matrice de transposition et matrice d’élimination Remarque Matrice de transposition Échanger 2 lignes lj et lk dans une matrice équivaut à la multiplier à gauche par la matrice de transposition T (lj , lk ) obtenue en échangeant les lignes lj et lk dans la matrice identité In. Matrice d 0 élimination : Posons −l(k+1)k = le coeff. de la lig. du pivotage à l’étape k Colonne k ↓ 1 ··· ···   ..   .  .     1..   Ek =  .  ligne k .....   → ligne (k + 1) .. −l(k+1)k.    ..  . 1  −lnk 1 Exemple Résoudre par la méthode de Guasse AX = b et déduire P1 , P2 , E1 , E2 , où     1 1 1 1 A =  3 −3 4  et b =  0 . 2 1 3 0 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss Solution On a   1 1 1 1 A|b ∼  3 −3 4 0  ↔ P1 AX = P1 b 2 1 3 0   l1 1 1 1 1 ∼ l2 − 3l1  0 −6 1 −3  ↔ E1 P1 AX = E1 P1 b l3 − 2l1 0 −1 1 −2   l1 1 1 1 1 ∼ l2 − 3l1  0 -6 1 −3  ↔ P2 E1 P1 AX = P2 E1 P1 b l3 − 2l1 0 −1 1 −2   l1 1 1 1 1 ∼ l2  0 −6 1 −3  ↔ E2 P2 E1 P1 AX = E2 P2 E1 P1 b l3 − 16 l2 0 0 5 6 − 32 Ensuite on résout le denier système avec la méthode de remonté. On a P1 = P2 = I3 , et     1 0 0 1 0 0 E1 =  −3 1 0 , E2 =  0 1 0 . −2 0 1 0 − 16 1 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Méthode d’élimination de Gauss Procédure d’élimination (suite) En résumé d’une façon général pour une matrice inversible de taille n, la procédure d’élimination de Gauss se traduit par AX = b ⇔ En−1 Pn−1...E1 P1 Ax = En−1 Pn−1...E1 P1 b ⇔ MAX = MB | {z } M où les Pi sont les matrices de permutation, et Ei sont les matrices d’élimination. Corollaire. Si on pose U = MA, alors det(A) = (−1)p det(U), où p est le nombre total des permutations de lignes effectuées. Remarque Pour n assez grand, le nombre d’opérations dédié à la méthode Gauss est équivalent à 2 3 3 n HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Factorisation et méthode LU HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Méthode de la factorisation LU (de Gauss) Le but ici c’est d’écrire la matrice A sous la forme d’un produit d’une matrices L tringulaire inférieure avec Ljj = 1 et U une matrice tringulaire supérieure. Dans le cas où det A 6= 0, la méthode de Guass est applicable et la matrice. U = En−1 Pn−1...E1 P1 est triangulaire supérieure. Alors, | {z } M A = M −1 U. Supposons maintenant qu’on a jamais échangé des lignes (ie : Pi = In ). Alors, M −1 = (E1 )−1 · · · (En−1 )−1. Puisque toutes les Ek sont tringulaires infériueres alors M −1 l’est aussi (car produit de matrices tringulaires inférieurs). On pose ainsi L = M −1 et on obtient A = LU. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Détermination de L et U En passant par la procédure d’élimination de Gauss : On obtient U et on montre que 1 ··· ··· 1 ···     ...  ..       .  ....          1.  ⇒ E −1 =    1.  Ek =  . k  ....  .... ..... .. −l(k+1)k .. l(k+1)k          ..  ..  . 1  . 1  −lnk 1 lnk 1 Puis  ···  1 0 0 ......   l ...   L =  21 ......   ... 0  ln1 ··· ln(n−1 1 Par identification : On peut déterminer L et U par identification et en utilisant les formes de L (triangulaire inférieure) et de U (triangulaire inférieure). HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Exemple (Factorisation LU) ! ! ! 1 1 1 l1 1 1 1 l1 1 1 1 A= 3 −3 4 ∼ l2 − 3l1 0 −6 1 ∼ l2 0 −6 1 2 1 3 l3 − 2l1 0 −1 1 l3 − 61 l2 0 0 5 6 La méthode de Gauss est applicable avec P1 = P2 = I3 , et     1 0 0 1 0 0 E1 =  −3 1 0  , E2 =  0 1 0 . −2 0 1 0 − 16 1     1 1 1 1 0 0 Ce qui implique U =  0 −6 1  et L =  3 1 0  5 1 0 0 6 2 6 1 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Condition nécessaire et suffisante pour la factorisation LU Théorème a11 · · · a1k   ...  Soit A = (aij ) et ∀k = 1,... , n, ∆k = ... . Alors, la matrice ak1 · · · akk A une unique factorisation LU de Gauss si et seulement si ∀k = 1,... , n, det(∆k ) 6= 0. Il y a deux cas particuliers de matrices admettant la factorisation LU, plus précisement A est symétrique définie positive ⇒ det(∆k ) > 0, ∀k = 1,... , n A est à diagonale strictement dominante ⇒ det(∆k ) 6= 0, ∀k = 1,... , n Remarque Une matrice carrée A est P dite à diagonale strictement dominante par lignes si ∀i, |aii | > |aij | , j, j6=i ou P par colonnes si ∀j, |ajj | > |aij | , i, i6=j HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Factorisation LU de crout et factorisation PLU Remarque (Factorisation LU de crout) Il y aune autre version de la factorisation LU dûe à crout, où on s’intéresse à écrire A sous la forme d’un produit d’une matrices L tringulaire inférieure et U une matrice tringulaire supérieure avec Ujj = 1. Théorème (Factorisation PLU) Soit A inversible. Alors, il existe une matrice de permutation P, une matrice L = (lij ) triangulaire inférieure avec des 1 sur la diagonale et une matrice U triangulaire supérieure telles que PA = LU. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Résolution d’un système à partir de la factorisation LU ou PLU Si l’on a à résoudre plusieurs systèmes avec la même matrice, la décomposition LU se révèle trés utile. On résout donc par Gauss une fois un système et on déduit la décomposition LU ou PLU. - Si A = LU , on résout les autres systèmes comme suite : 1) On trouve la solution de LY = b par la méthode de descente 2) Après on trouve X en utilisant la méthode de remontéé UX = Y car L et U sont triangulaires respectivement inférieur et supérieure. - Si PA = LU, il vient PAX = Pb. On résout dles autres systèmes comme suite : 1) LY = Pb 2) UX = Y Or L et U sont triangulaires, il suffit d’utiliser respectivement la méthode de remontée et la méthode de descente. Remarque Pour une matrice A, pour déterminer les matrices L, U et P, il suffit de taper : >> [L, U, P] = lu(A) HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode LU Exemple     1 1 1 1 Soient A =  1 3 5  et b =  0  2 8 15 0 On a det(∆1 ) = 1 6= 0, det(∆2 ) = 2 6= 0 et det(∆3 ) = 1 6= 0, donc A admet une factorisation LU.       1 1 1 l1 1 1 1 l1 1 1 1 A= 1 3 5  ∼ l2 − l1  0 2 4  ∼ l2  0 2 4  2 8 15 l3 − 2l1 0 6 13 l3 − 3l2 0 0 1     1 0 0 1 0 0 On en déduit E1 =  −1 1 0  , E2 =  0 1 0  et −2 0 1 0 −3 1     1 0 0 1 1 1 L =  1 1 0 , U =  0 2 4 . 2 3 1 0 0 1 Pour résoudre AX = b → on résout d’abord LY = b, on trouve Y = t (1, −1, 1). → Ensuite, on résout d’abord UX = Y , on trouve X = t ( 52 , − 52 , 1) HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky Factorisation et méthode de Cholesky HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky La factorisation et la méthode de Cholesky On considère ici des matrices réelles symétriques et définies positives. Définition On dit qu’une matrice symétrique A ∈ Mn (R) est définie positive si t XAX > 0, ∀X ∈ Rn \ {0Rn }. Exemple   2 1 La matrice A = est symétrique et définie positive 1 2 Théorème Soit A ∈ Mn (R). Alors, on a les équivalences suivantes : A est définie positive ⇔ SP(A) ⊂ R+∗ ⇔ det ∆k > 0 (critère de Sylvester). Remarque Remarquons que ces matrices admettent la factorisation LU, mais on va donner une ”meilleure” factorisation. Si une matrice A est symétrique et définie positive, alors aii ∈ R∗+. Pour le constater il suffit de prendre dans la définition X = ei (le ième vecteur de la base canonique). HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky Théorème (Factorisation de Cholesky) Soit A ∈ Mn (R) symétrique définie positive. Alors, il existe une matrice réelle triangulaire inférieure C = (cij ) (pas unique) telle que A = C t C qui est appellée factorisation de Cholesky. De plus, il existe une seule décomposition de Cholesky de A avec cii > 0, ∀i = 1,... , n. Remarque Inversement si la matrice inversible A admet une décomposition de Cholesky, alors elle est symétrique définie positive. Interêt Résolution des systèmes de la forme Au = b, Remplacer A par C. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky Calcul de C pour une matrice symétrique et définie positive En utilisant la factorisation LU Si on a déja calculé L et U. On prend (de la démonstration du dernier théorème !) C = LΛ avec (Λ est bien définie car uii > 0)  √ u11 0 ··· 0  ....   0..  Λ= .  .....  0  √ 0 ··· 0 unn. En utlisant l’identification On peut calculer C d’une manière pratique sans utiliser la factorisation LU. On pose C = (cij ) (C est triangulaire supérieure) et on suppose que A = C t C. n Alors, aij = (C t C )ij = P cik cjk k=1 Remarque Sous Matlab, pour une matrice A donnée la commande chol(A) donne t C >> chol(A) HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky Exemple de factorisation Exemple   1 1 1 La matrice A =  1 5 5  On a A = LU avec (exo) 1 5 14     1 0 0 1 1 1 L =  1 1 0 , U =  0 4 4 . 1 1 1 0 0 9 Par ailleurs, puisque A est symétrique et définie positive (car det(∆1 ) = 1 > 0, det(∆2 ) = 4 > 0, det(∆3 ) = 36 > 0), alors elle admet une factorisation de Cholesky C t C avec C = LΛ où   1 0 0 Λ =  0 2 0 , 0 0 3 ainsi   1 0 0 C = 1 2 0 . 1 2 3 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode de Cholesky Méthode de résolution de Cholesky etn nombre d’opération Méthode de Cholesky de résolution du système Soit A une matrice symétrique définie positive et Au = b avec b donné. Alors, on procède de la manière suivante : factorisation de Cholesky on résout CY = b et ensuite on résout t Cu = Y Le nombres d’opérations Pour n assez grand le nombre d’opérations élémentaires nécessaires pour résoudre un système par la méthode de Cholesky dans le cas où on utilise la deuxième méthode pour calculer C, est équivalent à 31 n3 opérations. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR Factorisation et méthode QR HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR Factorisation QR : But, Théorie et Intérêt But Le but est d’écrire toute matrice A sous la forme QR avec Q matrice orthogonale et R matrice triangulaire supérieure. Définition On dit qu’une matrice Q ∈ Mn (R) est orthogonale si Q −1 = t Q. Théorème Soit A ∈ Mn (R). On peut toujours factoriser A sous la forme QR. De plus on peut toujours s’arranger pour que tous les coefficients diagonaux de R soient tous positifs ou nul. Intérêt Résolution de systèmes linéaires. Sachant la décomposition QR d’une matrice A carrée ou non , la résolution d’un système AX = b revient à calculer t Qb et résoudre ensuite le système RX =t Qb. Si A est carrée, |det A| = |r11 ×... × rnn | car det(Q) = ±1. Calcul des Valeurs propres. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR La factorisation QR en utilisant Cholesky !(Première idée) Dans le cas des matrices carrées inversible En utilisant la factorisation de Cholesky ! Soit A ∈ Mn (R) inversible. 1) On pose B = t A A. Montrons que B admet une factorisation de Cholesky de la forme C t C avec C une matrice triangulaire inférieure. 2) Montrons que la matrice A(t C )−1 est orthogonale. 3) Déduisons que A admet une factorisation QR en posant Q = A(t C )−1 et R =t C. Remarque Ce procédé est coûteux en opérations donc pas pratique, néanmoins il présente une preuve que tout matrice inversible peut s’écrire sous la forme QR. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR Factorisation QR par le procédé d’orthonormalisation de Gram-Schmidt Ici on considère toujours le cas des matrices carrées inversible Théorème Soit A ∈ Mn (R) inversible. Alors, il existe une matrice orthogonale Q et une matrice triangulaire supérieure R à diagonale positive telles que A = QR. De plus cette décomposition est unique et peut être obtenue par le procédé d’orthonormalisation de Gram-Schmidt de la base formée des vecteurs colonnes de A, c’est à dire : Si A = (a1 ,..., an )  r11 = ka1 k2 , 1 q1 = r11 a1 , et    rim = hqi , am i , m = 2,... , n et i = 1,... , m − 1,  m−1 P am = am − e rim qi ,   i=1 1 rmm = keam k2 , qm = rmm am ,  e où h·, ·i désigne le produit scalaire. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR a1 a2 ··· an q1 q2 ··· qn......    ...... r11 r12 r1n  .. ···... ···..     ...  ...  0 r22.  ...  ... . r2n  .. ···.   =.. ···.    ......  ...  ... .... .......... ···... ···.         . .......    .......0 ··· 0 rnn.. ···... ···. | {z } | {z } | {z } R A Q r =< q1 , a2 >,    r11 = ka1 k2 ,  12  1 ae2 = a2 − r12 q1 , a1 −→ a2 −→ 1  q1 = a1 ,  r22 = kae2 k2 et q2 = ae2 , r11  r22 r =< q1 , a3 >, r23 =< q2 , a3 >,   13  ae3 = a3 − r13 q1 − r23 q2 , a3 −→  r = kae k et q = 1 ae ,  33 3 2 3 3 r33 r =< q1 , an >, · · · , rn−1n =< qn−1 , an >,   1n  aen = an − r1n q1 − · · · − rn−1n qn an −→ 1  rnn = kaen k2 et qn =  aen. rnn HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Factorisation et méthode QR Exemple et Défaut Exemple   1 1 Décomposer la matrice A = sous la forme QR 1 2 1 ! √ ! − √12 2 √32   1 1 √ = 2. 1 √1 1 1 2 √ 2 2 0 √ 2 Remarque La factorisation QR en utilisant Ckolesky est trés coûteuse en nombre d’opérations relativement aux autres méthodes Le procédée d’orthonormalisation de Gram-Schmidt de factorisation sous la forme QR est peu stable. il y a d’autres méthodes pour factoriser sous forme QR. Par exemple,une qui utilise ce qu’on appelle les rotations de Givensn, ou une autre dite de Householder qui est plus stable que le procédée dorthonormalisation de Gram-Schmidt. De plus, plus général, vu qu’elle traite tout type de matrice (carrée ou non carrée). Cette dernière méthode est décrite dans l’annexe de ce chapitre. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025 Conditionnement du problème et stabilité des méthodes Conditionnement du problème et stabilité des méthodes HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Choix du pivot pour éviter l’instabilité Pour appliquer la méthode, il suffit de prendre des éléments non nuls comme pivots. Pour illustrer l’influence des choix sur la qualité du résultat on consid‘ère l’exemple suivant : Problème : Erreurs d’arrondis lors de calcul sur une machine On choisit une mantisse à 3 chiffres décimaux et la base 10. Prenons l’exemple du système suivant : 10−4 u1 + u2 = 1  la solution exacte est ' (1.00010.., 0.99990) (1) u1 + u2 = 2 → On prend 10−4 comme pivot, le système (2) est équivalent à 10−4 u1 + u2 = 1 10−4 u1 + u2 = 1    u1 = 0, 4 4 pour ⇒ (1 − 10 ) u2 = 2 − 10 ⇒ −9990 u2 = −9990 u2 = 1 la machine Solution mauvaise. → On prend 1 comme pivot, le système (2) est équivalent à    u1 + u2 = 2 u1 + u2 = 2 u1 = 1, −4 −4 pour ⇒. (1 − 10 ) u2 = 1 − 2 × 10 ⇒ 0, 999 u2 = 0, 999 u2 = 1 la machine Solution acceptable. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Stratégie des choix En général les pivots petits en valeurs absolues sont mauvais. On peut alors développer des strategies dans le choix des pivots.. (k) Les principales Si on note Ak = (aik ) la matrice obtenue à l’étape k en appliquant l’élimination de Guass à A. 1) Stratégie du pivot partiel : On prend à l’étape k comme pivot (k) (k) aik = max alk. k≤l≤n 1) Stratégie du pivot total : On prend à l’étape k comme pivot (k) (k) aik = max alh. k≤l,h≤n HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Conditionnement Prenons l’exemple classique suivant (dû à R. S.Wilson) :       10 7 8 7 32 32, 1  7 5 6 5   23  0  22, 9  A=   , b=   et b =   . 8 6 10 9  33  33, 1  7 5 9 10 31 30, 9 Les solutions des systèmes AX = b et AX 0 = b 0 sont respectivement     1 9, 2  1  0  −12, 6  X =  1  et X =  4, 5 .    1 −1, 1 On remarque que     0, 1 −8, 2 0  −0, 1  0  13, −  b −b =  0, 1  et X − X =  0 − 3, 5     −0.1 2, 1 Donc, une petite perturbation sur le membre de droite a conduit un grand changement dans la solution. Chose qui semble être étrange à première vu. Pour expliquer ce phénomene on est conduit à parler de la notion de conditionnement. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Définition et propriétés Définition Soit A ∈ Mnn (R). Alors, pour une norme matricielle induite donnée k·k, le conditionnement de la matrice A est défini par cond(A) = kAk A−1. Proposition Pour toute matrice A ∈ Mnn (R) inversible, On a (1) cond(A) = cond(A−1 ), (2) cond(αA) = cond(A), ∀α ∈ R∗ , (3) cond(A) ≥ 1, λmax (4) cond2 (A) = (λmax et λmin sont la plus grande et la plus petite valeurs λmin singulières de A). (5) De plus, pour le cas particulier de la normé k·k2 , pour une matrice A est orthogonale, on a cond2 (A) = 1. Remarque Pour une matrice inversible A ∈ Mnn (R). Les valeurs sigulières d’une matrice A sont les racines carrées des valeurs propres de la matrice t A A. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Résultats importants Les deux résultats suivants expliquent le phénomène illustré par l’exemple ci-dessus : Théorème Soient A ∈ Mnn (R) et b, b 0 ∈ Rn avec b 6= 0.. Alors si X et X 0 sont solutions de AX = b et AX 0 = b 0 , on a kX 0 − X k kb 0 − bk ≤ cond(A). kX k kbk Théorème Soient A, A0 ∈ Mnn (R) et b ∈ Rn avec A 6= 0 et b 0 6= 0. Alors si X et X 0 sont solutions de AX = b et A0 X 0 = b, on a kX 0 − X k kA0 − Ak 0 ≤ cond(A). kX k kAk HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Notions et discussions Le conditionnement donne une idée de la sensibilité de la solution d’un système aux erreurs sur les données, à savoir sur la matrice A ainsi que sur le membre à droite b. Si le conditionnement est très élevé par rapport à 1 cond(A) >> 1, alors le système est dit mal conditionné, le conditionnement est très proche de 1, alors le système est dit bien conditionné. cond(A) = 1, alors le système est parfaitement bien conditionné donc si cond(A) = 1, c’est le cas des matrices orthogonales. Il est à signalé enfin que le conditionnement d’un système d’équations linéaires dépend seulement de la matrice du système. L’analyse de sensibilité de la solution aux perturbations sur les données demande le calcul de cond(A), mais son calcul demande le calcul de l’inverse de A qui est trop coûteux en complexité (de l’ordre de O(n3 )). Dans la pratique, on préfere calculer une approximation du cond(A), pour cela ils existent des algorithmes pour calculer l’inverse du cond(A) avec une complexité d’ordre O(n2 ). L’algorithme LAPACK est un exemple de ceux ci et la cammande rcond de matlab l’utilise. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Conditionnement du problème et stabilité des méthodes Exemple donnant une idée comment améliorer le conditionnement Soit la matrice suivante :   100 0 0 A= 0 10 0  0 0 1 On a Cond∞ (A) = 100. Mais, si on multiplie A à Gauche par la matrice D où  1  100 0 0 1 D= 0 10 0  0 0 1 on obtient Cond∞ (DA) = 1. Donc, au lieu de résoudre un système de la forme AX = b (système mal conditionné), on résoud le système équivalent DAX = Db qui est bien conditionné. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Exemple motivant Exemple motivant HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Exemple motivant Exemple motivant Une entreprise de fabrication de souvenirs veut organiser la production de trois types de souvenirs S1 , S2 et S3 sur trois machines M1 , M2 et M3 disponibles respectivement que pendant 5 heures, 3 heures et 4 heures, tels que - S1 nécessite 1 minute sur M1 , et 2 minutes sur M2 et sur M3 , - S2 nécessite 2 minutes sur M1 , et 1 minute sur M2 et sur M3 , - S3 nécessite 2 minutes sur M1 et sur M3 , et 1 minute sur M2.. Le directeur de la production veut déterminer x1 (resp. x2 et x3 ) le nombre de souvenirs de type S1 (resp. S2 et S3 ) que l’entreprise doit fabriquer pour utiliser tout le temps disponible sur les 3 machines. Le tableau suivant résume les données S1 S2 S3 Temps diponible (en minutes) Temps nécessaire sur M1 1 2 2 300 Temps nécessaire sur M2 2 1 1 180 Temps nécessaire sur M3 2 1 2 240 Ce problème peut être donc formalisé à l’aide d’un système linéaire AX = b, avec       1 2 2 x1 300 A= 2 1 1  , X =  x2  et b =  180 . 2 1 2 x3 240 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Exemple motivant 1- Appliquer la méthode de Gauss, pour trouver une solution au problème posé. 2- Déduire une factorisation PLU de la matrice du système. 3- Est ce que la matrice du système admet une factorisation de Cholesky ? Justifier. 4- On généralise au cas de 50 types de souvenirs sur 50 machines disponibles pour des temps limités. Supposons que vous disposez d’un ordinateur qui effectue un milliard d’opérations par seconde et que la matrice du système associé est inversible, a) estimer le temps théorique d’exécution en année (resp. en secondes !) nécéssaire pour résoudre le nouveau problème par la méthode de Cramer (resp de Gauss). b) Si la matrice du système est symétrique et que l’exécution de l’instruction ”chol” de matlab sur cette matrice se fait sans erreur, proposer un autre algorithme qui peut donner la solution en un temps plus réduit (à estimer en secondes). Justifier HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Exemple motivant 1- On a par la procédure d’élimination de Gauss :     1 2 2 300 l1 1 2 2 300 l1  2 1 1 180  l2 ∼  0 −3 −3 −420  l2 − 2l1 2 1 2 240 l3 0 −3 −2 −360 l3 − 2l1   1 2 2 300 l1 ∼ 0 −3 −3 −420  l2 − 2l1.Donc, AX = b ⇔ X =t (20 80 60) 0 0 1 −360 l3 − l1 2- On en déduit que     1 0 0 1 2 2 P = I3 , L= 2 1 0  et U =  0 −3 −3 . 2 1 1 0 0 1 3- A n’admet pas une factorisation de Cholesky car elle n’est pas définie positive (det ∆2 < 0). 4- a) Comme n = 50 est grand, l’algorithmes de Cramer (resp. de Gauss) demande l’équivalent de 50(51!) opérations (resp. 32 503 opérations), ainsi il 2 (50)3 50(51)! 48 demande ≈ 3600×24×365×10 9 ≈ 245 × 10 années (resp. ≈ 3 109 ≈ 0, 00008 s). b) Comme la matrice est symétrique et définie positive (car chol s’exécute sans erreur sur matlab), alors on propose la méthod de Cholesky qui demande 1 (50)3 ≈ 3 109 ≈ 0, 00004 s. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Annexe 1 : La factorisation QR par Householder Définition On appelle matrice ou réflexion de Householder relative à un vecteur v ∈ Rn r {0}, la matrice Hv définie par. 2 Hv = In − tv v tv. v Convention : On considéra que la matrice identité est de Householder. Exemple   0 −1 Si a =t (1, 1) alors Ha =. −1 0 Proposition On a Les matrices de Householder sont symétriques et orthogonales. La matrice Hv − In est de rang 1. Hλv = Hv tout scalair λ 6= 0. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder La factorisation QR par Householder (suite) Proposition Soient e1 le premier vecteur de la base canonique de Rn et a ∈ Rn r {0} avec Pn |ai | = 6 0. Alors, si on note v+ = a + kak2 e1 et v− = a − kak2 e1 , on a i=2 Hv+ a = − kak2 e1 et Hv− a = kak2 e1 , c’est à dire − kak2     0       Hv+ a =  ,  ..   .       0  + kak2     0        Hv− a =  .    ..     .  0 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Les étapes de la factorisation QR par Householder Étape 1 On pose A = A1 et on note a1 la premiere colonne de A. Si a1 = 0 on prend H1 = In , si non, on pose v1 = a1 + ka1 k2 e1 ou v1 = a1 − ka1 k2 e1 ) ce qui implique − ka1 k2 ka1 k2      0   0  Hv1 a =  ou Hv1 a =  .    .. .. .  .  0 0 On note alors H1 = Hv1 et A2 = H1 A1 par suite on obtient ∗ ∗ ∗  ...  0 ∗... ∗  A2 =   ......  ...  0 ∗... ∗ HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Étape 2 Notons maintenant a2 le vecteur de Rn−1 obtenu en prenant les (n − 1) dernières composantes de la deuxième colonne de A2. Si a2 = 0, on prend H2 = In. Si non, c-à-d : si a2 6= 0, on sait déterminer un vecteur v2 non nul tel que − ka2 k2 ka2 k2      0   0  Hv2 a =  ou Hv2 a =  .    .. .. .  .  0 0 On note alors   1 0 H2 = 0 Hv2. Ainsi, en posant A3 = H2 A2 = H2 H1 A1 , il vient ∗ ∗ ∗ ∗  ...  0 ∗ ∗... ∗    A2 =  0 0  ∗... ∗   ...... .... ..  0 0 ∗... ∗ HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Étape finale On peut itére ce procédé (n − 1) fois pour construire (n − 1) matrices H1 ,... , H(n−1) symétriques orthogonales telles que H(n−1)... H1 A1 soit une matrice triangulaire supérieure R = H(n−1)... H1 A. Donc, A = (H(n−1)... H1 )−1 R Comme les Hi sont symétriques orthogonales, on en déduit (H(n−1)... H1 )−1 = (H1 )−1... (Hn−1 )−1 t = H1... t Hn−1 = H1... Hn−1 Ainsi, si on pose Q = H1... Hn−1 on obtient A = QR avec R une matrice triangulaire supérieure et Q une matrice orthogonale (car Q −1 = (Hn−1 )−1... (H1 )−1 = t Hn−1... t H1 = t Q). HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Résultat général de la factorisation QR Ici on considère le cas général (matrices carrées ou non carrées). Soit A ∈ Mn,m (R) (avec n ≥ m). Alors, il existe une matrice orthogonale Q ∈ Mn,n (R) et une matrice trapézoı̈dale supérieure R ∈ Mn,m (R) dont les lignes sont nulles à partir de (m+1)ème c’est à dire sous la forme ∗ ∗ ∗... ∗    0 ∗ ∗... ∗  0 0 ∗... ∗    ........      ....     0 0 0... ∗     0 0 0... 0   ........  ....  0 0 0... 0 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Exemple de factorisation QR par Householder Déterminer la factorisation QR de la matrice − 46 − 43   1 5 5 A= 2 8 23 . −2 − 28 5 26 5 Étape 1 : On pose A1 = A et a1 = t (1 2 −2) et prenons (on travaille dans R3 ) v1 = a1 − ka1 k2 e1 t = (−2 2 − 2). Alors 2 H1 ≡ Hv1 ≡ I3 − tv v1 t v1 1 · v1   −2 2 = I3 −    2  (−2 2 − 2) −2 −2 (−2 2 − 2)  2  −2  1 2 2  3 3 −3 =  2 1 2 . 3 3 3 − 23 32 1 3 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Ainsi, A2 ≡ H1 A1  1 2 − 32 − 46 − 43   3 3 1 5 5 2 1 2  2 =  3 3 3 8 23  − 23 2 3 1 3 −2 − 28 5 26 5   3 6 9 =  0 − 36 5 27 5 . 48 114 0 5 5 Étape 2 : On pose a2 = t (− 36 5 48 5 ) et prenons (on travaille dans R2 ) v2 = a2 + ka2 k2 e1 t 24 48 = (5 5 ). Alors, 24   2 2 Hv 2 ≡ I2 − tv · v v2 t v2 = I2 −  24  5 48 ( 24 5 44 5 ) 2 2 5 ( 24 5 48 5 ) 5 48 5 3 − 45     2 1 = I2 − (1 2) = 5. 2 − 45 − 35   1 (1 2) 2 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Ainsi,   1 0 0 3 H2 =  0 5 − 54 . 0 − 45 − 35 et donc A3 ≡ H2 A2    1 0 0 3 6 9 3 4 =  0 5 −5  0 − 36 5 27 5  0 − 45 − 53 0 48 5 114 5   3 6 9 =  0 −12 −15 . 0 0 −18 On en déduit que A = QR avec   3 6 9 R ≡ H2 H1 A =  0 −12 −15  0 0 −18 Q = (H2 H1 )−1 = H1−1 H2 −1 = H1 H2  1 14 2  3 15 − 15 =  2 − 13 − 23 . 3 2 2 −3 15 − 11 15 HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder De plus, si par exemple le système AX = t (1 0 0) devient t Ax = (1 0 0) ⇔ QRx = t (1 0 0) ⇔ Rx = Q t (1 0 t 0)  1  3 14 ⇔ Rx =  15  2 − 15 71 47 1
⇔ x= 270 − 540 135. HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025  Annexe 1 : La factorisation QR par Householder Théorème Soit A ∈ Mn (R). On peut trouver une factorisation QR telle que tous les coefficients diagonaux de R sont tous positifs ou nul. Dans ce cas, si A est inversible, cette factorisation est unique. L’idée de la démonstration Soit A ∈ Mn (R). Alors, d’aprés ce qui précède il existe une matrice R 0 triangulaire supérieure et une matrice Q 0 orthogonale telles que A = Q 0 R 0. Notons la matrice D = diag (d1 ,... , dn ) avec ( rii si rii = 0, di = |rii | 1 si rii = 0. Posons Q = Q 0 D −1 et R = DR 0 , alors    A = QR, Q est orthogonale,    R est triangulaire supérieure, rii ≥ 0, ∀i ∈ {1,... n}.  HADDADOU Hamid Méthodes directes de résolution des systèmes linéaires-ESI 2024/2025