Metodi numerici per la risoluzione di Sistemi Lineari
by user
Comments
Transcript
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti Metodi Iterativi Metodi numerici per la risoluzione di Sistemi Lineari Stefano Berrone Dipartimento di Matematica tel. 011 0907503 [email protected] http://calvino.polito.it/~sberrone Laboratorio di modellazione e progettazione materiali Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Problemi Triangolari - Sostituzione in avanti Vogliamo risolvere il sistema lineare Ly = b con L triangolare inferiore (lower triangular) non singolare. b1 y1 l11 0 0 ... 0 l21 l22 0 ... 0 y2 b2 l31 l32 l33 ... 0 y3 b3 = .. .. .. .. . . . . yn ln1 ln2 ln3 ... lnn bn Osservazione Se L è non singolare, tutti gli elementi diagonali sono diversi da 0: det(L) = n Y lii i=1 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Esplicitando le equazioni: l11 y1 = b1 , l21 y1 + l22 y2 = b2 , .. . ln1 y1 + ln2 y2 + ... + lnn yn = bn . =⇒ y1 = b1 /l11 =⇒ y2 = (b2 − l21 y1 )/l22 =⇒ yi = (bi − i−1 X lij yj )/lii j=1 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodo di sostituzione in avanti (o forward substitution): for i=1:n y(i)=b(i); for j=1:i-1 y(i)=y(i)-L(i,j)*y(j); end y(i)=y(i)/L(i,i); end Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Un’alternativa che sfrutta la capacità di MATLAB di operare su sottomatrici è la seguente: y(1) =b(1)/L(1,1); for i=2:n y(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i); end Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Sostituzione all’indietro Vogliamo risolvere il sistema lineare Ux = y con U triangolare superiore (upper triangular) non singolare. u11 u12 u13 ... u1n 0 u22 u23 ... u2n .. .. . . 0 0 ... un−1,n−1 un−1,n 0 0 ... 0 unn x1 x2 x3 .. . = xn y1 y2 y3 .. . . yn Osservazione Di nuovo, dalla nonsingolarità di U segue che uii 6= 0∀i Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi u11 x1 + u12 x2 + u22 x2 + ... ... .. . + + u1n xn u2n xn = = y1 , y2 , un−1,n−1 xn−1 + un−1,n xn = yn−1 , unn xn = yn . =⇒ xn = yn /unn . =⇒ xn−1 = (yn−1 − un−1,n xn )/un−1,n−1 . =⇒ xi = (yi − n X uij xj )/uii . j=i+1 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodo di sostituzione all’indietro (o backward substitution): for i=n:-1:1 x(i)=y(i); for j=i+1:n x(i)=x(i)-U(i,j)*x(j); end x(i)=x(i)/U(i,i); end Sfruttando le funzioni ottimizzate della libreria BLAS: x(n)=y(n)/U(n,n); for i=n-1:-1:1 x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i); end Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Eliminazione Gaussiana La risoluzione di sistemi triangolari è quindi un problema particolarmente semplice da trattare. Idea: Nel caso di sistemi non triangolari, cercare di ricondursi ad uno o più sistemi equivalenti di forma triangolare. Dato un sistema lineare Ax = b le seguenti operazioni conducono ad un sistema equivalente (ovvero con la stessa soluzione): 1 scambiare due equazioni; 2 moltiplicare un’equazione per uno scalare; 3 sostituire un’equazione con una combinazione lineare della stessa e di un’altra equazione. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodo di eliminazione di Gauss: consiste nell’applicare ripetutamente l’operazione 3 in modo da: 1 trasformare il sistema assegnato in un sistema a matrice triangolare superiore 2 e quindi nel risolvere il sistema equivalente cosı̀ ottenuto con una backward substitution. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Scriviamo il sistema Ax = b esplicitamente: a11 x1 + a12 x2 + ... + a1n xn a21 x1 + a22 x2 + ... + a2n xn a31 x1 + a32 x2 + ... + a3n xn .. . an1 x1 + an2 x2 + ... + ann xn = b1 , = b2 , = b3 , = bn . Supponiamo a11 6= 0 r2 ← r2 + (− Stefano Berrone a21 ) ∗ r1 a11 Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting a21 x1 + a22 x2 + ... + a2n xn = b2 + 21 − aa11 ( a11 x1 + a12 x2 + ... + a1n xn = b1 ) 0x1 (2) (2) (2) (2) + a22 x2 + ... + a2n xn = b2 . (2) a22 = a22 − a12 a21 /a11 , ...,a2n = a2n − a1n a21 /a11 , (2) b2 = b2 − b1 a21 /a11 Nella seconda equazione non compare più la prima incognita. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Procediamo con tutte le ultime n − 1 equazioni: ri ← ri + mi1 ∗ r1 , i = 2, ..., n ai1 mi1 = − , i = 2, ..., n. a11 (1) a11 x1 0 0 =⇒ 0 (1) (1) aij = aij , bi (2) aij = (1) aij (1) (1) (1) (2) (2) (2) + a12 x2 + ... + a1n xn = b1 (2) (2) (2) + a22 x2 + ... + a2n xn = b2 (2) (2) (2) + a32 x2 + ... + a3n xn = b3 .. . + an2 x2 + ... + ann xn = bn = bi , i, j = 1, ..., n; (1) (2) + mi1 a1j , bi (1) = bi Stefano Berrone (1) + mi1 b1 , i, j = 2, ..., n. Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Seconda colonna: (2) supponiamo a22 6= 0. Eliminiamo x2 dalle ultime n − 2 equazioni trasformate: (2) ri ← ri + r2 ∗ mi2 , (1) a11 x1 0 0 =⇒ 0 (3) (2) mi2 = − (1) ai2 (2) i = 3, ..., n a22 (1) (1) (1) (3) (3) (3) + a12 x2 + a13 x3 + ... + a1n xn = b1 (2) (2) (2) (2) + a22 x2 + a23 x3 + ... + a2n xn = b2 (3) (3) (3) + 0 + a33 x3 + ... + a3n xn = b3 .. . + + an3 x3 + ... + ann xn = bn 0 (2) (3) aij = aij + mi2 a2j , bi (2) = bi Stefano Berrone (2) + mi2 b2 , i, j = 3, ..., n Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Dopo n − 1 passi: (1) (1) a11 x1 + a12 x2 (2) + a22 x2 0 0 + 0 0 + 0 (1) Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting (1) (1) (n) (n) + a13 x3 + ... + a1n xn = b1 (2) (2) (2) + a23 x3 + ... + a2n xn = b2 (3) (3) (3) + a33 x3 + ... + a3n xn = b3 .. . + 0 + ... + ann xn = bn Osservazione Il termine noto b ha subito esattamente tutte le trasformazioni subite dagli elementi della matrice A. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Definizione (pivot e moltiplicatori) (i) L’elemento a11 e i successivi elementi aii sono detti elementi pivot. I coefficienti mij sono detti moltiplicatori. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodo di eliminazione di Gauss e fattorizzazione LU In notazione matriciale, il k-esimo passo di eliminazione Gaussiana puó essere rappresentato come segue: A(k+1) = Mk A(k) con 1 1 mk+1,k Mk = .. . mn,k Stefano Berrone 1 .. . 1 Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Dopo n − 1 passi: A(n) = Mn−1 A(n−1) = Mn−1 Mn−2 A(n−2) = . . . = Mn−1 . . . M2 M1 A | {z } M b (n) = Mn−1 . . . M2 M1 b e A(n) x = b (n) con A(n) triangolare superiore ⇓ Ax = b ⇔ MAx = Mb Posto U = MA, si può scrivere Ux = Mb da cui ottengo la soluzione con sostituzione all’indietro. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Esempio Dati 4 4 8 A = 2 8 7 1 3 6 Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting 12 b=9 7 calcolare la soluzione con eliminazione gaussiana. 1 0 0 M1 = − 12 1 0 − 14 0 1 e 4 4 8 | 12 M1 ∗ (A|b) = 0 6 3 | 3 . 0 2 4 | 4 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Esempio 1 0 M2 = 0 1 0 − 31 e 0 0 1 4 4 8 | 12 M2 M1 ∗ (A|b) = 0 6 3 | 3 = (U|b (3) ). 0 0 3 | 3 Con sostituzione all’indietro: x3 = 1 x2 = (3 − 3 ∗ 1)/6 = 0 x1 = (12 − 4 ∗ 0 − 8 ∗ 1)/4 = 1 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Com’è fatta M? - triangolare inferiore, ma 1 0 0 m21 1 0 m31 m32 1 M 6= .. .. .. . . . mn1 mn2 mn3 Poniamo ... 0 . . . 0 . . . 0 .. .. . . ... 1 L = M −1 U = MA =⇒ LU = LMA = A Quindi il sistema di partenza equivale a LUx = b. Com’è fatta L? Suo ruolo? Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Struttura di L L triangolare inferiore. Il calcolo è meno complicato e costoso di quanto sembri a prima vista: non occorre costruire M e poi calcolarne l’inversa. 1 0 0 ... 0 −m21 1 0 ... 0 ... 0 L = −m31 −m32 1 .. . . . . −mn1 −mn2 . . . −mn,n−1 1 Osservazione L può quindi essere costruita durante l’eliminazione gaussiana semplicemente conservando gli elementi mij e cambiandogli il segno. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Fattorizzazione LU Abbiamo fattorizzato A = LU con L triangolare inferiore U triangolare superiore Utilizzo per risolvere il sistema lineare: Ax = b LUx = b L |{z} Ux = b y ⇓ Ly = b, Stefano Berrone Ux = y Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Il problema di partenza può essere affrontato in due fasi, risolvendo due sistemi a matrice triangolare: 1 prima determiniamo y risolvendo Ly = b; 2 una volta noto y , determiniamo x risolvendo Ux = y . Osservazione La risoluzione del sistema lineare Ly = b è del tutto equivalente ad effettuare sul vettore b le trasformazioni dell’eliminazione gaussiana. Quindi il vettore y calcolato come soluzione del sistema lineare Ly = b coincide con b (n) ! Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Meglio fattorizzazione LU + 2 sistemi lineari triangolari, o eliminazione gaussiana senza costruire esplicitamente L? Con l’eliminazione Gaussiana, b subisce le stesse trasformazioni di A. Su un solo sistema sono equivalenti. Se devo risolvere un secondo sistema lineare Ax = c con eliminazione gaussiana devo rieseguire tutto il procedimento per effettuare le trasformazioni su c e ottenere c (n) ; con fattorizzazione LU già calcolata devo solo risolvere i sistemi Ly = c e Ux = y . È veramente più conveniente fare cosı̀ piuttosto che rieffettuare ex-novo l’eliminazione gaussiana sul sistema Ax = c? Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Costo computazionale Per costo computazionale di un algoritmo intendiamo l’ordine di grandezza del numero di moltiplicazioni o divisioni effettuate, calcolato in funzione della dimensione n del problema (in questo caso l’ordine della matrice). Costo computazionale per calcolo fattorizzazione LU (coincidente con quello della sola eliminazione gaussiana): 3 n O moltiplicazioni. 3 Il costo computazionale soluzione di un sistema triangolare (superiore o inferiore indifferentemente): 2 n O moltiplicazioni. 2 Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Risolvendo un sistema lineare con fattorizzazione LU la maggior parte del costo risiede nella fattorizzazione e non nella soluzione dei sistemi triangolari ottenuti: il costo totale è infatti di O( n3 n3 ) + O(n2 ) ' O( ) 3 3 operazioni, essendo (per n grande) n2 trascurabile rispetto a n3 . Osservazione I metodi di fattorizzazione sono quindi particolarmente convenienti quando si debbano risolvere un certo numero di sistemi lineari in cui la matrice dei coefficienti è sempre la stessa e cambia solo il vettore dei termini noti. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Esempio Data la matrice Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting 4 4 8 A = 2 8 7 1 3 6 calcolarne la fattorizzazione LU. Basta recuperare i moltiplicatori: 4 4 8 U = 0 6 3 0 0 2 1 L = 12 1 4 Osservazione: det(A) = det(L) · det(U) = 1 · Stefano Berrone 0 0 1 0 1 3 1 Qn i=1 uii = 48. Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Osservazione L può essere memorizzata nella stessa matrice A via via che in questa vengono azzerati gli elementi sotto la diagonale. 4 4 8 4 4 4 8 2 8 7 −→ 1 6 3 −→ 1 2 2 1 1 1 3 6 4 2 4 4 Stefano Berrone 4 8 6 3 1 3 2 Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting L’eliminazione di Gauss / fattorizzazione LU può essere sempre portata a termine con successo? Definizione Sia A ∈ Rn×n ; definiamo Ak , con k = 1, ..., n, sottomatrice principale di testa di ordine k la matrice ottenuta intersecando le prime k righe e k colonne di A. Teorema Sia A una matrice di ordine n. Se Ak è non singolare per k = 1, ..., n − 1 allora esiste ed è unica la fattorizzazione LU di A. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Osservazione Le matrici a diagonale dominante per righe o per colonne e le matrici simmetriche definite positive soddisfano le ipotesi del teorema precedente. Osservazione Il teorema precedente non esclude che la matrice A sia singolare. Se anche A è singolare, ma det(Ak ) 6= 0, k = 1, ..., n − 1, esiste comunque la fattorizzazione LU, con la matrice U che avrà un elemento uii nullo. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting La Strategia del Pivoting Esempio Che sucede se voglio applicare eliminazione gaussiana a 0 1 A= ? 1 0 Occorre ricorrere alle matrici di permutazione. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Definizione Sia Pij la matrice ottenuta dalla matrice identità permutando la i-esima e la j-esima riga: 1 0 0 0 1 0 0 0 1 i 0 0 0 Pij = j 0 0 0 0 0 0 ... ... ... ... ... ... i 0 0 0 .. . 0 .. . 1 .. . 0 ... ... ... ... ... ... j 0 0 0 .. . 1 .. . 0 .. . 0 ... ... ... ... ... ... 0 0 0 0 . 0 1 Le matrici Pij sono dette matrici di permutazione semplice. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Data una matrice A la moltiplicazione a sinistra di A per Pij (Pij A) produce lo scambio della i-esima e j-esima riga di A; la moltiplicazione a destra di A per Pij (APij ) produce lo scambio della i-esima e j-esima colonna di A. Se moltiplichiamo più matrici di permutazione semplice fra loro, otteniamo una matrice (che indicheremo con P e che chiameremo matrice di permutazione) il cui effetto su una matrice A sarà quello di effettuare tutti gli scambi di righe (colonne) associate alle matrici di permutazione semplice. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Teorema Sia A una matrice di ordine n. Allora esiste una matrice di permutazione P per cui si può ottenere la fattorizzazione LU di A, cioè: PA = LU. Si osservi che, data A, vi possono essere diverse matrici di permutazione P tali che PA ammette fattorizzazione LU. Come individuare una matrice di permutazione? ⇓ Strategia di pivoting parziale. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Pivoting Parziale Se durante l’eliminazione di Gauss al passo k-esimo risulta (k) akk = 0, scambiamo due equazioni in modo che il nuovo (k) elemento pivot akk = 0 sia diverso da zero. La scelta è sempre tra le equazioni che seguono la k-esima. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Esempio Consideriamo il seguente sistema lineare: x1 + x2 + x3 = 1 x1 + x2 + 2x3 = 2 x1 + 2x2 + 2x3 = 1 Eliminazione gaussiana: m21 = −1 e m31 = −1 =⇒ x1 + x2 + x3 = 1 0 + 0 + x3 = 1 0 + x2 + x3 = 0 (2) a22 = 0! Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Esempio (segue) x1 + x2 + x3 = 1 0 + x2 + x3 = 0 ; 0 + 0 + x3 = 1 Questo è già in forma triangolare e si può ottenere per sostituzione all’indietro la soluzione: x3 = 1 x2 = −x3 = −1 x1 = 1 − x2 − x3 = 1 Abbiamo det(A) = − det(U) = −1; in generale, det(A) = (−1)s · det(U), dove s è il numero di permutazioni di riga effettuate. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Fattorizzazione PA = LU Cosa accade alla fattorizzazione LU con il pivoting? A(k+1) = Mk Pk A(k) n − 1 passi: A(n) = Mn−1 Pn−1 A(n−1) = ... = Mn−1 Pn−1 Mn−2 Pn−2 . . . M2 P2 M1 P1 A U = Mn−1 Pn−1 Mn−2 Pn−2 Mn−3 . . . M2 P2 M1 P1 A = Mn−1 M̄n−2 Pn−1 Pn−2 Mn−3 . . . M2 P2 M1 P1 A = Mn−1 M̄n−2 . . . M̄2 M̄1 Pn−1 Pn−2 . . . P2 P1 A | {z }| {z } M P = MPA Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting U = MPA M −1 U = PA Posto L = M −1 abbiamo determinato fattorizzazione LU di PA, i.e. fattorizzazione PA = LU. Fattorizzazione PA = LU per risolvere un sistema lineare: 1 Ly = Pb 2 Ux = y Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting In generale la strategia di permutare equazioni/righe viene usata anche se si incontrano elementi pivot diversi da zero per aumentare la stabilità dell’eliminazione gaussiana. Si può infatti dimostrare che per aumentare la stabilità è opportuno evitare moltiplicatori troppo grandi, poiché amplificano gli errori di arrotondamento, con conseguente instabilità dell’algoritmo. Per prevenire la crescita dei moltiplicatori è necessario prevenire elementi pivot troppo piccoli. Per raggiungere questo obiettivo possiamo applicare la strategia del pivoting parziale ammettendo scambi di righe anche se al (k) passo k l’elemento pivot akk non è nullo. Più precisamente, al passo k si determina l’elemento di modulo (k) massimo fra tutti gli elementi ark , con r ≥ k; una volta individuato l’indice r corrispondente all’elemento di modulo massimo si scambiano fra loro la riga k-esima e r -esima. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting In questo modo ci si assicura che ad ogni passo si sta usando l’elemento pivot più grande possibile. Questa strategia incrementa in modo significativo la stabilità dell’algoritmo. Osservazione Se la matrice A è a diagonale dominante per colonne è possibile dimostrare che la strategia del pivoting parziale è del tutto superflua in quanto non provoca alcun scambio di equazioni.Inoltre, se la matrice A è simmetrica e definita positiva, il metodo di Gauss senza pivoting risulta numericamente stabile, quindi non c’è bisogno di effettuare pivoting. Per tale classe di matrici esiste una diversa fattorizzazione molto più conveniente da usare (la fattorizzazione di Cholesky). Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Come costruire L P = Pn−1 . . . P1 U = A(n) L =? È facilmente ottenibile con l’eliminazione gaussiana tenendo conto delle permutazione anche nell’immagazzinamento dei moltiplicatori sotto la diagonale. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Esempio 1 4 8 A = 2 0 7 4 2 6 P1 = P13 : 4 2 6 P13 A = 2 0 7 1 4 8 4 2 6 A(2) = 12 −1 4 1 4 Stefano Berrone 7 2 13 2 Metodi numerici per la risoluzione di Sistemi Lineari Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Metodi diretti Metodi Iterativi Esempio (segue...) P2 = P23 : 2 6 P23 A(2) = 14 7 2 13 2 1 2 −1 (3) A 1 L = 14 1 2 0 1 − 72 4 4 = 1 4 1 2 0 4 2 0 U = 0 72 1 0 0 Stefano Berrone 2 7 2 − 27 6 6 4 13 2 41 7 13 2 41 7 P = P23 P13 0 0 1 = 1 0 0 0 1 0 Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Fattorizzazione di Cholesky Se A è simmetrica e definita positiva, essendo l’eliminazione gaussiana stabile senza pivoting, si può sfruttare la simmetria di A per ridurre il costo computazionale del calcolo della fattorizzazione usando la fattorizzazione di Cholesky anziché fattorizzazione LU. Teorema (Fattorizzazione di Cholesky) Sia A una matrice di ordine n simmetrica e definita positiva. Allora esiste un’unica matrice triangolare superiore R con elementi diagonali positivi tale che A = R T R. Costo computazionale per fattorizzazione di Cholesky: O( Stefano Berrone n3 ) 6 Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Problemi triangolari Eliminazione Gaussiana Fattorizzazione LU Pivoting Comandi MATLAB 1 [L,U,P]=lu(A) 2 R=chol(A) Attenzione! chol assume che A sia simmetrica definita positiva e usa solo la parte triangolare superiore! 3 \ Altre fattorizzazioni 1 Fattorizzazione QR (qr(A)): Q è ortogonale, R è triangolare superiore. Molto stabile ma anche più costosa della fattorizzazione PA = LU. Utile per altre applicazioni. Varie 1 norm(x), norm(A): calcolo di norme 2 cond(A): numero di condizionamento Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Metodi iterativi I metodi diretti consentono la costruzione della soluzione esatta, a meno degli errori di arrotondamento, in un numero finito di passi. Per sistemi con matrici dense i metodi diretti sono generalmente preferibili. Per sistemi con matrici sparse e di ordine elevato, i metodi diretti diventano impraticabili a causa del fill in ⇒ non può essere sfruttata la sparsità della matrice. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel I metodi iterativi lasciano inalterata la struttura di A e costruiscono una successione infinita di vettori {x (k) }k≥0 , che sotto opportune ipotesi converge alla soluzione x. Per i metodi iterativi dovremo discutere: 1 condizioni di applicabilità 2 costo computazionale (di ogni iterazione) 3 problemi connessi alla convergenza: convergenza si/no e velocità di convergenza. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Definizione Una successione di vettori {x (k) }k≥0 di Rn si dice convergente ad un vettore x ∈ Rn se esiste una norma per cui lim kx (k) − xk = 0, k→∞ ed in tal caso si pone lim x (k) = x. k→∞ Osservazione Grazie all’equivalenza delle norme su Rn tale definizione risulta indipendente dalla norma scelta in Rn ; se una successione di vettori converge in una norma di Rn , allora converge in tutte le norme di Rn . Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Osservazione La condizione di convergenza in una qualunque norma di Rn si traduce in una condizione di convergenza per componenti: lim x (k) = x k→∞ ⇐⇒ (k) lim x k→∞ i = xi , ∀i = 1, ..., n. Simili argomenti valgono anche per successioni di matrici {A(k) }k≥0 . Inoltre vale il seguente teorema: Teorema Sia A ∈ Rn×n . Allora indicando con ρ(A) il raggio spettrale di A lim Ak = 0 ∈ Rn×n ⇐⇒ ρ(A) < 1 k→∞ Ricordiamo che per ogni norma di matrice compatibile con una norma di vettore si ha ρ(A) ≤ kAk. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Generalità sui metodi iterativi Idea base: data una stima iniziale x (0) della soluzione del sistema lineare Ax = b, con A ∈ Rn×n , det(A) 6= 0, b ∈ Rn , si costruisce una successione di vettori {x (k) }k≥0 che converga (si spera) ad x, risolvendo ad ogni iterazione dei sistemi lineari molto più semplici di quello di partenza. Come? Consideriamo uno splitting della matrice A del tipo A = M + N, det(M) 6= 0. Allora si ha Ax = (M + N) x = b ⇐⇒ Mx = −Nx + b. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Supponiamo di essere all’iterazione k + 1, di avere quindi a disposizione x (k) e di voler calcolare x (k+1) . Definiamo x (k+1) in modo che valga Mx (k+1) = −Nx (k) + b Allora abbiamo x (k+1) = −M −1 Nx (k) + M −1 b che può essere scritto in forma compatta come x (k+1) = Bx (k) + c, avendo posto B = −M −1 N e c = M −1 b. La matrice B = −M −1 N è la matrice d’iterazione ed individua il particolare metodo; è fondamentale nello stabilire le condizioni e la velocità di convergenza. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Proprietà (consistenza) Se x (k) → x̃, allora x̃ soddisfa l’equazione x̃ = B x̃ + c e quindi il limite della successione x̃ è proprio la soluzione x di Mx = −Nx + b (cioè di Ax = b). Osservazione Al variare della stima iniziale x (0) si ottengono diverse successioni {x (k) }k≥0 . Definizione Un metodo iterativo si dice convergente se per ogni stima iniziale x (0) la successione {x (k) }k≥0 è convergente. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Studio della convergenza dei metodi iterativi Sotto quali condizioni la successione {x (k) }k≥0 generata da un metodo iterativo converge? Se la successione converge ad un vettore x, questo è senz’altro soluzione del sistema lineare. Definiamo l’errore e (k) commesso al passo k come la differenza fra l’approssimazione corrente e la soluzione: e (k) = x (k) − x. Chiaramente, x (k) → x ⇐⇒ e (k) → 0. Studiamo l’errore: e (k) = x (k) − x = Bx (k−1) + c − x = Bx (k−1) + c − (Bx + c) = B(x (k−1) − x) = Be (k−1) . Allora e (k) = Be (k−1) = B(Be (k−2) ) = B 2 e (k−2) = ... = B k e (1) = B k e (0) . Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Poiché e (k) = B k e (0) abbiamo lim x (k) = x ⇔ lim e (k) = 0 ⇔ lim B k e (0) = 0 ⇔ lim B k = 0 k→∞ k→∞ k→∞ k→∞ Quindi c’è convergenza se e solo se ρ(B) < 1. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Teorema Il metodo iterativo definito dalla matrice di iterazione B converge per ogni x (0) ∈ Rn se e solo se ρ(B) < 1. La condizione necessaria e sufficiente espressa dal teorema per la convergenza non è di facile verifica. Una condizione sufficiente è espressa dal seguente Teorema. Teorema Sia k · k una norma di matrice naturale. Se kBk < 1 il metodo iterativo converge. Dimostrazione: È sufficiente ricordare che per ogni norma naturale ρ(B) ≤ kBk Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Il raggio spettrale indica non solo se il metodo iterativo converge ma anche quanto velocemente converge. Intuitivamente, la velocitá di convergenza è una misura di quanto velocemente l’errore tende a zero. La velocitá di convergenza di un metodo iterativo è tanto maggiore quanto più il raggio spettrale della matrice di iterazione è vicino a 0; viceversa è tanto più lenta quanto più è vicino a 1. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Metodo di Jacobi Detto anche metodo delle sostituzioni simultanee. Consideriamo il seguente splitting di A: F D A = E + D + F, A= E dove E è la parte di A sotto la diagonale, 2 D è la diagonale di A, 3 F è la parte di A sopra la diagonale Prendiamo M = D, N = E + F . Matrice d’iterazione BJ : 1 BJ = −D −1 (E + F ) Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Data l’iterata corrente x (k) ∈ Rn , la x (k+1) ∈ Rn viene calcolata tramite il seguente algoritmo: for i = 1 : n P Pn (k+1) (k) (k) xi = bi − i−1 a x − a x /aii j=1 ij j j=i+1 ij j end In forma compatta, può essere scritto x (k+1) = D −1 b − Ex (k) − Fx (k) . Osservazione Se aii = 0 per qualche i, D non è invertibile. Si possono però preliminarmente riordinare le equazioni del sistema lineare in modo da portare in posizione ii un elemento non nullo. Ciò è sempre possibile se A è non singolare. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Metodo di Gauss-Seidel Consideriamo sempre lo splitting A = E + D + F, ma M = E + D e N = F . La matrice di iterazione: BGS = −(E + D)−1 F Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Data l’iterata corrente x (k) ∈ Rn , la x (k+1) ∈ Rn viene calcolata tramite il seguente algoritmo: for i = 1 : n P P (k+1) (k+1) (k) = bi − i−1 xi − nj=i+1 aij xj /aii j=1 aij xj end In forma compatta, può essere scritto x (k+1) = D −1 b − Ex (k+1) − Fx (k) . Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Test di arresto Un metodo iterativo fornisce una successione infinita di iterate x (k) . Questo diventerà un algoritmo solo se stabiliremo in modo preciso quando interrompere il procedimento di generazione di iterate. Dobbiamo ovvero stabilire un criterio di arresto delle iterazioni. 1 2 si arresta il processo iterativo quando risulta (k+1) x − x (k) ≤ tol, x (k+1) dove tol è una tolleranza relativa fissata sulla base di opportuni criteri di utilizzo della soluzione. oppure se (k+1) (k) − x ≤ tol, x dove tol è questa volta una tolleranza assoluta. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Una strada alternativa per definire un test d’arresto consiste nel controllo del residuo dell’equazione. Sia x (k+1) la soluzione ottenuta al passo k + 1, il residuo corrispondente è definito come r (k+1) = b − Ax (k+1) . Test di arresto relativo sul residuo: b − Ax (k+1) ≤ tol. kbk Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Osservazione Residuo ”piccolo” non garantisce errore ”piccolo”. (k) (k) (k) −1 −1 (k) e = x − x = x − A b = A Ax − b = −A−1 r (k) ≤ A−1 r (k) . Ricordando 1 kx k ≤ kAk kbk (v. conti su condizionamento) abbiamo (k) e (k) r −1 (k) k A k ≤ A = K (A) r kx k kbk kbk Quindi (k) r kbk ≤ tol =⇒ (k) e kx k Stefano Berrone ≤ K (A) (k) r kbk ≤ K (A)tol Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Teorema (Condizioni di convergenza per il metodo di Jacobi) Per matrici a predominanza diagonale stretta (per righe o per colonne indifferentemente) il metodo di Jacobi converge. Teorema (Condizioni di convergenza per il metodo di Gauss-Seidel) Per matrici a predominanza diagonale stretta (per righe o per colonne) il metodo di Gauss-Seidel converge. Teorema (Condizioni di convergenza per il metodo di Gauss-Seidel) Per matrici simmetriche definite positive il metodo di Gauss-Seidel converge. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Confronto delle proprietà di convergenza del metodo di Jacobi e di Gauss-Sidel In generale la convergenza di uno dei metodi non implica la convergenza dell’altro e viceversa. Tuttavia quando entrambi convergono la velocità di convergenza del metodo di Gauss-Seidel è generalmente superiore. Teorema (Teorema di Stein-Rosenberg) Sia A ∈ Rn×n con aij ≤ 0, ∀i 6= j e aii > 0, allora si verifica uno e uno solo dei seguenti risultati: 1 0 < ρ(BGS ) < ρ(BJ ) < 1, 2 1 < ρ(BJ ) < ρ(BGS ), 3 ρ(BGS ) = ρ(BJ ) = 0, 4 ρ(BGS ) = ρ(BJ ) = 1. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari Metodi diretti Metodi Iterativi Generalità Metodo di Jacobi Metodo di Gauss-Seidel Test di arresto Condizioni di convergenza per Jacobi e Gauss-Seidel Teorema Sia A ∈ Rn×n una matrice tridiagonale con elementi diagonali non nulli, allora ρ(BGS ) = ρ2 (BJ ), cioè il metodo di Gauss-Seidel ed il metodo di Jacobi convergono o divergono simultaneamente e il tasso asintotico di convergenza del metodo di Gauss-Seidel è doppio di quello del metodo di Jacobi. Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari