LU senza pivoting

Guardando attentamentela definizione di $\underline{g}_i$ si nota che viene effettuata una divisione e quindi un'ulteriore precondizione per la fattibilità del problema è che $a_{i\:i}^i\neq 0\; \forall i=1\ldots n$ e questa condizione è equivalente a $\det (A_k) \neq 0 \; \forall k=1 \ldots n$ dove con $A_k$ viene indicato il minore principale di ordine $k$ della matrice A.

[x,A]=simpleLU(A,b) - Data la matrice A tale che $\det A_{k} \neq 0\quad k=1\ldots n$ ed il vettore b, fattorizza la matrice A nel prodotto LU, con L matrice triangolare inferiore con diagonale composta da tutti 1 ed U matrice triangolare superiore. Restituisce in x la soluzione del sistema $LU\underline{x}=\underline{b}$ e memorizza nella parte triangolare superiore di A la parte significativa di U, nella parte inferiore esclusa la diagonale invece inserisce la parte significativa (senza la diagonale) della matrice L.

%SIMPLELU
%[x,LU]=SIMPLELU(A,b)
%Pre: A matrice n x n tale che tutti i suoi minori principali
%     sono non singolari,
%     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.
% 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).
% Le matrici L,U sono tali che A=L*U
%
% See also PLU, QRHOUSE
function [x,A]=simpleLU(A,b)
n=length(b);
for i=1:n-1
   for j=i+1:n
      if A(i,i)==0
         error('Matrice non fattorizzabile LU')
      end
      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)+eye(n);
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

>> inv(A)*b

ans =

  -1.20731707317073
   2.38414634146343
  -1.14634146341463
   0.21951219512195

>> simpleLU(A,b)

ans =

  -1.20731707317073
   2.38414634146341
  -1.14634146341463
   0.21951219512195