Risoluzione di sistemi lineari Lucia Gastaldi Dipartimento di Matematica, http://dm.ing.unibs.it/gastaldi/ Laboratorio 30 Gennaio 2006 Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Outline 1 Matrici particolari Matrici predefinite Matrici strutturate 2 Risoluzione di sistemi lineari Il comando per risolvere i sistemi lineari Fattorizzazione 3 Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Matrici predefinite Matrici strutturate Matrici predefinite ones A=ones(m,n) produce la matrice A di dimensioni mxn, i cui elementi sono uguali ad 1. zeros B=zeros(m,n) produce la matrice B di dimensioni mxn, i cui elementi sono uguali ad 0. eye I=eye(N) produce la matrice identità I di dimensioni NxN. C=eye(m,n) produce la matrice C di dimensioni mxn, che ha 1 sulla diagonale principale e 0 fuori. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Matrici predefinite Matrici strutturate Matrici diagonali e diagonali di una matrice diag - da vettori a matrici Sia V un vettore di N componenti. diag(V,K) è una matrice quadrata di ordine N+ABS(K) che ha gli elementi di V sulla diagonale K-esima. se K= 0 se K> 0 se K< 0 è la diagonale principale si trova sopra la diagonale principale si trova sotto la diagonale principale diag(V)è uguale a diag(V,0) e pone V sulla diagonale principale. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Matrici predefinite Matrici strutturate Matrici diagonali e diagonali di una matrice diag - estrazione di una diagonale da una matrice Sia X una matrice. diag(X,K) è un vettore colonna formato dagli elementi della diagonale K-esima di X. diag(X) è la diagonale principale di X. diag(diag(X)) è la matrice diagonale che ha la stessa diagonale principale di X. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Matrici predefinite Matrici strutturate Esercizio Costruire il vettore colonna V di dimensione 5 che ha tutti gli elementi pari a 6. Costruire la matrice A che ha V sulla seconda diagonale sopra la diagonale principale. Sommare alla matrice A la matrice identità. Creare la matrice tridiagonale di dimensione 8 che ha sulla diagonale principale i numeri interi da 1 a 8, sulla diagonale superiore tutti gli elementi pari a -2 e sulla diagonale inferiore tutti gli elementi pari a -4. spy Il comando spy permette di visualizzare la struttura di una matrice. Con un puntino blu sono contrassegnati gli elementi diversi da zero. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Matrici predefinite Matrici strutturate Matrici triangolari tril L=tril(A) è la parte triangolare inferiore della matrice A. L=tril(A,k) è una matrice che contiene gli elementi di A che si trovano sulla diagonale k-esima e sotto di essa. Sopra la diagonale k-esima gli elemnti sono tutti nulli. triu U=triu(A) è la parte triangolare superiore della matrice A. U=triu(A,k) è una matrice che contiene gli elementi di A che si trovano sulla diagonale k-esima e sopra di essa. Sotto la diagonale k-esima gli elemnti sono tutti nulli. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Sistemi lineari Problema Data la matrice A ∈ Rm×n , ed il vettore b ∈ Rm , trovare il vettore x ∈ Rn tale che Ax = b. Tre casi possibili: Sistemi quadrati, m = n. Sistemi sovradeterminati, m > n. Sistemi sottodeterminati, m < n. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Come risolvere un sistema lineare con MATLAB La risoluzione del sistema lineare si ottiene usando il simbolo di divisione: backslash \. x = A\b indica la soluzione del sistema Ax = b essendo b un vettore colonna di m componenti. L’operatore backslash usa algoritmi differenti per trattare diversi tipi di matrici: Permutazioni di matrici triangolari. Matrici simmetriche e definite positive. Matrici quadrate, non singolari m = n. Sistemi rettangolari sovradeterminati m > n. Sistemi rettangolari sottodeterminati m < n. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Metodo di sostituzione all’indietro Sia U ∈ Rn×n una matrice triangolare superiore; se uii 6= 0 per i = 1, . . . , n, la soluzione del sistema lineare Ux = b è data da: xn = bn unn n X 1 bi − uij xj xi = uii per i = n − 1, . . . , 1. j=i+1 Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Metodo di sostituzione in avanti Sia L ∈ Rn×n una matrice triangolare inferiore; se `ii 6= 0 per i = 1, . . . , n, la soluzione del sistema lineare Lx = b è data da: x1 = b1 `11 i−1 X 1 bi − `ij xj xi = `ii per i = 2, . . . , n. j=1 Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Esercizi Esercizio: metodo di sostituzione in avanti Scrivere una function function x=sost avanti(L,b) che risolva un sistema lineare con matrice triangolare inferiore L e termine noto b. Esercizio: metodo di sostituzione all’indietro Scrivere una function function x=sost indietro(U,b) che risolva un sistema lineare con matrice triangolare superiore U e termine noto b. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Algoritmo di eliminazione di Gauss for k = 1, . . . , n − 1 for i = k + 1, . . . , n (k) aik mik = (k) akk for j = k + 1, . . . , n (k+1) (k) (k) aij = aij − mik akj end (k+1) (k) (k) bi = bi − mik bk end end Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Fattorizzazione LU Teorema Costruiamo le seguenti matrici: 1 0 m21 1 L= ··· ··· mn1 · · · (1) (1) a a12 · · · ··· 0 11 (2) ··· 0 U = 0 a22 · · · ··· ··· ··· ··· ··· mnn−1 1 0 ··· 0 Allora LU = A. La fattorizzazione senza pivoting è realizzata dalla function lu gauss con la seguente sintassi: [L,U]=lu gauss(A) Lucia Gastaldi Sistemi lineari (1) a1n (2) a2n ··· (n) ann . Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Fattorizzazione LU senza pivoting Esercizio Scrivere un file di tipo script che raccoglie i seguenti comandi essendo memorizzati nel workspace una matrice A ed un vettore noto b: calcola la fattorizzazione LU di A con la function lu gauss; risolve il sistema triangolare inferiore Ly=b richiamando la propria function sost avanti; risolve il sistema triangolare superiore Ux=y richiamando la propria function sost indietro. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Fattorizzazione LU senza pivoting Esercizio Usare lo script per 1 2 A= 0 1 risolvere il sistema lineare Ax = b con: −2 −1 0 0 1 8+α 1 9+α b= 1 2 1 1 0 1 0 1 ed α = −1, 0, 10−15 . Calcolare l’errore relativo rispetto alla soluzione x = (1, 0, 1, −1)T . Calcolare anche la matrice B = A − LU. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Il comando per risolvere i sistemi lineari Fattorizzazione Uso della function lu Data la matrice A ∈ Rn×n , la function lu fornisce il risultato della fattorizzazione nelle seguenti forme: [L,U]=lu(A): fornisce le matrici L e U in modo che LU=A. [L,U,P]=lu(A): fornisce le matrici L, U e P in modo che LU=PA. Verificare il comportamento di lu sulla matrice: 1 2 −1 1 . A= 3 2 2 −1 0 Usare la fattorizzazione LU con pivoting per calcolare la fattorizzazione della matrice dell’esercizio precedente per i diversi valori di α. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi Numero di condizionamento Definizione K (A) = kA−1 k kAk si dice numero di condizionamento della matrice A. Teorema Si consideri il sistema lineare Ax = b. Siano δA e δb perturbazioni di A e di b rispettivamente e sia x + δx la soluzion del sistema lineare: (A + δA)(x + δx) = b + δb. Allora vale la seguente maggiorazione: kδxk K (A) ≤ kxk 1 − K (A)kδAk/kAk Lucia Gastaldi kδAk kδbk + kAk kbk Sistemi lineari . Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi Le funzioni MATLAB per il calcolo del numero di condizionamento Funzione cond rcond condest norm Significato numero di condizionamento stima del reciproco del numero di condizionamento in norma 1 stima del numero di condizionamento in norma 1 norma di vettore o di matrice Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi norm Se A è un array (matrice o vettore), il comando norm(A) calcola la norma 2 di A. Il comando norm(A,p) calcola la norma p di A come segue: Valore 1 2 ’fro’ inf Norma norma 1 norma 2 norma di Frobenius norma ∞ Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi cond Data una matrice A ∈ Rn×n il suo numero di condizionamento in norma 2 si calcola con il comando: c=cond(A) Per calcolare il numero di condizionamento in norma 1, ∞ o di Frobenius si usa il comando c=cond(A,p) con Valore 1 2 ’fro’ inf Norma norma 1 norma 2 norma di Frobenius norma ∞ Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi rcond Il comando c=rcond(A) calcola una stima del reciproco del numero di condizionamento di A in norma 1. Se rcond(A) è vicino a 1, allora la matrice è ben condizionata. Se rcond(A) è vicino a 0, allora la matrice è mal condizionata. Confrontato con cond, rcond è più efficiente ma meno affidabile. rcond rcond dà una misura della distanza della matrice dall’insieme delle matrici singolari. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi Numero di condizionamento Dati 1 2 −1 2 4.0000001 0 5 2 , A= −2 −4.0000001 1 −6 0 0.0000001 0 3.1 4 11.0000001 , b= −11.0000001 3.1000001 calcolare la soluzione esatta del sistema Ax = b. Introdotta le perturbazioni r1 = 10−8 [−1, −1, −1, −1]T ed r2 = 10−8 [−1, 0, 1, 0]T al termine noto, si considerino le soluzioni del sistema Ax̂1 = b + r1 e Ax̂2 = b + r2 . Calcolare l’errore relativo commesso nei due casi, e confrontarlo con la perturbazione relativa del termine noto. Calcolare il numero di condizionamento della matrice e verificare se la stima teorica è soddisfatta. Lucia Gastaldi Sistemi lineari Matrici particolari Risoluzione di sistemi lineari Analisi degli errori Numero di condizionamento Calcolo del numero di condizionamento con Matlab Esercizi Esempio di matrice mal condizionata Si consideri il sistema lineare Ax = b con A ∈ Rn×n matrice di Hilbert di elementi aij = 1 , i +j −1 i = 1, . . . , n, e b ∈ Rn tale che la soluzione del sistema sia x = (1, . . . , 1)T . Per costruire la matrice di Hilbert usare il comando hilb(n). Calcolare con MATLAB la fattorizzazione LU con pivoting e risolvere il sistema al variare di n. Sia x̂ la soluzione calcolata. Riportare l’errore relativo E = kx − x̂k/kxk in un grafico in scala semilogaritmica per diversi valori di n. Riportare sullo stesso grafico anche la norma della matrice R = LU − PA e del residuo kb − Ax̂k/kbk. Per calcolare le norme usare il comando norm. Lucia Gastaldi Sistemi lineari
Scarica
