Decomposizione LU: differenze tra le versioni
Contenuto cancellato Contenuto aggiunto
Nessun oggetto della modifica |
fix vari |
||
| (79 versioni intermedie di 44 utenti non mostrate) | |||
Riga 1:
In [[algebra lineare]] una '''decomposizione LU''', o '''decomposizione LUP''' o '''decomposizione di Doolittle''' è una [[fattorizzazione]] di una [[
== Idea intuitiva ==
Si immagini di dover risolvere un sistema di equazioni lineari come un puzzle complesso. La decomposizione LU è come scomporre questo puzzle in due puzzle più semplici da risolvere in sequenza.
Una [[matrice]] <math>A</math> può essere "spezzata" in due parti più managevoli:
* una matrice triangolare inferiore <math>L</math> (''Lower'') - ha tutti zeri sopra la diagonale;
* una matrice triangolare superiore <math>U</math> (''Upper'') - ha tutti zeri sotto la diagonale.
È come dire: invece di risolvere direttamente <math>Ax = b</math>, risolviamo prima <math>Ly = b</math> (facile, perché <math>L</math> è triangolare), poi <math>Ux = y</math> (anche questo facile, perché <math>U</math> è triangolare).
=== Perché è utile? ===
La decomposizione LU è vantaggiosa per tre motivi principali:
Riutilizzo computazionale: se si deve risolvere lo stesso sistema con diversi termini noti <math>b</math>, si fa la decomposizione una volta sola e poi si riutilizza. È come assemblare una volta per tutte gli attrezzi giusti, poi usarli per molti lavori simili.
Efficienza numerica: risolvere due sistemi triangolari è molto più veloce che risolvere direttamente il sistema originale. Le matrici triangolari si risolvono "a cascata" dall'alto verso il basso o viceversa.
Stabilità: con il ''pivoting'' (scambi di righe), il metodo diventa numericamente stabile anche per matrici "difficili".
== Come funziona l'algoritmo (versione semplificata) ==
L'algoritmo è sostanzialmente l'[[algoritmo di Gauss|eliminazione gaussiana]], ma invece di dimenticare i passaggi intermedi, vengono inlcusi come parte dell'output:
'''Passo 1''': Trasformare la matrice <math>A</math> in triangolare superiore <math>U</math> usando l'eliminazione gaussiana.
'''Passo 2''': Includere tutti i "moltiplicatori" usati nell'eliminazione in una matrice triangolare inferiore <math>L</math>.
'''Passo 3''': Se servono scambi di righe per la stabilità numerica, tenerne traccia con una matrice di permutazione <math>P</math>.
Risultato: <math>PA = LU</math>, dove <math>P</math> rappresenta gli scambi di righe.
=== Esempio ===
Si consideri il sistema:
:<math>\begin{cases}
2x + 3y = 8 \\
4x + 7y = 18
\end{cases}</math>
In forma matriciale:
:<math>\begin{pmatrix} 2 & 3 \\ 4 & 7 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 8 \\ 18 \end{pmatrix}.</math>
La decomposizione LU produrrebbe:
* <math>L = \begin{pmatrix} 1 & 0 \\ 2 & 1 \end{pmatrix}</math> (il moltiplicatore 2 viene "ricordato");
* <math>U = \begin{pmatrix} 2 & 3 \\ 0 & 1 \end{pmatrix}</math> (matrice triangolare superiore).
== Applicazioni pratiche ==
La decomposizione LU non è solo teoria accademica, ma ha applicazioni concrete:
[[Ingegneria strutturale]]: analisi di telai e travi, dove lo stesso sistema strutturale deve essere analizzato sotto carichi diversi.
[[Circuito elettrico|Circuiti elettrici]]: simulazione di reti elettriche complesse, dove la topologia del circuito resta fissa ma cambiano i valori di corrente e tensione.
[[Computer grafica]]: trasformazioni geometriche e ''rendering'' 3D, dove le stesse matrici di trasformazione vengono applicate a molti punti.
Analisi finanziaria: modelli di portafoglio e ''pricing'' di derivati, dove la struttura di correlazione tra ''asset'' rimane stabile.
Simulazioni fisiche: [[fluidodinamica]], [[meccanica dei solidi]], dove le [[Equazione differenziale|equazioni differenziali]] vengono discretizzate in grandi sistemi lineari che cambiano nel tempo ma mantengono una struttura simile.
== Definizione ==
Sia <math>A</math> una [[matrice invertibile]]. Allora <math>A</math> può essere decomposta come:
:<math>PA = LU</math>
dove
== Idea principale ==
La decomposizione <math>LU</math> è simile all'[[algoritmo di Gauss]]. Nell'eliminazione gaussiana si prova a risolvere l'equazione matriciale:
:<math>Ax = b.</math>
Il processo di eliminazione produce una matrice triangolare superiore <math>U</math> e trasforma il vettore <math>b</math> in <math>b'</math>
:<math>
Poiché
Durante la decomposizione <math>LU</math>,
:<math>Ax = LUx = b,</math>
così si può riutilizzare la decomposizione se si vuole risolvere lo stesso sistema per un differente <math>b</math>.
Nel caso più generale, nel quale la fattorizzazione della matrice comprende anche l'utilizzo di scambi di riga nella matrice, viene introdotta anche una [[matrice di permutazione]] <math>P</math>, ed il sistema diventa:
:<math>\begin{cases}
Ly=Pb\\
Ux=y.
\end{cases}</math>
La risoluzione di questo sistema permette la determinazione del vettore <math>x</math> cercato.
== Algoritmo ==
Applicando delle serie di [[matrici elementari|trasformazioni elementari di matrice]] (cioè moltiplicazioni di <math>A</math> a sinistra) si costruisce una matrice triangolare superiore <math>U</math> che parte da <math>A</math>. Questo metodo è chiamato [[algoritmo di Gauss|metodo di Gauss]]. Queste [[matrici elementari|trasformazioni elementari di matrice]] sono tutte delle [[Trasformazione lineare|trasformazioni lineari]] di tipo combinatorio (il terzo tipo elencato nella voce "[[Matrice elementare#Combinazione lineare|Matrice elementare]]"). Si supponga che <math>T</math> sia il prodotto di <math>N</math> trasformazioni <math>T_N \dots T_2 T_1=T</math>, allora la matrice triangolare superiore è:
:<math>TA=T_N\dots T_2 T_1 A =\colon U.</math>
L'inversa della matrice <math>T</math> è
:<math>T^{-1}=T_1^{-1}T_2^{-1}\dots T_N^{-1}.</math>
Come l'[[algoritmo di Gauss]] usa solo la terza forma dei tre tipi di [[matrici elementari|trasformazioni elementari di matrice]] rendendo <math>A</math> triangolare superiore, si può dedurre che tutte le <math>T_i^{-1}</math> sono triangolari inferiori (vedi [[matrici elementari|trasformazioni elementari di matrice]]). Essendo un prodotto di <math>T_i^{-1}</math> anche:
:<math>T^{-1}=T_1^{-1}T_2^{-1}\dots T_N^{-1}=\colon L</math>
è triangolare inferiore.
:<math>LU=T^{-1}TA=A</math>
== Applicazioni ==
=== Matrici inverse ===
La fattorizzazione <math>PA = LU</math> viene anche usata per calcolare la [[matrice inversa]] di <math>A</math>. Infatti:
:<math>PA = LU</math>
:<math>A = P^{-1}LU,</math>
da cui:
:<math>A^{-1} = (P^{-1}LU)^{-1} = U^{-1}L^{-1}P.</math>
=== Determinante ===
Un altro utilizzo di questa decomposizione è per il calcolo del determinante della matrice <math>A</math>. Infatti:
:<math>A = P^{-1}LU</math>
implica
:<math>\det A = \det(P^{-1}LU),</math>
quindi per il [[teorema di Binet]]:
:<math>\det A = \det(P^{-1})\det L\det U.</math>
Sapendo che il determinante di una [[matrice di permutazione]] vale <math>1</math> se questa corrisponde ad un [[Numeri pari e dispari|numero pari]] di [[Permutazione|permutazioni]], <math>-1</math> altrimenti, e che <math>\det A^{-1} = \frac{1}{\det A}</math>, si ottiene che:
:<math>\det A = (-1)^s\det L\det U,</math>
dove <math>s</math> indica il numero di scambi di riga effettuati nel processo (implicitamente indicati nella matrice <math>P</math>). Sapendo poi che il determinante di una matrice triangolare è dato dal prodotto dei termini sulla sua [[diagonale principale]], si ha che:
:<math>\det A=\left(-1\right)^s\prod_{i=1}^{n}{u_{ii}l_{ii}},</math>
dove i termini <math>u_{ij}</math> e <math>l_{ij}</math> indicano l'elemento nella riga <math>i</math> e colonna <math>j</math> rispettivamente delle matrici <math>U</math> e <math>L.</math> Ma sapendo anche che i termini sulla diagonale principale di <math>L</math> sono tutti <math>1</math>, quindi si può infine concludere:
:<math>\det A=\left(-1\right)^s\prod_{i=1}^{n}u_{ii}.</math>
=== Codice in C ===
<syntaxhighlight lang="c">
/* INPUT: A - vettore di puntatori alle righe della matrice quadrata di dimensione N
* Tol - Callore della tolleranza minima per identificare quando la matrice è prossima a degenerare
* OUTPUT: La matrice A è cambiata, contiene sia le matrici L-E e U tali che A = (L-E)+U e P*A = L*U
* La matrice di permutazione non è salvata in memoria come una matrice, ma in un vettore
di interi di dimensione N+1
* contenente gli indici delle colonne in cui la matrice ha come elementi "1". L'ultimo elemento P[N]=S+N,
* dove S è il numero delle righe scambiate necessarie per il calcolo del determinante, det(P)=(-1)^S
*/
int LUPDecompose(double **A, int N, double Tol, int *P) {
int i, j, k, imax;
double maxA, *ptr, absA;
for (i = 0; i <= N; i++)
P[i] = i; //Matrice di permutazione unaria, P[N] inizializzato con N
for (i = 0; i < N; i++) {
maxA = 0.0;
imax = i;
for (k = i; k < N; k++)
if ((absA = fabs(A[k][i])) > maxA) {
maxA = absA;
imax = k;
}
if (maxA < Tol) return 0; //La matrice è degenerata
if (imax != i) {
//pivoting P
j = P[i];
P[i] = P[imax];
P[imax] = j;
//pivoting delle righe di A
ptr = A[i];
A[i] = A[imax];
A[imax] = ptr;
//Conteggio dei pivot partendo da N per il calcolo del determinante
P[N]++;
}
for (j = i + 1; j < N; j++) {
A[j][i] /= A[i][i];
for (k = i + 1; k < N; k++)
A[j][k] -= A[j][i] * A[i][k];
}
}
return 1; //Decomposizione conclusa
}
/* INPUT: A,P matrici riempite nella funzione LUPDecompose; b - vettore; N - dimensione
* OUTPUT: x - vettore soluzione di A*x=b
*/
void LUPSolve(double **A, int *P, double *b, int N, double *x) {
for (int i = 0; i < N; i++) {
x[i] = b[P[i]];
for (int k = 0; k < i; k++)
x[i] -= A[i][k] * x[k];
}
for (int i = N - 1; i >= 0; i--) {
for (int k = i + 1; k < N; k++)
x[i] -= A[i][k] * x[k];
x[i] = x[i] / A[i][i];
}
}
/* INPUT: A,P matrici riempite nella funzione LUPDecompose; N - dimensione
* OUTPUT: IA è l'inversa della matrice iniziale
*/
void LUPInvert(double **A, int *P, int N, double **IA) {
for (int j = 0; j < N; j++) {
for (int i = 0; i < N; i++) {
if (P[i] == j)
IA[i][j] = 1.0;
else
IA[i][j] = 0.0;
for (int k = 0; k < i; k++)
IA[i][j] -= A[i][k] * IA[k][j];
}
for (int i = N - 1; i >= 0; i--) {
for (int k = i + 1; k < N; k++)
IA[i][j] -= A[i][k] * IA[k][j];
IA[i][j] = IA[i][j] / A[i][i];
}
}
}
/* INPUT: A,P matrici riempite nella funzione LUPDecompose; N - dimensione.
* OUTPUT: La funzione restituisce il valore del determinante della matrice iniziale
*/
double LUPDeterminant(double **A, int *P, int N) {
double det = A[0][0];
for (int i = 1; i < N; i++)
det *= A[i][i];
if ((P[N] - N) % 2 == 0)
return det;
else
return -det;
}
</syntaxhighlight>
=== Codice in C# ===
<syntaxhighlight lang="c#">
public class SystemOfLinearEquations
{
public double[] SolveUsingLU(double[,] matrix, double[] rightPart, int n)
{
// decomposition of matrix
double[,] lu = new double[n, n];
double sum = 0;
for (int i = 0; i < n; i++)
{
for (int j = i; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[i, k] * lu[k, j];
lu[i, j] = matrix[i, j] - sum;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[j, k] * lu[k, i];
lu[j, i] = (1 / lu[i, i]) * (matrix[j, i] - sum);
}
}
// lu = L+U-I
// Calcolo delle soluzioni di Ly = b
double[] y = new double[n];
for (int i = 0; i < n; i++)
{
sum = 0;
for (int k = 0; k < i; k++)
sum += lu[i, k] * y[k];
y[i] = rightPart[i] - sum;
}
// Calcolo delle soluzioni di Ux = y
double[] x = new double[n];
for (int i = n - 1; i >= 0; i--)
{
sum = 0;
for (int k = i + 1; k < n; k++)
sum += lu[i, k] * x[k];
x[i] = (1 / lu[i, i]) * (y[i] - sum);
}
return x;
}
}
</syntaxhighlight>
=== Codice in Matlab ===
<syntaxhighlight lang="matlab">
function x = SolveLinearSystem(A, b)
n = length(b);
x = zeros(n, 1);
y = zeros(n, 1);
% Decomposizione della matrice per mezzo del metodo Doolittle
for i = 1:1:n
for j = 1:1:(i - 1)
alpha = A(i,j);
for k = 1:1:(j - 1)
alpha = alpha - A(i,k)*A(k,j);
end
A(i,j) = alpha/A(j,j);
end
for j = i:1:n
alpha = A(i,j);
for k = 1:1:(i - 1)
alpha = alpha - A(i,k)*A(k,j);
end
A(i,j) = alpha;
end
end
% A = L+U-I
% calcolo delle soluzioni di Ly = b
for i = 1:1:n
alpha = 0;
for k = 1:1:i
alpha = alpha + A(i,k)*y(k);
end
y(i) = b(i) - alpha;
end
% calcolo delle soluzioni di Ux = y
for i = n:(-1):1
alpha = 0;
for k = (i + 1):1:n
alpha = alpha + A(i,k)*x(k);
end
x(i) = (y(i) - alpha)/A(i, i);
end
end
</syntaxhighlight>
== Bibliografia ==
* {{en}} Bunch, James R.; Hopcroft, John (1974), "Triangular factorization and inversion by fast matrix multiplication", ''Mathematics of Computation'' 28: 231–236, ISSN 0025-5718.
* {{en}} Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford (2001), ''Introduction to Algorithms'', MIT Press and McGraw-Hill, ISBN 978-0-262-03293-3.
* {{en}} Golub, Gene H.; Van Loan, Charles F. (1996), ''Matrix Computations (3rd ed.)'', Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
* {{en}} Horn, Roger A.; Johnson, Charles R. (1985), ''Matrix Analysis'', Cambridge University Press, ISBN 0-521-38632-2. See Section 3.5.
* {{en}} Householder, Alston S. (1975), ''The Theory of Matrices in Numerical Analysis'', New York: Dover Publications, MR 0378371.
* {{en}} Okunev, Pavel; Johnson, Charles R. (1997), ''Necessary And Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix'', arXiv:math.NA/0506382.
* {{en}} Poole, David (2006), ''Linear Algebra: A Modern Introduction (2nd ed.)'', Canada: Thomson Brooks/Cole, ISBN 0-534-99845-3.
* {{en}} Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Section 2.3", ''Numerical Recipes: The Art of Scientific Computing (3rd ed.)'', New York: Cambridge University Press, ISBN 978-0-521-88068-8.
* {{en}} Trefethen, Lloyd N.; Bau, David (1997), ''Numerical linear algebra'', Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.
== Voci correlate ==
* [[Analisi numerica]]
* [[Decomposizione di Cholesky]]
* [[Decomposizione di Crout]]
* [[Decomposizione di una matrice]]
* [[Decomposizione QR]]
* [[Matrice triangolare]]
== Altri progetti ==
{{interprogetto}}
== Collegamenti esterni ==
* {{Collegamenti esterni}}
* {{en}}[http://www.math-linux.com/spip.php?article51 LU decomposition] on ''Math-Linux''.
* {{en}}[https://web.archive.org/web/20070609042305/http://math.fullerton.edu/mathews/n2003/LUFactorMod.html Module for LU Factorization with Pivoting], Prof. J. H. Mathews, California State University, Fullerton
* {{en}}[http://numericalmethods.eng.usf.edu/topics/lu_decomposition.html LU decomposition] at ''Holistic Numerical Methods Institute''
* {{en}}[https://www.mathworks.com/help/techdoc/ref/lu.html LU matrix factorization]. MATLAB reference.
* {{cita web|1=http://sole.ooz.ie/|2=WebApp descriptively solving systems of linear equations with LU Decomposition|lingua=en|urlmorto=sì|accesso=10 aprile 2014|urlarchivio=https://archive.is/20110425223706/http://sole.ooz.ie/|dataarchivio=25 aprile 2011}}
* {{en}}[https://web.archive.org/web/20081212221215/http://www.bluebit.gr/matrix-calculator/ Matrix Calculator], bluebit.gr
* {{en}}[https://web.archive.org/web/20081020061504/http://www.uni-bonn.de/~manfear/matrix_lu.php LU Decomposition Tool], uni-bonn.de
* {{en}}[http://demonstrations.wolfram.com/LUDecomposition/ LU Decomposition] by Ed Pegg, Jr., The Wolfram Demonstrations Project, 2007.
== Voci correlate ==
* [[Analisi numerica]]
* [[Decomposizione di Cholesky]]
* [[Decomposizione di Crout]]
* [[Decomposizione di una matrice]]
* [[Decomposizione QR]]
* [[Matrice triangolare]]
== Altri progetti ==
{{interprogetto}}
== Collegamenti esterni ==
* {{Collegamenti esterni}}
* {{en}}[http://www.math-linux.com/spip.php?article51 LU decomposition] on ''Math-Linux''.
* {{en}}[https://web.archive.org/web/20070609042305/http://math.fullerton.edu/mathews/n2003/LUFactorMod.html Module for LU Factorization with Pivoting], Prof. J. H. Mathews, California State University, Fullerton
* {{en}}[http://numericalmethods.eng.usf.edu/topics/lu_decomposition.html LU decomposition] at ''Holistic Numerical Methods Institute''
* {{en}}[https://www.mathworks.com/help/techdoc/ref/lu.html LU matrix factorization]. MATLAB reference.
* {{cita web|1=http://sole.ooz.ie/|2=WebApp descriptively solving systems of linear equations with LU Decomposition|lingua=en|urlmorto=sì|accesso=10 aprile 2014|urlarchivio=https://archive.is/20110425223706/http://sole.ooz.ie/|dataarchivio=25 aprile 2011}}
* {{en}}[https://web.archive.org/web/20081212221215/http://www.bluebit.gr/matrix-calculator/ Matrix Calculator], bluebit.gr
* {{en}}[https://web.archive.org/web/20081020061504/http://www.uni-bonn.de/~manfear/matrix_lu.php LU Decomposition Tool], uni-bonn.de
* {{en}}[http://demonstrations.wolfram.com/LUDecomposition/ LU Decomposition] by Ed Pegg, Jr., The Wolfram Demonstrations Project, 2007.
{{Portale|matematica}}
[[Categoria:Decomposizione matriciale]]
| |||