Fattorizzazione LU con pivoting parziale

La possibilità che $\exists i \vert a_{i\:i}^i =0\:i=1 \ldots n$ non dà noia più di tanto, inquanto siamo sicuri che una colonna non può essere nulla ( $\det(A)\neq 0$) e quindi con un piccolo accorgimento è possibile applicare la fattorizzazione LU ad ogni matrice non singolare. Infatti al passo i-esimo basta permutare l'eventuale elemento nullo con un altro non nullo e riflettere questa permutazione su tutta la riga (anche sul vettore dei termini noti) per ripristinare la situazione ideale. Esiste un accorgimento per migliorare il condizionamento del problema: ad ogni passo permutare l'elemento pivot ( $a_{i\:i}^i$) con l'elemento di massimo modulo sulla stessa colonna, da qui viene il nome pivoting parziale. Per permutare una riga con un'altra basta moltiplicare la matrice A a sinistra per una matrice di permutazione che è la matrice identica con due righe scambiate, le quali sono quelle che volevamo scambiare in A. Questo procedimento viene applicato ad ogni passo ed in questo modo si ottiene una fattorizzazione $PA=LU$

[x,A,p]=plu(A,b) - Data la matrice $A$ con determinante non nullo e $\underline{b}$ vettore dei termini noti, la funzione fattorizza $A$ come il prodotto di matrici $LU$ con $U$ triangolare superiore, $L$ triangolare inferiore con elementi diagonali uguali ad 1 e risolve il sistema $LU\underline{x}=P\underline{b}$ con $P$ opportuna matrice di permutazione. La funzione restituisce anche il vettore p contenente i vari indici della permutazione finale delle righe di $A$. La matrice L viene posta nella parte triangolare inferiore di A meno la diagonale principale che deve contenere tutti 1, la matrice U viene posta nela parte triangolare superiore di A.

%PLU
%[x,LU,p]=PLU(A,b)
%Pre: A matrice n x n non singolare, b vettore di lunghezza n
%
% La funzione risolve il sistema lineare di n equazioni in n 
% incognite
% Ax=b con il metodo della fattorizzazione LU con pivoting 
% parziale.
% Restituisce:
% - in x la soluzione del sistema,
% - in LU una matirce che contiene nella parte triangolare 
%   superiore la matrice U 
% e nella parte strettamente 
%   triangolare inferiore gli elementi significativi 
% della
%   matrice L (in pratica L è una matrice che ha la parte 
%   strettamente triangolare
% inferiore uguale ad LU e tutti 1
%   sulla diagonale principale).
% - nel vettore p la permutazione applicata alla matrice A. 
%   La matrice P
% può essere ricostruita permutando le colonne
%   della matrice identità secondo 
%l'ordine descritto dal 
%   vettore p.
% Le matrici L,U,P sono tali che P*A=L*U
%
% See also SIMPLELU, QRHOUSE
function [x,A,p]=PLU(A,b)
n=length(b);
p=1:n; % vettore contenente la permutazione finale
for i=1:n-1
   [m,k]=max(abs(A(i:n,i)));
   % k e' l'indice del massimo elemento in modulo
   % la la sua posizione e' relativa alla dimensione del
   % vettore, quindi si deve aggiustare tale valore
   k=k+i-1;
   % esecuzione di vari swap su
   %%% righe matrice A
   tmp=A(i,:); A(i,:)=A(k,:); A(k,:)=tmp;
   %%% vettore dei termini noti
   tmp=b(i); b(i)=b(k); b(k)=tmp;
   %%% vettore permutazioni
   tmp=p(i); p(i)=p(k); p(k)=tmp;
   for j=i+1:n
      A(j,i)=A(j,i)/A(i,i);
      A(j,i+1:n)=A(j,i+1:n)-A(j,i)*A(i,i+1:n);
   end;
end;
L=tril(A,-1);
for c=1:n
    A(c,c)=1;
end
U=triu(A);
y=solveL(L,b);
x=solveU(U,y);
return

>> A

A =

     1     2     3     4
     5     6     7     8
     9    10    32   354
    65    78    98    54

>> b

b =

     1
     2
    54
     7

>> plu(A,b)

ans =

  -1.20731707317072
   2.38414634146340
  -1.14634146341463
   0.21951219512195

>> inv(A)*b

ans =

  -1.20731707317073
   2.38414634146343
  -1.14634146341463
   0.21951219512195