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
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
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 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
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(
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')
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
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
Finalmente, risolvo il sistema:
La formula ricorsiva è:
L’algoritmo corrispondente è:
%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
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
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
Risolvo il sistema:
La formula ricorsiva è:
E l’algoritmo corrispondente è:
%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
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
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);
end
end
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)
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);
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
-1 2
0 -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'-b
ans =
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 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:
>> a=hilb(5);
>> x=ones(5,1);
>> b=a*x;
>> [l,u]=elleu(a);
>> y=tri_inf(l,b);
>> x1=tri_sup(u,y);
>> 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 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.
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\b
x=
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)
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.
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
0
u=
7.0000 8.0000 9.0000
0 0.8571 1.7143
0
0 0.0000
Notare che l contiene la
matrice triangolare L
con 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)
Per forzare la soluzione di un sistema lineare mediante
fattorizzazione QR devo quindi dare i comandi:
>> [q,r]=qr(a);
>> x=r\q'*b;
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:
-Vorrei
a=hilb(n)
studiare come varia il residuo
-r b=rand(n,1)
= ||Ax-b||, in funzione del condizionamento
- risolve il sistema a*x=b
della matrice A. Scelgo come matrice mal
- calcola la norma del residuo a*x-b
condizionata
la
matrice
di
Hilbert,
mentre
per
- calcola il condizionamento di a
b- stampa
assegno
vettore
numeri casualied il residuo
un un
grafico
con ildicondizionamento
in funzione di n
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')
Ottengo questi risultati
Effetto del condizionamento
sulla precisione della
soluzione
Scrivo un programma che calcoli, per n da 1 a 100:
- a=hilb(n)
Imposto
un problema Ax=b di cui conosco la
- xexa=ones(n,1) (xexa=soluzione esatta)
soluzione
esatta.
numericamente
il
- b=a*xexa
(b =Risolvo
termine noto
corrispondente a xexa)
- risolve il calcolando
sistema a*x=b la
(x=soluzione
approssimata)
problema,
soluzione
stimata y.
- calcola
l’errore l’errore
relativo fra||x-y||.
x e xexa
Quindi
calcolo
- calcola il condizionamento di a
Mi -aspetto
se uso per eAl’errore
una
stampa unerrori
grafico grandi
con il condizionamento
relativo
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 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')
Ottengo questi risultati
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
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 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
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
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
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 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