Algebra lineare

G. Puppo

Algebra lineare

Numero di condizionamentoSistemi triangolariFattorizzazione LUSoluzione di sistemi lineari in MatlabEffetto del numero di

condizionamento

Numero di condizionamento

Per calcolare il numero di condizionamento, Matlab usa l’istruzione cond(a,p), dove:

p = 1 è per la norma 1p = 2 è per la norma 2 (default)p = inf è per la norma infinito

Andamento del numero di condizionamento

Costruiamo un programma che calcoli il numero di condizionamento al variare di N per matrici N per N assegnate.In particolare il programma deve:avere un ciclo su Ncalcolare le matrici richieste per ogni Ncalcolare il numero di condizionamento

corrispondente per ogni matrice, in ognuna delle norme 1, 2 e inf

stampare un grafico con i risultati

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 inffor n=1:nmax a=rand(n); con1(n) = cond(a,1); con2(n) = cond(a,2); coninf(n) = cond(a,inf);endplot(log10(con2))hold onplot(log10(con1),'g')plot(log10(coninf),'m')(log10(

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).

Matrice tridiagonale tridiag(n)

Matrici di Hilbert

La scala verticale è logaritmica

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

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') returnendfor i=1:n for j=(i+1):n if a(i,j) ~=0 display('A non e'' triangolare inf') return end endend

Termino il controllo dei dati

% Controlla le dimensioni di Bif length(b) ~= n display('B non e'' compatibile') returnend% Controlla se A e' singolare:for i=1:n if abs( a(i,i)) < eps*norm(b) display('A e'' quasi singolare') return endend

Finalmente, risolvo il sistema:

%Risolve il sistemax(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

La formula ricorsiva è:

L’algoritmo corrispondente è:

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') returnendfor i=1:n for j=1:i-1 if a(i,j) ~=0 display('A non e'' triangolare sup') return end endend

Controllo che B sia compatibile e che A non sia singolare

% Controlla le dimensioni di Bif length(b) ~= n display('B non e'' compatibile') returnend% Controlla se A e' singolare:for i=1:n if abs( a(i,i)) < eps*norm(b) display('A e'' quasi singolare') return endend

Risolvo il sistema:

%Risolve il sistemax(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

La formula ricorsiva è:

E l’algoritmo corrispondente è:

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') returnend

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

…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); endend

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:

>> [l,u]=elleu(a);>> y=tri_inf(l,b);>> x=tri_sup(u,y)

Calcolare la fattorizzazione LU;Risolvere il sistema triangolare inferiore Ly=b;Risolvere il sistema triangolare superiore Ux=y;

Esempio

Risolvo il sistema Ax=b, per:

>> a=[2 -1 0; -1 2 -1; 0 -1 2];>> aa = 2 -1 0 -1 2 -1 0 -1 2>> b=ones(3,1);

Calcola la fattorizzazione LU:

>> [l,u]=elleu(a)l = 1.0000 0 0 -0.5000 1.0000 0 0 -0.6667 1.0000u = 2.0000 -1.0000 0 0 1.5000 -1.0000 0 0 1.3333

Osservo che: >> l*uans = 2 -1 0 -1 2 -1 0 -1 2

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'-bans = 1.0e-015 * 0 0.2220 -0.4441

Quindi, il residuo è piccolo, anche se non è zero. La soluzione è accettabile

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

L’effetto del numero di condizionamento si vede meglio con questoesempio:

>> a=hilb(5);>> x=ones(5,1);>> b=a*x;>> [l,u]=elleu(a);>> y=tri_inf(l,b);>> x1=tri_sup(u,y);

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

La differenza ||x-x1|| è molto più grande del residuo

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 lasoluzione 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 diwarning, ma calcola ugualmente la soluzione.

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)

Esempio

>> a=[2 -1 0; -1 2 -1; 0 -1 2];>> b=ones(3,1);>> x=a\bx = 1.5000 2.0000 1.5000

Attenzione!Matlab usa sia / che \ ma i due operatori hanno un effetto diverso.Provare per credere!

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.

Fattorizzazione LU

Matlab calcola la fattorizzazione LU di una matrice con il comando:

[l,u] = lu(a)

Oppure posso usare tre argomenti in output::

[l,u,p]=lu(a)

In questo caso, l=L, u=U, p=P.

In questo caso, l contiene anche gli scambi di riga effettuati.Se la fattorizzazione è PA = LU, allora l=P-1 L.

Esempio

>> 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 0u = 7.0000 8.0000 9.0000 0 0.8571 1.7143 0 0 0.0000

Notare che l contiene lamatrice triangolare Lcon le righe scambiate,cioè l=P-1 L

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;

Fattorizzazione QR

Matlab esegue la fattorizzazione A=QR mediante il comando:

>> [q,r]=qr(a);>> x=r\q'*b;

>>[q,r] = qr(a)

Per forzare la soluzione di un sistema lineare mediante fattorizzazione QR devo quindi dare i comandi:

Esercizio

Effetto del numero di condizionamento

Variazione del residuo per matrici mal condizionate

Effetto del condizionamento sulla precisione della soluzione numerica

Stima del rango

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

%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)endplot(log10(c))hold onplot(log10(r+eps),'m')

Script res_hilb.m

Ottengo questi risultati

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.

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 condizionamentonmax=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);endplot(log10(c))hold onplot(log10(err+eps),'m')

Ottengo questi risultati

Stima del rango di una matricePer stimare il rango di una matrice posso:

Per matrici mal condizionate il secondo metodo è più affidabile.La tolleranza assegnata ha un ruolo fondamentale

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

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

Questa possibilità di scelta è implementata con la functionNARGIN, 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 RANGOcalcola la tolleranza, TOL=EPS*NORM(A,1)altrimenti TOL non viene calcolata, ma viene usato il valoreassegnato in input

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)=0if 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; endend

Provo ora a calcolare il rango delle matrici di Hilbert.Uso le seguenti istruzioni:

for i=1:100r(i)=rango(hilb(i));end >> plot(r)

Notare che tutte le matrici di Hilbert sono non singolari,quindi dovrei avere r(i)=i

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

Rango delle matrici di Hilbert, con TOL=1e-40

Ora tutte le matrici hanno rango pieno

Apparentemente, i risultati ottenuti con TOL=1e-40 sonopiù corretti di quelli ottenuti con TOL=EPS*NORM(A,1).Calcoliamo però qual è l’accuratezza della fattorizzazione QRche 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 dicirca 1e-16, perché la matrice R non è abbastanza accurata