1 / 49

Algebra lineare

Algebra lineare. G. Puppo. Algebra lineare. Numero di condizionamento Sistemi triangolari Fattorizzazione LU Soluzione di sistemi lineari in Matlab Effetto del numero di condizionamento. Numero di condizionamento.

oshin
Download Presentation

Algebra lineare

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Algebra lineare G. Puppo

  2. Algebra lineare • Numero di condizionamento • Sistemi triangolari • Fattorizzazione LU • Soluzione di sistemi lineari in Matlab • Effetto del numero di condizionamento

  3. Numero di condizionamento • Per calcolare il numero di condizionamento, Matlab usa l’istruzione cond(a,p), dove: • p = 1 è per la norma 1 • p = 2 è per la norma 2 (default) • p = inf è per la norma infinito

  4. Andamento del numero di condizionamento • avere un ciclo su N • calcolare le matrici richieste per ogni N • calcolare il numero di condizionamento corrispondente per ogni matrice, in ognuna delle norme 1, 2 e inf • stampare un grafico con i risultati Costruiamo un programma che calcoli il numero di condizionamento al variare di N per matrici N per N assegnate. In particolare il programma deve:

  5. function plot_con(nmax) %Calcola il numero di condizionamento per le matrici %RAND(N), con N=1:NMAX, confrontando su un grafico %i risultati ottenuti in norma 1, norma 2, e norma inf for n=1:nmax a=rand(n); con1(n) = cond(a,1); con2(n) = cond(a,2); coninf(n) = cond(a,inf); end plot(log10(con2)) hold on plot(log10(con1),'g') plot(log10(coninf),'m')(log10(

  6. Proviamo altre matrici ... • Cambiamo le matrici calcolate nella function plot_con • Sostituiamo a rand(n) le matrici tridiagonali, fornite dalla function tridiag(n) • Sostituiamo le matrici di Hilbert, hilb(n), che hanno un condizionamento ~exp(n).

  7. Matrice tridiagonale tridiag(n)

  8. Matrici di Hilbert La scala verticale è logaritmica

  9. Sistemi triangolari inferiori Scrivo una function che calcola la soluzione del sistema Ax=b, nel caso in cui A è triangolare inferiore. In una function che deve essere utilizzata da diversi utenti è bene controllare che i dati forniti siano coerenti. In questo caso controllo che: - A sia quadrata - A sia triangolare inferiore - A sia non singolare - Il vettore b abbia le dimensioni giuste

  10. Controllo che la matrice A sia ben impostata function x=tri_inf(a,b) %TRI_INF(A,B) Risolve il sistema triangolare inferiore AX=B. % Se A non e' quadrata, o A non e' triangolare inferiore, % stampa un messaggio di errore % Controlla le dimensioni di A: [n,m] = size(a); if m ~= n display('A non e'' quadrata') return end for i=1:n for j=(i+1):n if a(i,j) ~=0 display('A non e'' triangolare inf') return end end end

  11. Termino il controllo dei dati % Controlla le dimensioni di B if length(b) ~= n display('B non e'' compatibile') return end % Controlla se A e' singolare: for i=1:n if abs( a(i,i)) < eps*norm(b) display('A e'' quasi singolare') return end end

  12. Finalmente, risolvo il sistema: La formula ricorsiva è: %Risolve il sistema x(1) = b(1)/a(1,1); for i=2:n sum=b(i); for j=1:(i-1) sum = sum - a(i,j)*x(j); end x(i) = sum/a(i,i); end L’algoritmo corrispondente è:

  13. Risolvo un sistema triangolare superiore. Controllo che A sia ben impostata: function x=tri_sup(a,b) %TRI_SUP(A,B) Risolve il sistema triangolare superiore AX=B. % Se A non e' quadrata, o A non e' triangolare superiore, % stampa un messaggio di errore % Controlla le dimensioni di A: [n,m] = size(a); if m ~= n display('A non e'' quadrata') return end for i=1:n for j=1:i-1 if a(i,j) ~=0 display('A non e'' triangolare sup') return end end end

  14. Controllo che B sia compatibile e che A non sia singolare % Controlla le dimensioni di B if length(b) ~= n display('B non e'' compatibile') return end % Controlla se A e' singolare: for i=1:n if abs( a(i,i)) < eps*norm(b) display('A e'' quasi singolare') return end end

  15. Risolvo il sistema: La formula ricorsiva è: %Risolve il sistema x(n) = b(n)/a(n,n); for i=n-1:-1:1 sum=b(i); for j=i+1:n sum = sum - a(i,j)*x(j); end x(i) = sum/a(i,i); end E l’algoritmo corrispondente è:

  16. Calcola la fattorizzazione LU di una matrice function [l,u] = elleu(a) %ELLEU(A) Calcola la fattorizzazione LU di A, senza pivoting %Sintassi: [L,U]=ELLEU(A) % esce con un messaggio di errore, se trova un pivot < % EPS*NORM(A) nor=norm(a); % Controlla le dimensioni di A: [n,m] = size(a); if m ~= n display('A non e'' quadrata') return end

  17. Calcola la fattorizzazione LU, riscrivendo su A: l=-eye(n); for k=1:n-1 if abs(a(k,k)) < eps*nor display('Pivot troppo piccolo') return end m(k+1:n,k) = -a(k+1:n,k)/a(k,k); a(k+1:n,k+1:n) = a(k+1:n,k+1:n)+m(k+1:n,k)*a(k,k+1:n); end

  18. …infine, immagazzina i risultati nelle matrici L e U: % Immagazzina i risultati nelle matrici l e u: l = -m; u = zeros(n); for i=1:n for j=i:n u(i,j) = a(i,j); end end

  19. Risoluzione di un sistema lineare Con le tre functions appena costruite, posso risolvere un sistema lineare. Devo impostare una matrice quadrata A e un vettore B di termini noti. Per risolvere il sistema Ax=B, devo: • Calcolare la fattorizzazione LU; • Risolvere il sistema triangolare inferiore Ly=b; • Risolvere il sistema triangolare superiore Ux=y; >> [l,u]=elleu(a); >> y=tri_inf(l,b); >> x=tri_sup(u,y)

  20. Esempio Risolvo il sistema Ax=b, per: >> a=[2 -1 0; -1 2 -1; 0 -1 2]; >> a a = 2 -1 0 -1 2 -1 0 -1 2 >> b=ones(3,1);

  21. Calcola la fattorizzazione LU: >> [l,u]=elleu(a) l = 1.0000 0 0 -0.5000 1.0000 0 0 -0.6667 1.0000 u = 2.0000 -1.0000 0 0 1.5000 -1.0000 0 0 1.3333 Osservo che: >> l*u ans = 2 -1 0 -1 2 -1 0 -1 2

  22. Risolvo i due sistemi triangolari >> y=tri_inf(l,b) y = 1.0000 1.5000 2.0000 >> x=tri_sup(u,y) x = 1.5000 2.0000 1.5000 Osservo che >> a*x'-b ans = 1.0e-015 * 0 0.2220 -0.4441 Quindi, il residuo è piccolo, anche se non è zero. La soluzione è accettabile

  23. Provo un altro esempio con a=hilb(5), b=ones(5,1). Ottengo: >> a=hilb(5); >> b=ones(5,1); >> [l,u]=elleu(a); >> y=tri_inf(l,b); >> x=tri_sup(u,y); >> norm(a*x'-b) ans = 3.1776e-014 Quindi il residuo è aumentato, infatti: >> cond(a) ans = 4.7661e+005

  24. L’effetto del numero di condizionamento si vede meglio con questo esempio: Scelgo b, in modo da conoscere la soluzione esatta, x. Calcolo la soluzione x1 risolvendo il sistema lineare. In aritmetica esatta, avrei x-x1=0. In aritmetica floating point: >> norm(x-x1')/norm(x) ans = 1.5531e-012 >> a=hilb(5); >> x=ones(5,1); >> b=a*x; >> [l,u]=elleu(a); >> y=tri_inf(l,b); >> x1=tri_sup(u,y); La differenza ||x-x1|| è molto più grande del residuo

  25. Risoluzione di sistemi lineari con Matlab Matlab risolve il sistema lineare Ax=b con il comando: x=A\b. Se A è quadrata, questo comando implica i seguenti passi: - Calcola la fattorizzazione LU con pivoting sulle colonne - Risolvi i due sistemi triangolari Se A è rettangolare (sistemi sovradeterminati), Matlab calcola la soluzione nel senso dei minimi quadrati, cioè: - Calcola la fattorizzazione QR con pivoting sulle colonne - Risolve il sistema triangolare superiore Se il condizionamento di A è elevato, stampa un messaggio di warning, ma calcola ugualmente la soluzione.

  26. Risoluzione di sistemi lineari con Matlab 2 Nella versione 7, è stata inserita una nuova function per risolvere sistemi lineari con struttura particolare x=linsolve(a,b)

  27. Esempio >> a=[2 -1 0; -1 2 -1; 0 -1 2]; >> b=ones(3,1); >> x=a\b x = 1.5000 2.0000 1.5000 Attenzione! Matlab usa sia / che\ ma i due operatori hanno un effetto diverso. Provare per credere!

  28. Attenzione! Matlab calcola una soluzione anche quando il sistema non ammette soluzione. Nell’esempio seguente, a è singolare: >> a=[1 2 3; 4 5 6; 7 8 9]; >> b=[1; 1; 0]; >> x=a\b; Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-018. >> norm(a*x-b) ans = 5.0990 RCOND è il reciproco del numero di condizionamento: se RCOND è piccolo, la matrice è mal condizionata.

  29. Fattorizzazione LU Matlab calcola la fattorizzazione LU di una matrice con il comando: [l,u] = lu(a) In questo caso, l contiene anche gli scambi di riga effettuati. Se la fattorizzazione è PA = LU, allora l=P-1 L. Oppure posso usare tre argomenti in output:: [l,u,p]=lu(a) In questo caso, l=L, u=U, p=P.

  30. Esempio Notare che l contiene la matrice triangolare L con le righe scambiate, cioè l=P-1 L >> a=[1 2 3; 4 5 6; 7 8 9]; >> [l,u]=lu(a) l = 0.1429 1.0000 0 0.5714 0.5000 1.0000 1.0000 0 0 u = 7.0000 8.0000 9.0000 0 0.8571 1.7143 0 0 0.0000

  31. Per forzare la soluzione di un sistema tramite fattorizzazione LU, devo quindi dare i comandi: >> [l,u]=lu(a); >> y=l\b; >> x=u\y;

  32. Fattorizzazione QR Matlab esegue la fattorizzazione A=QR mediante il comando: >>[q,r] = qr(a) Per forzare la soluzione di un sistema lineare mediante fattorizzazione QR devo quindi dare i comandi: >> [q,r]=qr(a); >> x=r\q'*b;

  33. Esercizio

  34. Effetto del numero di condizionamento • Variazione del residuo per matrici mal condizionate • Effetto del condizionamento sulla precisione della soluzione numerica • Stima del rango

  35. Variazione del residuo per matrici mal condizionate Scrivo un programma che calcoli, per n da 1 a 100: - a=hilb(n) - b=rand(n,1) - risolve il sistema a*x=b - calcola la norma del residuo a*x-b - calcola il condizionamento di a - stampa un grafico con il condizionamento ed il residuo in funzione di n Vorrei studiare come varia il residuo r = ||Ax-b||, in funzione del condizionamento della matrice A. Scelgo come matrice mal condizionata la matrice di Hilbert, mentre per b assegno un vettore di numeri casuali

  36. Script res_hilb.m %Esegue un grafico del residuo, in funzione di n, %per sistemi lineari del tipo Hilb(n)*x=rand(n,1) nmax=100; for n=1:nmax a=hilb(n); b=rand(n,1); x=a\b; c(n)=cond(a); r(n)=norm(a*x-b) end plot(log10(c)) hold on plot(log10(r+eps),'m')

  37. Ottengo questi risultati

  38. Effetto del condizionamento sulla precisione della soluzione Scrivo un programma che calcoli, per n da 1 a 100: - a=hilb(n) - xexa=ones(n,1) (xexa=soluzione esatta) - b=a*xexa (b = termine noto corrispondente a xexa) - risolve il sistema a*x=b (x=soluzione approssimata) - calcola l’errore relativo fra x e xexa - calcola il condizionamento di a - stampa un grafico con il condizionamento e l’errore relativo Imposto un problema Ax=b di cui conosco la soluzione esatta. Risolvo numericamente il problema, calcolando la soluzione stimata y. Quindi calcolo l’errore ||x-y||. Mi aspetto errori grandi se uso per A una matrice mal condizionata, come, di nuovo, la matrice di Hilbert.

  39. Script hilb_exa.m % Calcola la differenza fra soluzione esatta e soluzione % approssimata, per i sistemi lineari hilb(n)*x=b % e ne stampa il grafico, confrontando con il condizionamento nmax=100; for n=1:nmax a=hilb(n); xexa=ones(n,1); b=a*xexa; x=a\b; c(n)=cond(a); err(n)=norm(x-xexa)/norm(xexa); end plot(log10(c)) hold on plot(log10(err+eps),'m')

  40. Ottengo questi risultati

  41. Stima del rango di una matrice Per stimare il rango di una matrice posso: • Calcolare la fattorizzazione LU e contare gli elementi U(i,i) che hanno modulo maggiore di una tolleranza assegnata • Calcolare la fattorizzazione QR e contare gli elementi R(i,i) che hanno modulo maggiore di una tolleranza assegnata Per matrici mal condizionate il secondo metodo è più affidabile. La tolleranza assegnata ha un ruolo fondamentale

  42. Listato della function rango.m function r=rango(a,tol) % Calcola il rango della matrice A, usando % il metodo QR % Sintassi: R=RANGO(A) Calcola il rango di A. Usa EPS % se |r(i,i)| <EPS*norm(a,1) => r(i,i)=0 % R=RANGO(A,TOL) Calcola il rango di A, % se |r(i,i)| <TOL => r(i,i)=0 Questa function può essere chiamata in due modi: rango(a) usa una tolleranza di default: TOL=EPS*NORM(a,1) rango(a,tol) usa la tolleranza TOL fissata in input

  43. Questa possibilità di scelta è implementata con la function NARGIN, che conta il numero di argomenti in input: if nargin < 2 tol=eps*norm(a,1); end Quindi, se NARGIN ritorna il valore 1, la function RANGO calcola la tolleranza, TOL=EPS*NORM(A,1) altrimenti TOL non viene calcolata, ma viene usato il valore assegnato in input

  44. function r=rango(a,tol) % Calcola il rango della matrice A, usando % il metodo QR % Sintassi: R=RANGO(A) Calcola il rango di A. Usa EPS % se |r(i,i)| <EPS*norm(a,1) => r(i,i)=0 % R=RANGO(A,TOL) Calcola il rango di A, % se |r(i,i)| <TOL => r(i,i)=0 if nargin < 2 tol=eps*norm(a,1); end [q,rr]=qr(a); r=0; [m,n]=size(a); for i=1:n if abs(rr(i,i))>tol r=r+1; end end

  45. Provo ora a calcolare il rango delle matrici di Hilbert. Uso le seguenti istruzioni: for i=1:100 r(i)=rango(hilb(i)); end >> plot(r) Notare che tutte le matrici di Hilbert sono non singolari, quindi dovrei avere r(i)=i

  46. Rango delle matrici di Hilbert, con TOL=EPS*NORM(A,1)

  47. Rango delle matrici di Hilbert, con TOL=1e-40 Ora tutte le matrici hanno rango pieno

  48. Apparentemente, i risultati ottenuti con TOL=1e-40 sono più corretti di quelli ottenuti con TOL=EPS*NORM(A,1). Calcoliamo però qual è l’accuratezza della fattorizzazione QR che stiamo calcolando, per esempio per HILB(50): >> a=hilb(50); >> [q,r]=qr(a); >> norm(q*r-a) ans = 5.7905e-016 Quindi non ha senso usare un valore di TOL più basso di circa 1e-16, perché la matrice R non è abbastanza accurata

More Related