Metode Numerice - Curs 3 PDF

METODE NUMERICE – curs 3 Cap. 2 Sisteme determinate de ecuaţii algebrice liniare 2.1 Formularea problemei  Metodele numerice de rezolvare a unui sistem determinat de ecuaţii algebrice liniare metode directe: metode iterative...

METODE NUMERICE – curs 3 Cap. 2 Sisteme determinate de ecuaţii algebrice liniare 2.1 Formularea problemei  Metodele numerice de rezolvare a unui sistem determinat de ecuaţii algebrice liniare metode directe: metode iterative (indirecte):  sistem echivalent, direct rezolvabil prin mijloace  construirea unui şir de aproximaţii pentru elementare → eliminarea progresivă a soluţia sistemului, convergent la aceasta: necunoscutelor (eliminare gaussiană) 𝑘→∞ 𝑥 [𝑘] 𝑥  prin transformări elementare de echivalenţă, se 𝑘≥0 aduce matricea A la anumite forme tipice: u 11   u 1n   l11 0  0  [k]  0 lim || x − x ||  = 0         k → U= L=          0      calculele se opresc la un index de iterare  0  0 u nn  l  n1   l nn  [s], când este îndeplinită condiţia: [s ] [s −1] matrice superior matrice inferior || x −x ||   impus triunghiulară triunghiulară  procedura de transformare se numeşte triangularizare  METODE NUMERICE – curs 3 2.2 Rezolvarea sistemelor prin triangularizare directă 2.2.1 Principiul metodei  se consideră sistemul de n ecuaţii algebrice liniare cu n necunoscute: 𝐴 ∙ 𝑥 = 𝑏, 𝐴 ∈ ℝ𝑛×𝑛 , 𝑏 ∈ ℝ𝑛×1  notaţie: A = [a ij ]1i  n 1 j n  matricile: [a 11 ] , …, [a ij ]1i n −1 - submatrici principale (minori principali directori) ale lui A 1 j n −1 Teoremă: Dacă matricea 𝐴 ∈ ℝ𝑛×𝑛 are toate submatricele principale inversabile (nesingulare), atunci există matricele 𝐿, 𝐷, 𝑈 ∈ ℝ𝑛×𝑛 astfel încât: A =LDU unde L este o matrice inferior triunghiulară, D este o matrice diagonală şi U este o matrice superior triunghiulară.  METODE NUMERICE – curs 3 Observaţii: 1. uzuale sunt factorizările: 𝐴 = 𝐿 ∙ 𝑈 (factorizare L-U); 2. demonstraţia teoremei enunţate este constructivă, constituind însuşi algoritmul de descompunere L-U a matricei A; 3. algoritmul de descompunere → procedeul de eliminare gaussiană a necunoscutelor - matricea A este adusă la forma superior triunghiulară în urma unui şir de transformări de asemănare → se “acumulează” într-o matrice inferior triunghiulară, cu elementele de pe diagonala principală egale cu 1 (descompunere Doolittle); 4. considerând descompunerea L-U a matricei sistemului A, rezolvarea sistemului implică două subetape: a) substituţie înainte: rezolvarea sistemului 𝐿 ∙ 𝑦 = 𝑏 - din aproape în aproape: 𝑦1 din prima ecuaţie → 𝑦2 din a doua ecuaţie → … 𝑦𝑛 din ultima ecuaţie (𝑦 = 𝑦𝑖 𝑖=1,𝑛 ) b) substituţie inversă: rezolvarea sistemului U∙ 𝑥 = 𝑦 - din aproape în aproape: 𝑥𝑛 din ultima ecuaţie → 𝑥𝑛−1 din penultima ecuaţie → … 𝑥1 din prima ecuaţie (𝑥 = 𝑥𝑖 𝑖=1,𝑛 )  METODE NUMERICE – curs 3 Propoziţie: Dacă matricea A admite o descompunere L-U, atunci această descompunere este unică.  procedura de triangularizare directă necesită un număr de operaţii în virgulă mobilă de 3 ordinul lui 𝑛 ൗ3;  dacă matricea A este simetrică (𝐴 = 𝐴𝑇 ) şi pozitiv definită (∀𝑥 ∈ ℝ𝑛×1 , 𝑥 ≠ 0𝑛×1 , 𝑥 𝑇 ∙ 𝐴 ∙ 𝑥 > 0 şi → 𝑥 𝑇 ∙ 𝐴 ∙ 𝑥 = 0 ⇔ 𝑥 = 0𝑛×1 ), A se descompune sub forma 𝐴 = 𝐿 ∙ 𝐿𝑇 (descompunerea Cholesky). 2.2.2 Procedura de triangularizare directă a unei matrice  algoritmul triangularizării directe: atribuie 𝐴1 ⟵ 𝐴 pentru 𝑘 = 1, 𝑛 − 1 execută  determinare matrice 𝑀𝑘 astfel încât 𝐴𝑘+1 = 𝑀𝑘 ∙ 𝐴𝑘 să aibă elementele: [𝑘+1] [𝑘+1] [𝑘] 𝑎𝑖,𝑘 = 0, 𝑖 = 𝑘 + 1, … , 𝑛 şi 𝑎𝑖,𝑗 = 𝑎𝑖,𝑗 , 𝑖 = 1, … , 𝑛; 𝑗 = 1, … , 𝑘 − 1 atribuie 𝐴𝑘+1 ⟵ 𝑀𝑘 ∙ 𝐴𝑘  atribuie 𝑈 ⟵ 𝐴𝑛  METODE NUMERICE – curs 3  algoritmul parcurge (n – 1) etape în care se zerorizează elementele unei coloane, situate sub diagonala principală, lăsând nemodificate coloanele transformate la etapele anterioare [𝑘]  notaţie: 𝜉 = 𝑐𝑘 𝐴𝑘 ; 𝑐𝑘 𝐴𝑘 - coloana k a matricii 𝐴𝑘 = 𝑎𝑖,𝑗 𝑖,𝑗=1,…,𝑛  = [1   k  k +1   n ]T = [a 1[ k,k]  a [kk,k] a [kk+]1,k  a [nk,k] ]T 𝑇  𝑀𝑘 se construieşte încât 𝑀𝑘 ∙ 𝜉 = 𝜉1 ⋯ 𝜉𝑘 0 ⋯ 0  se consideră vectorul de multiplicatori: 𝑚𝑘 = 0 ⋯ 0 𝜇𝑘+1,𝑘 ⋯ 𝜇𝑛,𝑘 𝑇 Definiţie: Matricea 𝑀𝑘 se numeşte matrice de transformare elementară de ordin n şi indice k sau matrice Gauss şi este definită prin: T Mk = In − mk  ek unde 𝑒𝑘 este coloana k a matricii unitate de ordin n, 𝐼𝑛. CALCUL NUMERIC – curs 3  proprietăţi ale matricilor Gauss: sunt inferior triunghiulare; elementele diagonalei principale sunt egale cu 1; sunt nesingulare, având inversa: M −k1 = I n + m k  e k T  efectul aplicării matricei 𝑀𝑘 asupra vectorului 𝜉: T T M k   = (I n − m k  e k )   =  − m k  e k   =  − m k   k = [ 1   k  k +1 −  k +1,k   k   n −  n ,k   k ]T trebuie zerorizate  presupunând 𝜉𝑘 ≠ 0, se aleg multiplicatorii Gauss:  i,k =  i /  k , i = k + 1,..., n pivot M k   = [1   k 0  0]T  METODE NUMERICE – curs 3 Observaţii: 1. În practică, pe calculator, etapa k descrisă mai sus se poate realiza testând condiţia: | a [kk,k] |  [𝑘] în loc de a verifica 𝑎𝑘,𝑘 ≠ 0, unde 𝜀 este o constantă impusă, de valoare mică sau foarte mica (de exemplu, epsilonul-maşină); aceasta se realizează datorită faptului că, dacă în aritmetica exactă pivotul este nul, în aritmetica virgulei mobile, datorită erorilor de calcul, această situaţie [k] este echivalentă cu: | a k ,k |  ; 2. Când pivotul este în modul mai mic sau egal cu 𝜀, eliminarea gaussiană eşuează; aceasta corespunde situaţiei când matricea iniţială A are submatricea principală de ordin k singulară, deci conform teormei enunţate anterior, descompunerea L-U a matricei A nu există.  efectul aplicării transformării 𝑀𝑘 asupra celorlaltor coloane ale matricei 𝐴𝑘 : 𝑇 - se consideră un vector 𝜂 ≠ 𝜉, cu 𝜂 = 𝜂1 ⋯ 𝜂𝑛 M k   = [1  k k +1 −  k +1,k  k  n −  n,k  k ]T METODE NUMERICE – curs 3 Concluzii: 1. Matricea 𝑀𝑘 lasă nemodificate primele k – 1 coloane ale matricei 𝐴𝑘 :  = c j (A k ) = [*  * 0  0  0]T , j  k, j = 1,..., k − 1 1 j j+1 k n M k  c j (A k ) = [*  * 0  0 0 −  k +1,k  0  0 −  n ,k  0]T 1 j j+1 k k +1 n = c j (A k ) 2. Matricea 𝑀𝑘 transformă coloana k a matricei 𝐴𝑘 zerorizând liniile 𝑘 + 1, ⋯ , 𝑛; 3. Matricea 𝑀𝑘 transformă coloanele 𝑘 + 1, ⋯ , 𝑛 ale matricei 𝐴𝑘 în liniile 𝑘 + 1, ⋯ , 𝑛:  = c j (A k ), j = k + 1,..., n M k  c j (A k ) = [*  * a [kk+]1, j −  k +1,k  a [kk, ]j  a [nk, ]j −  n ,k  a [kk, ]j ]T (notaţia * semnificând faptul că elementele implicate rămân nemodificate) METODE NUMERICE – curs 3  tabloul general al transformărilor: M n −1    M 2  M 1  A = U A = M1−1  M −21    M −n1−1  U A =LU L n −1 L = M 1−1  M −21    M −n1−1 = I n +  m k  e k T k =1 Observaţie: 1. Matricea L este inferior triunghiulară unitate şi conţine în fiecare coloană, sub elementul unitar de pe diagonala principală, subvectorii Gauss; 2. Prima sub-etapă de rezolvare a sistemului 𝐴 ∙ 𝑥 = 𝑏 este substituţia înainte aplicată sistemului de ecuaţii 𝐿 ∙ 𝑦 = 𝑏 : y = L−1  b = M n −1    M 2  M1  b  METODE NUMERICE – curs 3 Concluzie: La triangularizarea simplă, unde elementele matricei A se modifică corespunzător relaţiei: a [ijk +1] = a [ijk ] −  ik  a [kjk ] , i = k + 1,..., n; j = k,..., n multiplicatorii 𝜇𝑖𝑘 pot avea, în principiu, orice valoare - dacă aceste valori sunt mari sau foarte mari, atunci: pot apare fenomenele de omitere catastrofală şi/sau de neutralizare a termenilor; [𝑘] amplifică erorile prezente în termenii 𝑎𝑘𝑗. triangularizarea simplă este instabilă numeric  stabilizarea algoritmului triangularizării unei matrici → multiplicatori Gauss cu modul subunitar  METODE NUMERICE – curs 3 𝟐 −3 5 Exemplu: 𝐴 = 𝐴1 = 4 7 −3 1 9 4 iteraţia 1: 0 0 1 0 0 𝟐 −3 5 𝑚1 = 4Τ𝟐 = 2 𝑀1 = −2 1 0 𝐴2 = 𝑀1 ∙ 𝐴1 = 0 𝟏𝟑 −13 1Τ𝟐 0.5 −0.5 0 1 0 10.5 1.5 iteraţia 2: 0 0 1 0 0 𝑚2 = 0 = 0 𝑀2 = 0 1 0 10.5Τ𝟏𝟑 0.8077 0 −0.8077 1 𝟐 −3 5 𝐴3 = 𝑀2 ∙ 𝐴2 = 0 𝟏𝟑 −13 = 𝑈 0 0 12 1 0 0 𝐿= 𝑀1−1 ∙ 𝑀2−1 = 2 1 0 0.5 0.8077 1 METODE NUMERICE – curs 3 2.3 Rezolvarea sistemelor prin triangularizare cu pivotare parţială  principiu: la pasul k se caută pivotul printre elementele din coloana k, pornind de la elementul de pe diagonala principală în jos, alegându-se elementul care are cea mai mare valoare în modul: | a [i kk ],k |= max{| a [i,kk] |} =  k k i  n ⎯k Ak = 0 ⎯ ik ⎯n | k  dacă 𝑖𝑘 ≠ 𝑘 → se permută liniile k şi 𝑖𝑘 cu ajutorul matricii de permutare de linii 𝑃𝑘 Pk  A k A k +1 = M k  (Pk  A k ) |  i,k | 1, i = k + 1,..., n METODE NUMERICE – curs 3  proprietăţi matrice de permutare de linii, Pk: det(Pk ) = −1; Pk = Pk−1 Teoremă: Dacă matricea A este nesingulară, atunci există o matrice numită matrice generală de permutare de linii, astfel încât: P  A = L'  U în care U este o matrice superior triunghiulară şi Lˈ este o matrice inferior triunghiulară unitate cu elementele | l i, j | 1, i  j  Demonstraţia → constructivă, constituind însuşi algoritmul de triangularizare cu pivotare parţială a unei matrici METODE NUMERICE – curs 3  algoritmul de triangularizare cu pivotare parţială: atribuie 𝐴1 ⟵ 𝐴 pentru 𝑘 = 1, 𝑛 − 1 execută [𝑘] [𝑘] [𝑘] determinare 𝜉𝑘 ← 𝑎𝑖𝑘 ,𝑘 unde 𝑎𝑖𝑘 ,𝑘 = max 𝑎𝑖,𝑘 𝑖=𝑘,𝑛 dacă (𝑖𝑘 ≠ 𝑘) atunci determină 𝑃𝑘 altfel atribuie 𝑃𝑘 ← 𝐼𝑛  calcul 𝑃𝑘 ∙ 𝐴𝑘 determinare matrice 𝑀𝑘 astfel încât 𝐴𝑘+1 = 𝑀𝑘 ∙ 𝑃𝑘 ∙ 𝐴𝑘 să aibă elementele: [𝑘+1] [𝑘+1] [𝑘] 𝑎𝑖,𝑘 = 0, 𝑖 = 𝑘 + 1, … , 𝑛 şi 𝑎𝑖,𝑗 = 𝑎𝑖,𝑗 , 𝑖 = 1, … , 𝑛; 𝑗 = 1, … , 𝑘 − 1 atribuie 𝐴𝑘+1 ⟵ 𝑀𝑘 ∙ 𝑃𝑘 ∙ 𝐴𝑘  atribuie 𝑈 ⟵ 𝐴𝑛  METODE NUMERICE – curs 3  tabloul general al transformărilor: M n −1  Pn −1    M 2  P2  M 1  P1  A = U - folosind 𝑃𝑘 = 𝑃𝑘−1 (M n −1  Pn −1    M 1  P1  P1  P2    Pn −1 )  (Pn −1    P2  P1 )  A = U −1 P 𝐿′ matrice generală de permutare de linii P  A = L'  U L' = Pn −1    P2  M1−1  P2  M −21    Pn −1  M −n1−1 matrice inferior triunghiulară unitate având în fiecare coloană, sub diagonala principală, subvectori Gauss cu liniile permutate METODE NUMERICE – curs 3  etapele rezolvării sistemului 𝐴 ∙ 𝑥 = 𝑏, A ∈ ℝ𝑛×𝑛 , 𝑏 ∈ ℝ𝑛×1 : a) factorizarea L_U a matricii A prin triangularizare cu pivotare parţială: P  A = L'  U b) calculul vectorului: 𝑐 = 𝑃 ∙ 𝑏 c) rezolvare sistem 𝐿′ ∙ 𝑦 = 𝑐 prin substituţie înainte d) Rezolvare sistem 𝑈 ∙ 𝑥 = 𝑦 prin substituţie inversă PAx =Pb L'  U  x = P  b P  A = L'  U c=Pb L'  y = c y=Ux  Observaţie: Dacă algoritmul de triangularizare cu pivotare parţială eşuează, în sensul că pivotul găsit la o anumită etapă [k] este nul sau foarte mic în modul, aceasta corespunde situaţiei când în aritmetica reală primele k coloane ale matricei A sunt liniar dependente.