Soluzione di un sistema di equazioni lineari

9 maggio 2008

1 Soluzione di un sistema di equazioni lineari

Si consideri il sistema lineare di equazioni in incognite

\begin{matrix}_{1 1} _{1} _{1 2} _{2}  _{1} _{}_{1} \\ _{2 1} _{1} _{2 2} _{2}  _{2} _{}_{2} \\ _{3 1} _{1} _{3 2} _{2}  _{3} _{}_{3} \\  \\ _{1} _{1} _{2} _{2}  _{}_{}_{} \end{matrix}(1)

Definita la matrice

; ; ; ; ; \begin{matrix}_{1 1} & _{1 2} & & _{1} \\ _{2 1} & _{2 2} & & _{2} \\ _{3 1} & _{3 2} & & _{3} \\ & & & \\ _{1} & _{2} & & _{} \end{matrix} ; ; ; ; ;

e i vettori colonna delle incognite e dei termini noti

; ; ; ; ; \begin{matrix}_{1} \\ _{2} \\ _{3} \\ \\ _{} \end{matrix} ; ; ; ; ; ; ; ; ; ; \begin{matrix}_{1} \\ _{2} \\ _{3} \\ \\ _{} \end{matrix} ; ; ; ; ;

il sistema di equazioni (1) si può riscrivere come

dove è l’usuale prodotto matrice-vettore per cui ( )_{}_{1}^{}_{}_{}con ( )_{} la -sima componente del vettore .

1.1 Risolubilità di un sistema di equazioni lineari

Soluzione di un sistema di equazioni lineari

Ci si pone il problema se esista e sia unica la soluzione del sistema lineare . A tal proposito distinguiamo alcuni casi 1. : in tal caso la matrice è quadrata di ordine e la soluzione ^{}�

esiste ed è unica se e solo se la matrice è invertibile (che è equivalente a dire che il determinante della matrice A è diverso da 0); per calcolare in Matlab il determinante di una matrice si usa la riga di comando det(A);

2.  oppure  : questa volta la matrice è rettangolare e per la discus-

sione generale si considera il teorema di Rouchè -Capelli, descritto in

     http://it.wikipedia.org/wiki/Teorema_di_Rouch%C3%A9-Capelli

Per il momento, a meno di ipotesi contrarie, considereremo solo matrici quadrate. Inoltre sottintenderemo come d’uso il nel prodotto tra matrici con matrici e matrici con vettori. Così scriveremo invece di e invece di .

1.2 Metodo numerico di risoluzione del sistema lineare, con la fattorizzazione LU

Supponiamo che esistano due matrici quadrate di ordine , diciamo ed , rispettivamente triangolare inferiore e superiore, tali che . Questa scomposizione di A in fattori, è nota come fattorizzazione LU ( sta per lower, sta per upper). Ricordiamo che una matrice quadrata ( _{}) è triangolare inferiore se del tipo

; ; ; ; ; \begin{matrix}_{1 1} & _{1 2} & & _{1} \\ 0 & _{2 2} & & _{2} \\ 0 & 0 & & _{3} \\ & 0 & & \\ 0 & 0 & 0 & _{} \end{matrix} ; ; ; ; ;

mentre è triangolare inferiore se ha la forma

; ; ; ; ; \begin{matrix}_{1 1} & 0 & 0 & 0 \\ _{2 1} & _{2 2} & 0 & 0 \\ _{3 1} & _{3 2} & & 0 \\ & & & \\ _{1} & _{2} & & _{} \end{matrix} ; ; ; ; ;

Ciò significa che tutte le componenti rispettivamente sotto e sopra la diagonale sono uguali a 0.

Essendo  e contemporaneamente  , necessariamente  .
Posto  , si ha  . Se ^{}è la soluzione di questo sistema lineare di tipo
triangolare allora  se e solo se  ^{}. Quindi una volta risolto il sistema

lineare di tipo triangolare superiore ^{}, si è calcolata pure la soluzione del sistema . Visivamente si risolvono di seguito

\begin{matrix} \\  \end{matrix}

In altre parole, se è una matrice invertibile e , si può calcolare l’unica soluzione del sistema lineare risolvendo due sistemi lineari di tipo triangolare. Osserviamo che essendo triangolare inferiore, il vettore ^{}può essere facilmente ricavato con il metodo delle sostituzioni in avanti. In modo simile, essendo triangolare superiore, noto ^{}, il vettore ^{}può essere facilmente ricavato con il metodo delle sostituzioni all’indietro.

Notiamo che non tutte le matrici posseggono una fattorizzazione . Ad esempio per la matrice (di permutazione)

0 1 1 0

non esistono , con le proprietà sopra citate tali che .

Comunque, qualora esista tale fattorizzazione si può provare che senza ulteriori richieste, esistono infinite fattorizzazioni . Sia ad esempio _{1} , _{2} ^{1} , per 0. Se , per _{1} _{1} e _{1} _{2} abbiamo pure _{1} _{1} con _{1}, _{1} che sono rispettivamente triangolare inferiore e superiore.

Si dimostra che se ammette una fattorizzazione , essa è unica se si assume che abbia tutte le componenti diagonali _{}uguali a 1 (fattorizzazione di Doolittle) o, alternativamente, che abbia tutte le componenti diagonali _{}uguali a 1 (fattorizzazione di Crout). Vediamo quale delle due alternative è implementata in Matlab. Eseguiamo il comando help LU dalla shell di Matlab/Octave. Otteniamo:

>> help LU
 LU     LU factorization.
    [L,U] = LU(X) stores an upper triangular matrix in U and a
    "psychologically lower triangular matrix" (i.e. a product
    of lower triangular and permutation matrices) in L, so
    that X = L*U. X can be rectangular.
    [L,U,P] = LU(X) returns lower triangular matrix L, upper
    triangular matrix U, and permutation matrix P so that
    P*X = L*U.
    LU(X), with one output argument, returns the output from
    LAPACK'S DGETRF or ZGETRF routine.
    LU(X,THRESH) controls pivoting in sparse matrices, where
    THRESH is a pivot threshold in [0,1].  Pivoting occurs
    when the diagonal entry in a column has magnitude less than
    THRESH times the magnitude of any sub-diagonal entry in that
    column.  THRESH = 0 forces diagonal pivoting.  THRESH = 1 is
    the default.
    See also LUINC, QR, RREF.
>>

Non dice molto, eccetto che la sintassi è

[L,U,P] = LU(X)

In questo caso con matrice di permutazione (vedremo dopo il significato). Proviamo quale esempio

>> A=[ 7 8 9; 1 2 3; 4 5 6]
A =
     7     8     9
     1     2     3
     4     5     6
>> [L, U, P]=lu(A)
L =
    1.0000         0         0
    0.1429    1.0000         0
    0.5714    0.5000    1.0000

U =
    7.0000    8.0000    9.0000
         0    0.8571    1.7143
         0         0    0.0000
P =
     1     0     0
     0     1     0
     0     0     1
>>

Come previsto la matrice è triangolare inferiore, triangolare superiore e il semplice comando L*U mostra che in effetti

>> L*U
ans =
     7     8     9
     1     2     3
     4     5     6
>>

Visto che è la matrice identica allora che evidentemente coincide con .

Osserviamo che la diagonale di ha componenti uguali a 1. Supporremo quindi sia _{}1 per 1 .

1.2.1 Un implementazione Matlab del metodo LU

Vediamo un implementazione in Matlab di quanto visto, per risolvere il problema , dove al solito supporremo quadrata e non singolare. Come detto precedentemente si tratta di calcolare, se esiste, la fattorizzazione della matrice per poi risolvere i problemi rispettivamente in avanti e all’indietro

\begin{matrix} \\  \end{matrix}

Definiamo le routines

function [x]=sostit_indietro(U,b)
% Risoluzione di un sistema lineare Ux=b con le sostituzioni
% all'indietro: versione per righe.
n = size(U,1);
x = zeros(n,1);
x(n,1) = b(n)/U(n,n);
for k=(n-1):-1:1
  x(k,1)=(b(k) - sum(U(k,k+1:n).*x(k+1:n,1)'))/U(k,k);
end;

function [x]=sostit_avanti(L,b)
% Risoluzione di un sistema lineare Lx=b con le sostituzioni in
% avanti: versione per righe.
n = size(L,1);
x = zeros(n,1);
x(1,1) = b(1)/L(1,1);
for k=2:n
  x(k,1)=(b(k) - sum(L(k,1:k-1).*x(1:k-1,1)'))/L(k,k);
end;

Qualora possegga una fattorizzazione , le matrici ed possono essere ottenute tramite la seguente funzione Matlab/Octave

%  Fattorizzazione LU "sul posto" della matrice A (A=L*U) con la
%  sequenza "kji"
%
%  "sul posto" significa che i fattori L ed U sono ancora contenuti in A .
%
function [A] = fattLUkji(A)
n=size(A,1);
for k=1:n-1
  A(k+1:n,k)=A(k+1:n,k)/A(k,k);
  for j=k+1:n,
    for s=k+1:n

      A(s,j)=A(s,j)-A(s,k)*A(k,j) ;
    end,
  end
end

Come detto nell’help, fattorizzazione sul posto vuol dire che durante la fattorizzazione , l’immagazzinamento parziale delle matrici ed non è eseguito in due ulteriori matrici ed , ma nello spazio precedentemente usato rispettivamente nella parte triangolare inferiore e superiore della matrice . Questo era importante anni fa per poter risolvere sistemi lineari di non piccole dimensioni pur essendo ai limiti della memoria del computer.

Una variante più complessa e stabile è quella vista con il comando lu di Matlab/Octave.

Come commentato nel codice, ed si ottengono dalla matrice in output prendendo la parte triangolare inferiore di e la parte triangolare superiore (imponendo che gli elementi diagonali di siano uguali a 1).

Vediamo un esempio

>> A=hilb(3)
A =
    1.0000    0.5000    0.3333
    0.5000    0.3333    0.2500
    0.3333    0.2500    0.2000
>> b=[1; 0; 0]
b =
     1
     0
     0
>> % RISOLVIAMO IL PROBLEMA Ax=b CON IL METODO LU.
>> [A1] = fattLUkji(A)
A1 =
    1.0000    0.5000    0.3333
    0.5000    0.0833    0.0833
    0.3333    1.0000    0.0056
>>
>> % L E' LA MATRICE CHE SI OTTIENE DALLA PARTE TRIANGOLARE
>> % INFERIORE, TOGLIENDO LA DIAGONALE E AGGIUNGENDO LA
>> % MATRICE IDENTICA.
>> % tril(A) E' LA MATRICE TRIANGOLARE INFERIORE ESTRATTA
>> % ED INCLUDE LA DIAGONALE.
>>
>> L=tril(A1)-diag(diag(A1))+eye(3);
>>
>> % triu(A) E' LA MATRICE TRIANGOLARE SUPERIORE ESTRATTA
>> % ED INCLUDE LA DIAGONALE.
>> U=triu(A1);
>>
>> A
A =

    1.0000    0.5000    0.3333
    0.5000    0.3333    0.2500
    0.3333    0.2500    0.2000
>> L*U
ans =
    1.0000    0.5000    0.3333
    0.5000    0.3333    0.2500
    0.3333    0.2500    0.2000
>>
>> [y]=sostit_avanti(L,b);
>> [x]=sostit_indietro(U,y);
>> A*x-b
ans =
  1.0e-015 *
         0
   -0.8882
    0.8882
>> x
x =
    9.0000
  -36.0000
   30.0000
>>

Ricordiamo che triu e tril estraggono la parte triangolare superiore e inferiore di una matrice. Vediamo un esempio

>> A=hilb(3)
A =
    1.0000    0.5000    0.3333
    0.5000    0.3333    0.2500
    0.3333    0.2500    0.2000
>> triu(A)
ans =
    1.0000    0.5000    0.3333
         0    0.3333    0.2500
         0         0    0.2000
>> tril(A)
ans =
    1.0000         0         0
    0.5000    0.3333         0
    0.3333    0.2500    0.2000
>>

1.3 Tecnica del pivoting

Proviamo il seguente esperimento

>> A=[1 2 3; 4 5 6; 7 8 9]

A =
     1     2     3
     4     5     6
     7     8     9
>> [L,U,P]=lu(A)
L =
    1.0000         0         0
    0.1429    1.0000         0
    0.5714    0.5000    1.0000

U =
    7.0000    8.0000    9.0000
         0    0.8571    1.7143
         0         0    0.0000

P =
     0     0     1
     1     0     0
     0     1     0
>>

notiamo che Matlab/Octave forniscono una fattorizzazione con che non è la matrice identica. Viene da chiedersi quali siano i vantaggi. In molti casi questo è dovuto alla maggiore stabilità dell’algoritmo ma a volte per problemi strettamente matematici: come anticipato, non tutte le matrici posseggono una fattorizzazione . Ad esempio la matrice (detta di permutazione)

0 1 1 0

non possiede la fattorizzazione LU (lo si verifichi manualmente).

Se proviamo la fattorizzazione LU in Matlab (senza pivoting) otteniamo

>> A=[0 1; 1 0]
A =
     0     1
     1     0
>> [L,U]=lu(A)
L =

     0     1
     1     0
U =
     1     0
     0     1
>> L*U
ans =
     0     1
     1     0
>>

ma evidentemente la matrice non è triangolare inferiore.

D’altra parte da un’occhiata veloce al codice di fattLUkji si capisce subito che un primo problema è l’assegnazione

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

in cui il termine ( ) al denominatore può essere nullo. Se la matrice iniziale è invertibile, visto che le trasformazioni apportate non cambiano il valore assoluto del determinante (cosa non ovvia), al passo -simo è sicuramente non nullo uno dei termini al di sotto della diagonale e nella -sima colonna, cioè

( ) ( 1 ) ( )

Se così non fosse si potrebbe provare che la matrice non è invertibile, il che è assurdo in quanto abbiamo supposto det 0. Visivamente saremmo in una condizione del tipo

; ; ; ; ; \begin{matrix}_{1 1} & _{1 2} & _{1 3} & & _{1} \\ 0 & _{2 2} & _{2 3} & & _{2} \\ 0 & 0 & 0 & & _{3} \\ & & & & \\ 0 & 0 & 0 & & _{} \end{matrix} ; ; ; ; ;

ed essendo le prime tre righe linearmente dipendenti il determinante di è necessariamente 0.

Sappiamo quindi che se anche ( ) 0, esiste almeno un indice 1 tale che ( ) 0. Tra tutti gli indici tali che ( ) 0, si scelga per motivi di stabilità ^{}per cui è massimo ( ) . A questo punto si inverta la riga -sima di con la ^{}-sima e si proceda con l’algoritmo della fattorizzazione LU. Così facendo il processo della fattorizzazione LU può essere portato a termine ma si verifica che invece di si ha con matrice di permutazione cioè del tipo ^{1} ove al solito ^{}è la trasposta della matrice .

Nota la fattorizzazione  , se  allora  e quindi  .
Posto ˆ  , abbiamo  ˆ che può essere risolto come abbiamo visto prima

via sostituzioni all’indietro e all’avanti.

Vediamone un esempio che ben illustra il metodo LU con pivoting. Si consideri il sistema lineare (tratto da [5], p.127)

1 2 _{2} _{3} 4 _{4} 13
2 _{1} 0 _{2} 4 _{3} 3 _{4} 28
4 _{1} 2 _{2} 2 _{3} _{4} 20
3 _{1} _{2} 3 _{3} 2 _{4} 6

Riscritto matricialmente per

; ; ; \begin{matrix}1 & 2 & 1 & 4 \\ 2 & 0 & 4 & 3 \\ 4 & 2 & 2 & 4 \\ 3 & 1 & 3 & 2 \end{matrix} ; ; ;

e

; ; ; 13 28 20 6 ; ; ;

si cerca l’unico ^{}tale che ^{}.

Consideriamo un codice che implementa in Matlab tale metodo LU (con pivoting). Una versione più sofisticata è usata da Matlab nelle sue funzioni di algebra lineare e la descrizione che segue ha solo interesse didattico.

function [x,L,U,P]=lupivot(A,b)
%------------------------------------
% INPUT:
% A matrice.
% b termine noto.
%
% OUTPUT:
% x soluzione del sistema Ax=b.
% A matrice elaborata dalla funzione.
% r vettore degli scambi.
%-------------------------------------
[n,n]=size(A);
x=zeros(n,1);
y=zeros(n,1);
c=zeros(1,n);
r=1:n;
for p=1:n-1
    % Calcola il pivot relativo alla "p"-sima colonna.
    [max1,j]=max(abs(A(p:n,p)));
    % Scambia righe "p" e "j".
    c=A(p,:);

    A(p,:)=A(j+p-1,:);
    A(j+p-1,:)=c;
    d=r(p);
    r(p)=r(j+p-1);
    r(j+p-1)=d;
    if A(p,p) == 0
       fprintf('\n \t [WARNING]: A NON INVERTIBILE');
       break
    end
    for k=p+1:n
        mult=A(k,p)/A(p,p);
        A(k,p)=mult;
        A(k,p+1:n)=A(k,p+1:n)-mult*A(p,p+1:n);
    end
end
% Calcola "y".
y(1)=b(r(1));
for k=2:n
    y(k)=b(r(k))-A(k,1:k-1)*y(1:k-1);
end
% Calcola "x".
x(n)=y(n)/A(n,n);
for k=n-1:-1:1
    x(k)=(y(k)-A(k,k+1:n)*(k+1:n))/A(k,k);
end
% Calcoliamo P, L, U.
U=(tril(A'))'; L=A-U+eye(size(A));
P=zeros(size(A));
for k=1:n
    P(k,r(k))=1;
end

Dalla shell di Matlab digitiamo

>> A=[1 2 1 4; 2 0 4 3; 4 2 2 4; -3 1 3 2];
>> b=[13; 28; 20; 6];
>> [x,L,U,P]=lupivot(A,b)

Si ottiene

>> [x,L,U,P]=lupivot(A,b)

x =
    1.6364
   -6.2727
    2.0909
    5.4545

L =
    1.0000         0         0         0
   -0.7500    1.0000         0         0
    0.5000   -0.4000    1.0000         0
    0.2500    0.6000   -0.4583    1.0000

U =
    4.0000    2.0000    2.0000    4.0000
         0    2.5000    4.5000    5.0000
         0         0    4.8000    3.0000
         0         0         0    1.3750

P =
     0     0     1     0
     0     0     0     1
     0     1     0     0
     1     0     0     0

Se digitiamo >>P*A-L*U otteniamo la matrice nulla, e quindi . Certamente è una matrice di permutazione e lo conferma il fatto che se digitiamo >>P*P' otteniamo la matrice identica e quindi l’inversa di è proprio la sua trasposta. Vediamo se in effetti è piccolo. Se fosse ^{}allora sarebbe ^{}0. Scriviamo >>b-A*x e otteniamo

>> b-A*x
ans =
  1.0e-014 *
    0.1776
         0
         0
         0

La soluzione è esatta eccetto per la prima componente in cui viene fatto un errore di circa 10^{15}. Come già detto, potevamo illustrare il problema della risoluzione dei sistemi lineari basandoci su un appropriato comando Matlab, ma non avremo capito la struttura della implementazione (la corrispondente funzione Matlab non è leggibile come file). Comunque,

>> A=[1 2 1 4; 2 0 4 3; 4 2 2 4; -3 1 3 2];
>> b=[13; 28; 20; 6];
>> x1=A\b
x1 =
    1.6364
   -6.2727
    2.0909
    5.4545
>> b-A*x1
ans =
  1.0e-014 *
   -0.3553
         0
   -0.3553
    0.1776
>>

Usando il backslash solo per risolvere sistemi triangolari, il problema può essere risolto come segue

>> A=[1 2 1 4; 2 0 4 3; 4 2 2 4; -3 1 3 2];
>> b=[13; 28; 20; 6];
>> [L,U,P]=lu(A)
L =
    1.0000         0         0         0
   -0.7500    1.0000         0         0
    0.5000   -0.4000    1.0000         0
    0.2500    0.6000   -0.4583    1.0000
U =
    4.0000    2.0000    2.0000    4.0000
         0    2.5000    4.5000    5.0000
         0         0    4.8000    3.0000
         0         0         0    1.3750
P =
     0     0     1     0
     0     0     0     1
     0     1     0     0
     1     0     0     0
>> % PA=LU ALLORA LUx = PAx = Pb.

>> b1=P*b;
>>
>> % PRIMO SISTEMA TRIANGOLARE.
>> y=L\b1;
>>
>> % SECONDO SISTEMA TRIANGOLARE.
>> x=U\y
x =
    1.6364
   -6.2727
    2.0909
    5.4545
>>
>> % CONFRONTO TRA A*x E b.
>> A*x
ans =
   13.0000
   28.0000
   20.0000
    6.0000
>> b
b =
    13
    28
    20
     6
>>

2 Esercizio: fattorizzazione di Cholesky

Data una matrice di ordine simmetrica definita positiva (cioè _{}_{}e ^{} 0 se 0) si dimostra che esiste una matrice triangolare superiore, tale che e _{}1 per 1 . Supposto , : ( ) si dimostra che 1 (cf. [3, p.524], [4, p.278], [7], [8, p.93], [14], [2]).

Si dimostra che è definita positiva se e solo se i suoi autovalori sono positivi. Gli autovalori di si possono calcolare mediante il comando eig. Ad esempio la matrice simmetrica

5 11 11 25

è definita positiva in quanto

>> A=[5 11; 11 25]
A =
     5    11
    11    25
>> eig(A)
ans =

Figure 1: Andrè Louis Cholesky (1875-1918).

0.1339
29.8661
>>

Eseguire una funzione Matlab che verifichi se una matrice generica sia simmetrica e in tal caso esegua la fattorizzazione di Cholesky. Verificare che la quantità

norm(A-R'*R,inf)

sia piccola. Ricordiamo che come si vede dall’help

>> help norm
NORM   Matrix or vector norm.
For matrices...
NORM(X) is the largest singular value of X, max(svd(X)).
NORM(X,2) is the same as NORM(X).
NORM(X,1) is the 1-norm of X, the largest column sum,
                = max(sum(abs(X))).
NORM(X,inf) is the infinity norm of X, the largest row sum,
                = max(sum(abs(X'))).
NORM(X,'fro') is the Frobenius norm, sqrt(sum(diag(X'*X))).
NORM(X,P) is available for matrix X only if P is 1, 2, inf or 'fro'.
For vectors...
NORM(V,P) = sum(abs(V).^P)^(1/P).
NORM(V) = norm(V,2).
NORM(V,inf) = max(abs(V)).
NORM(V,-inf) = min(abs(V)).
See also COND, RCOND, CONDEST, NORMEST.

 Overloaded methods
    help ss/norm.m
    help lti/norm.m
    help frd/norm.m
    help idmodel/norm.m

Se la norma di una matrice di ordine è vicina a 0 allora 0_{}, dove 0_{}indica la matrice nulla di ordine . Si testi la bontà del programma relativamente alla matrice

5 11 11 25

3 Esercizio facoltativo: sulla matrice di Hilbert

Data una matrice il suo numero di condizionamento (in norma 2) è

2 ^{1} _{2}

Sapendo che la norma 2 di una matrice si calcola con il comando norm(A,2), calcolare il numero di condizionamento della matrice di Hilbert ( _{}) \frac{1}{1} con 1 5. Confrontare il risultato con cond(hilb(5)) che vale approssimativamente 4 7661 005. Riprovare lo stesso esperimento per 25.

4 Cenni storici

Il metodo di eliminazione gaussiana probabilmente era noto prima di Gauss. Nel molto interessante [9], si riferisce fosse già noto ai Babilonesi nel 300 a.c. e ai cinesi nel 200 a.c. Per capire come i Babilonesi risolvessero questi problemi si consulti [1, p.16]. In [13], si asserisce sia parte del Chapter Eight, Rectangular Arrays, del testo cinese Jiuzhang suanshu anche noto come The Nine Chapters on the Mathematical Art [12], in cui si studiano 18 problemi con 2 o 5 equazioni. La prima referenza a tale libro è datata 179ac.

Successivamente è presente in modo primitivo in certi studi di Cardano, Leibniz, McLaurin, Bezout e Laplace. In particolare si afferma:

Il termine determinante fu introdotto per primo da Gauss nel Disquisitiones arithmeticae (1801) dedicato allo studio di forme quadratiche. Egli usò tale termine perchè il determinante determina le proprietà della forma quadratica. Comunque il concetto è diverso da quello nostro. Nello stesso lavoro Gauss dispone i coefficienti delle sue forme quadratiche in schieramenti rettangolari. Egli descrive la moltiplicazione di matrice (che pensa come composizione non raggiungendo quindi il concetto di algebra matriciale) e lŠinversa di una matrice nel particolare contesto degli schieramenti di coefficienti delle forme quadratiche. LŠeliminazione Gaussiana, che apparve dapprima nel testo cinese Nove Capitoli dellŠArte Matematica scritto nel 200 a.C., fu usata da Gauss nel suo lavoro dedicato allo studio dellŠorbita

Figure 2: Carl Friedrich Gauss (1777-1855).

dellŠasteroide Pallas. Usando le osservazioni di Pallas fra il 1803 ed il 1809, Gauss ottenne un sistema di sei equazioni lineari in sei incognite. Gauss fornì un metodo sistematico per risolvere tali equazioni che è precisamente lŠeliminazione Gaussiana sulla matrice dei coefficienti.

5 Online

Per ulteriori dettagli sulla fattorizzazione LU si consultino i siti

http://it.wikipedia.org/wiki/Fattorizzazione_LU
http://en.wikipedia.org/wiki/LU_decomposition
http://en.wikipedia.org/wiki/Gaussian_elimination
http://en.wikipedia.org/wiki/The_Nine_Chapters_on_the_Mathematical_Art
http://math.fullerton.edu/mathews/n2003/BackSubstitutionMod.html
http://www.math.tamu.edu/~dallen/masters/egypt_babylon/babylon.pdf
http://it.wikipedia.org/wiki/Gauss
http://it.wikipedia.org/wiki/Decomposizione\_di\_Cholesky
http://it.wikipedia.org/wiki/Andr\%C3\%A9-Louis_Cholesky

References

[1] G.D. Allen, http://www.math.tamu.edu/d̃allen/masters/egypt_babylon/babylon.pdf.

[2] K. Atkinson, An Introduction to Numerical Analysis, Wiley, 1989.

[3] C. Brezinski, http://www.rehseis.cnrs.fr/calculsavant/Textes/biographied’andr.html.

[4] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[5] J.H. Mathews e K.D. Fink, Numerical Methods using Matlab, Prentice Hall, 1999.

[6] The MathWorks Inc., Numerical Computing with Matlab,

    http://www.mathworks.com/moler.

[7] A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.

[8] G. Rodriguez, Algoritmi Numerici, Pitagora Editrice, (2008).

 [9] Wikipedia, http://www.na.iac.cnr.it/even/matrici.htm.
[10] Wikipedia, http://it.wikipedia.org/wiki/Fattorizzazione_LU.
[11] Wikipedia, http://en.wikipedia.org/wiki/LU_decomposition.
[12] Wikipedia,  http://en.wikipedia.org/wiki/The_Nine_Chapters_on
    _the_Mathematical_Art.
[13] Wikipedia, http://en.wikipedia.org/wiki/Gaussian_elimination.
[14] Wikipedia, http://it.wikipedia.org/wiki/Decomposizione_di_Cholesky.
[15] Wikipedia, http://it.wikipedia.org/wiki/Andr%C3%A9-Louis_Cholesky.
[16] Wikipedia, http://it.wikipedia.org/wiki/Gauss.