Laboratorio 4 : L’Eliminazione di Gauss il Metodo LU

L”eliminazione gaussiana (Carl Friedrich Gauss, 1777-1855), noto anche come algoritmo di riduzione per le righe, è un algoritmo per risolvere sistemi quadrati di equazioni lineari. Come avete visto a lezione si riduce ad una sequenza di operazioni eseguite sulla corrispondente matrice di coefficienti.

Alcune estensioni di questo metodo possono essere utilizzate anche per calcolare il rango di una matrice, il determinante di una matrice quadrata e l”inversa di una matrice invertibile.

Come abbiamo visto nel secondo Laboratorio sul Condizionamento di un sistema lineare la precisione ottenuta in uscita da questa procedura. Ricordiamo che per valutare il numero di condizionamento di una matrice \(A\), MATLAB mette a disposizione due comandi:

• cond: Numero di condizione rispetto all’inversione. cond(A) restituisce il numero di condizionamento rispetto alla norma 2 (che ricordiamo essere il rapporto tra valore singolare più grande di \(A\) e di quello più piccolo);

• rcond: stimatore del reciproco del numero di condizionamento in norma 1. Utilizza lo stimatore implementato in LAPACK;

• condest: calcola un limite inferiore \(C\) per il numero di condizionamento in norma 1 di una matrice quadrata \(A\). Accetta come seconda opzione un parametro intero positivo \(T\) che rapprenta il numero di colonne utilizzate per l’iterazione dell’algoritmo di stima sottostante. L’aumentare del numero di colonne fornisce di solito una stima migliore del condizionamento, al prezzo di aumentarne il costo computazionale. Il valore predefinito è \(T = 2\), che fornisce quasi sempre una stima corretta entro un fattore 2.

Si provi ad eseguire il seguente comando:

A = gallery('grcar',8);
[cond(A,1) 1/rcond(A) condest(A)]
[cond(A,1) 1/rcond(A) condest(A)]

Come illustra questo esempio, condest non restituisce necessariamente lo stesso risultato ad ogni chiamata, poiché fa uso della funzione rand al suo interno (si tratta di uno stimatore stocastico).

I comandi per le fattorizzazioni in Matlab

Se \(A\) è una matrice \(n\) per \(n\) su un campo \(\mathbb{K}\) (\(A \in M_n(\mathbb{K})\)), allora \(A\) si dice «avere una fattorizzazione \(LU\)» se esiste una matrice triangolare inferiore \(L \in M_n(\mathbb{K})\) e una matrice triangolare superiore \(U \in M_n(\mathbb{K})\) tale che

\[A = LU.\]

Questo è particolarmente utile, poiché, se \(A = LU\), allora ogni soluzione \(\mathbf{x}\) che sia soluzione di:

\[U\mathbf{x} = \mathbf{y},\]

per una soluzione \(\mathbf{y}\) di

\[L\mathbf{y} = \mathbf{b},\]

e se \(L\) è invertibile, allora \(\mathbf{y}\) è unicamente determinata; e se \(U\) è invertibile, allora ogni soluzione \(\mathbf{y}\) di \(L\mathbf{y} = \mathbf{b}\) produrrà una soluzione \(\mathbf{x}\) di \(A\mathbf{x} = \mathbf{b}\) via \(U\mathbf{x} = \mathbf{y}\).

Pericolo

In generale non tutte le matrici \(A\) ammettono una fattorizzazione \(LU\) di questa forma, si consideri ad esempio la matrice:

\[\begin{split}A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} l_{11} & 0 \\ l_{12} & l_{22} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} \\ 0 & u_{22} \end{bmatrix} = LU,\end{split}\]

questa richiederebbe di avere gli elementi diagonali \(u_{ii}\), \(l_{ii}\), \(i=1,2\), diversi da \(0\), ma questo contraddirebbe \(0 = a_{11} = l_{11}u_{11}\).

Una condizione necessaria e sufficiente affinché una fattorizzazione di questa forma esista per una matrice \(A \in M_n(\mathbb{K})\), \(\mathbb{K} = \mathbb{R},\mathbb{C}\), è che

\[\begin{split}\operatorname{rank}(A_{11}) + k \geq \operatorname{rank}([A_{11}\,A_{12}]) + \operatorname{rank}\left(\begin{bmatrix}A_{11}\\A_{21}\end{bmatrix}\right),\end{split}\]

per \(k = 1,\ldots,n-1\) e

\[\begin{split}A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \qquad A_{11} \in M_k(\mathbb{K}).\end{split}\]

Quando esiste la decomposizione \(LU\) non è unica (le combinazioni di \(L\) e \(U\) per una data \(A\) sono infinite). Per determinare un algoritmo univoco è necessario porre determinati vincoli su \(L\) o \(U\) (Tabella 9).

Varianti dell’algoritmo \(LU\).

Tabella 9 Varianti dell’algoritmo \(LU\).

Nome	Vincoli
Doolittle	\(L_{ii} = 1\), \(i=1,\ldots,n\)
Crout	\(U_{ii} = 1\), \(i=1,\ldots,n\)
Choleski	\(U = L^T\)

Uso del comando lu in MATLAB

MATLAB fornisce il comando lu per calcolare la fattorizzazione \(LU\) di una matrice. Questa funzione è ottimizzata, stabile numericamente e offre diverse modalità d’uso a seconda delle esigenze.

La forma più semplice del comando è:

[L,U] = lu(A)

che restituisce una matrice triangolare superiore U e una matrice triangolare inferiore permutata L tali che A = L*U. Questa forma è utile per computazioni rapide ma non fornisce esplicitamente la matrice di permutazione.

Per ottenere la classica fattorizzazione \(PA = LU\), si utilizza:

[L,U,P] = lu(A)

Questa forma restituisce:

• L: matrice triangolare inferiore con diagonale unitaria

• U: matrice triangolare superiore

• P: matrice di permutazione tale che P*A = L*U

La matrice P può essere restituita in due formati:

[L,U,P] = lu(A,'matrix')  % P è una matrice (default)
[L,U,P] = lu(A,'vector')  % P è un vettore di indici

Nel formato vettoriale, P contiene gli indici di permutazione tali che A(P,:) = L*U.

A = [3 -1  4;
  -2  0  5;
   7  2 -2];
   
% Fattorizzazione base
[L,U] = lu(A);
norm(A - L*U)  % Verifica: dovrebbe essere piccolo

% Con matrice di permutazione
[L,U,P] = lu(A);
norm(P*A - L*U)  % Verifica: dovrebbe essere ~0

% Formato vettoriale
[L,U,p] = lu(A,'vector');
norm(A(p,:) - L*U)  % Equivalente alla forma matriciale

Il comando lu di MATLAB implementa automaticamente:

• Pivoting parziale per garantire stabilità numerica

• Ottimizzazioni per matrici dense e sparse

• Scaling automatico quando necessario

• Gestione efficiente della memoria

Per questa ragione, è generalmente preferibile usare lu piuttosto che implementazioni personalizzate, a meno che non si abbiano esigenze didattiche o di ricerca particolari.

Fattorizzazione per matrici sparse

Definizione 2

Una matrice sparsa o array sparso è una matrice in cui la maggior parte degli elementi è zero. Per dare una definizione rigorosa riguardo alla proporzione di elementi di valore zero affinché una matrice sia qualificata come sparsa è tipicamente necessario considerare la sorgente del problema che la ha generata, e.g., la discretizzazione di una PDE. Un criterio comune è che il numero di elementi diversi da zero sia approssimativamente uguale al numero di righe o colonne.

Al contrario, se la maggior parte degli elementi è diversa da zero, la matrice è considerata densa.

Per matrici sparse, MATLAB utilizza una struttura dati ottimizzata che memorizza solo gli elementi non nulli e le loro posizioni. Questo approccio permette di:

• Ridurre significativamente l’uso di memoria

• Accelerare le operazioni che coinvolgono molti zeri

• Gestire matrici di dimensioni molto grandi che sarebbero impraticabili in formato denso

% Esempio: matrice sparsa tridiagonale
n = 1000;
A_sparse = spdiags(ones(n,3)*[1 -2 1], -1:1, n, n);

% Informazioni sulla sparsità
fprintf('Elementi totali: %d\n', n^2);
fprintf('Elementi non nulli: %d\n', nnz(A_sparse));
fprintf('Percentuale di sparsità: %.2f%%\n', ...
  100*(1 - nnz(A_sparse)/n^2));

% Confronto di memoria
whos A_sparse  % Mostra l'uso di memoria

% Conversione da densa a sparsa e viceversa
A_full = full(A_sparse);  % Converte in formato denso
A_sp = sparse(A_full);    % Converte in formato sparso

Visualizzare una matrice sparsa.

Il comando spy in MATLAB è utilizzato per visualizzare la struttura di una matrice sparsa. Mostra un grafico in cui ogni punto rappresenta un elemento non nullo della matrice, permettendo di avere una rapida comprensione della distribuzione degli elementi non nulli. La sintassi di base è:

spy(A_sparse)

Dove A_sparse è la matrice sparsa che si desidera visualizzare. Questo comando è particolarmente utile per analizzare la sparsità e la struttura delle matrici in applicazioni ingegneristiche e scientifiche.

Esempi comuni di matrici sparse di interesse ingegneristico includono:

• Matrici tridiagonali: derivanti dalla discretizzazione di equazioni differenziali con il metodo delle differenze finite (es. problemi di conduzione termica, diffusione)

• Matrici a banda: tipiche dell’analisi strutturale con elementi finiti per travi, telai e strutture reticolari

• Matrici di connettività: rappresentanti reti elettriche, idrauliche o di trasporto mediante grafi

• Laplaciani discreti: dalla discretizzazione di operatori differenziali in problemi di campo (elettromagnetismo, fluidodinamica)

Per matrici sparse, MATLAB offre una versione più efficiente con doppio pivoting:

[L,U,P,Q] = lu(A)

che calcola la fattorizzazione \(PAQ = LU\), dove sia P che Q sono matrici di permutazione. Questa variante minimizza il fill-in (creazione di nuovi elementi non nulli) durante la fattorizzazione, rendendo il calcolo significativamente più efficiente per matrici sparse di grandi dimensioni.

Soluzione diretta del sistema senza passare per la fattorizzazione

La funzione per risolvere i sistemi lineari è la funzione mldivide \. Si tratta di una delle funzioni più elaborate di MATLAB, il comando A\B è la divisione matriciale di A in B, moralmente è la stessa operazione di INV(A)*B, ma è calcolata in modo diverso, stabile ed efficiente.

Se \(A\) è una matrice \(N\) per \(N\) e \(\mathbf{b}\) è un vettore colonna con \(N\) componenti, o una matrice con più di queste colonne, allora X = A\b è la soluzione dell’equazione A*X = B.

Avvertimento

Viene stampato un messaggio di avviso se A è scalata male o malcondizionata, ad esempio

x = hilb(15)\ones(15,1);

restituisce

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND =  5.460912e-19.

Implementare la fattorizzazione LU (facoltativo)

Come primo esercizio implementeremo la fattorizzazione \(LU\) nella forma di Doolittle, che è quella più strettamente legata all’algoritmo di eleminazione di Gauss che conoscete sin dal corso di algebra lineare.

Se scriviamo l’uguaglianza \(A = LU\) con il vincolo aggiuntivo per cui \(L_{ii} = 1\), \(i=1,\ldots,n\), possiamo ricavare che i termini della matrice \(U\) sono dati da

(6)\[\begin{split}\forall\,j\quad & \begin{array}{ll} i = 1, & u_{i,j} = a_{i,j},\\ i > 1, & u_{i,j} = \displaystyle a_{i,j} - \sum_{k=1}^{i-1} l_{ik} u_{k,j}, \end{array}\end{split}\]

mentre quelli della matrice \(L\) sono dati da:

(7)\[\begin{split}\forall\,i\quad & \begin{array}{ll} j = 1, & l_{ij} = \frac{a_{ij}}{u_{jj}}, \\ j > 1, & l_{ij} = \displaystyle \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} }{u_{jj}}. \end{array}\end{split}\]

Esercizio

Si usino le relazioni (6) e (7) per implementare la versione di Doolittle dell’algoritmo \(LU\) per una matrice quadrata. Un prototipo della funzione è dato da:

function [L,U] = doolittlelu(A)
%%DOOLITTLELU Calcola la fattorizzazione LU della matrice A con i vincolo di
%Doolittle.

% Controllo degli input:

% Pre-allocazione memoria:
L = eye(n,n);
U = zeros(n,n);
%!! Questa potrebbe essere anche:
% L = zeros(n,n);
% U = zeros(n,n);
% Dipende da come decidete di utilizzare le formule...

% Utilizzo delle formule per l_{ij} e u_{ij}.

end

• Si faccia un controllo dell’input: la matrice \(A\) è quadrata?

• Si utilizzino operazioni vettorizzate per le somme che compaiono in (6) e (7)

Per verificare la buona riuscita dell’implementazione si può usare il seguente problema di test.

%% Test della versione di Doolittle dell'algoritmo LU  

clear; clc;

A = [ 3 -1  4;
     -2  0  5;
      7  2 -2 ];
[L,U] = doolittlelu(A);
fprintf("|| A - LU || = %1.2e\n",norm(A - L*U,2));

% Vediamone gli ingressi:
format rat
disp(A)
disp(L)
disp(U)
disp(L*U)
format short

Suggerimento

Guardando alla formulazione delle equazioni (6) e (7) potremmo osservare che in realtà possiamo sovrascrivere le entrate della matrice \(A\) con la fattorizzazione \(LU\). Se avete concluso la prima parte pensate a come è possibile realizzare questa versione.

Un suggerimento più estensivo

Possiamo sviluppare il ciclo di istruzioni sostanzialmente in due ordine, in un primo procediamo alla fattorizzazioni colonna per colonna, ovvero:

for k=1:n
    L(k:n,k) = A(k:n,k) / A(k,k); % Dividiamo per il Pivot
    U(k,1:n) = A(k,1:n);          % Aggiorniamo le entrate della U
    A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);
end
U(:,end) = A(:,end);

Nel caso dell”esempio vediamo:

\[\begin{split}k = 0, \qquad L = \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right) \qquad U = \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right), \qquad A = \left(\begin{array}{ccc} 3 & -1 & 4\\ -2 & 0 & 5\\ 7 & 2 & -2 \end{array}\right)\end{split}\]

\[\begin{split}k = 1, \qquad L^{(1)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ -\frac{2}{3} & 0 & 0\\ \frac{7}{3} & 0 & 0 \end{array}\right), \qquad \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right) \qquad A^{(1)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & \frac{13}{3} & -\frac{34}{3} \end{array}\right), \\\end{split}\]

\[\begin{split}k = 2, \qquad L^{(2)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ \frac{7}{3} & -\frac{13}{2} & 0 \end{array}\right), \qquad U^{(2)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & 0 \end{array}\right), \qquad A^{(2)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & \frac{77}{2} \end{array}\right), \\\end{split}\]

\[\begin{split}k = 3, \qquad L^{(3)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ \frac{7}{3} & -\frac{13}{2} & 1 \end{array}\right), \qquad U^{(3)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & \frac{77}{2} \end{array}\right), \qquad A^{(3)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & \frac{77}{2} \end{array}\right)\end{split}\]

La seconda alternativa è quella di eseguire questo ciclo procedendo lungo le righe della matrice

for i = 1:1:n
    for j = 1:(i - 1)
        L(i,j) = (A(i,j) - L(i,1:(j - 1))*U(1:(j - 1),j)) / U(j,j);
    end
    j = i:n;
    U(i,j) = A(i,j) - L(i,1:(i - 1))*U(1:(i - 1),j);
end

Nel caso dell”esempio vediamo:

\[\begin{split}i = 0, \qquad L = \left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right), \qquad U = \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right), \qquad A = \left(\begin{array}{ccc} 3 & -1 & 4\\ -2 & 0 & 5\\ 7 & 2 & -2 \end{array}\right),\end{split}\]

\[\begin{split}i = 1, \qquad L^{(1)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right), \qquad U^{(1)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right),\end{split}\]

\[\begin{split}i = 2, \qquad L^{(2)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ 0 & 0 & 1 \end{array}\right), \qquad U^{(2)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & 0 \end{array}\right)',\end{split}\]

\[\begin{split}i = 3, \qquad L^{(3)} = \left(\begin{array}{ccc} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ \frac{7}{3} & -\frac{13}{2} & 1 \end{array}\right), \qquad U^{(3)} = \left(\begin{array}{ccc} 3 & -1 & 4\\ 0 & -\frac{2}{3} & \frac{23}{3}\\ 0 & 0 & \frac{77}{2} \end{array}\right)',\end{split}\]

L”unica differenza tra le due implementazioni è l’ordine, per colonna e per riga rispettivamente, in cui facciamo le operazioni.

La fattorizzazione con pivoting

Consideriamo il seguente esempio numerico in cui andiamo a confrontare gli errori assoluti e relativi,

\[E_{\text{abs}} = \| A - LU\|_2, \qquad E_{\text{rel}} = \frac{\| A - LU\|_2}{\| A\|_2},\]

tra la matrice \(A\) ed il prodotto dei termini \(LU\) ottenuto con l’algoritmo implementato nell’esercizio precedente.

clear; clc; close all;

sizes = 100:100:1500;
N = length(sizes);
maxerr = zeros(N,1);
relerr = zeros(N,1);
for i = 1:N
   A = rand(sizes(i),sizes(i));
   [L,U] = doolittlelu(A);
   [Lm,Um] = lu(A);
   normA = norm(A,2);
   maxerr(i) = norm(A - L*U,2);
   relerr(i) = maxerr(i)/normA;
   fprintf("Dimensione %d\n\t || A - LU || = %1.2e\n",sizes(i),maxerr(i));
   fprintf("\t || A - LU ||/||A|| = %1.2e\n",relerr(i));
end

figure(1)
loglog(sizes,maxerr,'o-',sizes,relerr,'x-','LineWidth',2);
xlabel('Dimensione');
legend({'Errore Assoluto','Errore Relativo'},...
    'Location','best','FontSize',14);

Quello che osserviamo nella figura è una crescita di più di un’ordine di grandezza nella costruzione fatta dalla fattorizzazione \(LU\) accompagnato da un certo numero di oscillazioni Fig. 5.

Fig. 5 Errori assoluti e relativi tra una matrice \(A\) ad entrate casuali uniformemente distribuite in \([0,1]\) e la relativa fattorizzazione \(LU\) calcolata senza pivoting.

Ci aspetteremo che almeno l”errore relativo si comportasse meglio. Quello che sta accadendo è che la scelta dei pivot affidata alla forma originale della matrice sta introducendo dell’instabilità numerica che ci fa perdere in precisione relativa più di quanto sarebbe lecito aspettarsi guardando ai numeri di condizionamento delle operazioni in questione.

Per migliorare la stabilità dell’algoritmo e fare sì che la fattorizzazione \(LU\) si possa calcolare per una generica matrice (ovvero anche quando si incontrano pivot nulli) possiamo introdurre una strategia di scambio delle righe della matrice detta, per l’appunto, di pivoting parziale.

Formalmente, questo vuol dire che otterremo una fattorizzazione della forma

\[P A = L U,\]

in cui \(P\) è una matrice di permutazione. Strutturalmente una matrice di permutazione non è nient’altro che una matrice identità di cui abbiamo permutato le righe. Vediamo un veloce esempio:

format rat
H = hilb(5)
P = eye(5)
P = P([5:-1:1],:)
P*H

che costruisce le matrici

\[\begin{split}H = \left(\begin{array}{ccccc} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5}\\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8}\\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{array}\right), \quad P = \left(\begin{array}{ccccc} 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 \end{array}\right), \quad PH = \left(\begin{array}{ccccc} \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9}\\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8}\\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\ 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{array}\right),\end{split}\]

e ci mostra anche come possiamo permutare due (o più) righe di una matrice MATLAB. Abbiamo ora tutti gli ingredienti necessari a costruire la nostra versione fatta in casa della fattorizzazione LU con pivoting in MATLAB.

Esercizio

Si costruisca una funzione ludecomp per costruire la fattorizzazione \(LU\) con pivoting parziale per una matrice \(A\). Un prototipo della funzione è quindi:

function [L,U,P] = ludecomp(A)
%%LUDECOMPP Fattorizzazione LU con pivoting parziale.

% Contrllo dell'input:

% Allocazione della memoria:
L = zeros(n);
U = zeros(n);
P = eye(n);

% Fattorizzazione con Pivoting
end

• Si faccia un controllo dell’input: la matrice \(A\) è quadrata?

• Si utilizzino operazioni vettorizzate per le somme che compaiono in (6) e (7).

Suggerimento

Se il ciclo esterno è per \(k=1,\ldots,n\) possiamo individuare il pivot facendo:

[~,r] = max(abs(A(k:end,k)));
r = n-(n-k+1)+r; % Trasliamo l'indice in questione

Qual è l’effetto della traslazione di indice?

Possiamo verificare l’implementazione della versione con pivoting parziale facendone un paragone con la versione dell’esercizio precedente con il seguente codice:

sizes = 100:100:1500;
N = length(sizes);
maxerr = zeros(N,1);
relerr = zeros(N,1);
maxerrm = zeros(N,1);
relerrm = zeros(N,1);
for i = 1:N
   A = rand(sizes(i),sizes(i));
   [L,U] = doolittlelu(A);
   [Lm,Um,P] = ludecomp(A);
   normA = norm(A,2);
   maxerr(i) = norm(A - L*U,2);
   relerr(i) = maxerr(i)/normA;
   maxerrm(i) = norm(P*A - Lm*Um,2);
   relerrm(i) = maxerrm(i)/normA;
   fprintf("Dimensione %d\n \t || A - LU || = %1.2e\n",sizes(i),maxerr(i));
   fprintf("\t || A - LU ||/||A|| = %1.2e\n",relerr(i));
   fprintf("\t || P*A - LU ||/||A|| = %1.2e\n",maxerrm(i));
   fprintf("\t || P*A - LU ||/||A|| = %1.2e\n",relerrm(i));
end

figure(2)
loglog(sizes,maxerr,'o-',sizes,relerr,'x-',...
    sizes,maxerrm,'o-',sizes,relerrm,'x-','LineWidth',2);
xlabel('Dimensione');
legend({'Errore Assoluto','Errore Relativo',...
    'Errore Assoluto (Pivoting)','Errore Relativo (Pivoting)'},...
    'Location','best','FontSize',14);

Con questa esecuzione possiamo osservare che abbiamo mantenuto un errore relativo che è all’incirca della precisione di macchina, migliorando sensibilmente le performance rispetto alla versione che utilizzava l’ordinamento di partenza.

Fig. 6 Errori assoluti e relativi tra una matrice \(A\) ad entrate casuali uniformemente distribuite in \([0,1]\) e la relativa fattorizzazione \(LU\) calcolata senza e con pivoting.

Calcolo del determinante

Dopo aver calcolato la fattorizzazione \(LU\) di una matrice \(A\) è possibile utilizzarla anche per calcolare il determinante della matrice \(A\). Si può fare sfruttando il teorema di Binet che ricordiamo collega il prodotto fra matrici quadrate con il determinante:

\[\det(AB) = \det(A)\det(B), \qquad \forall\,A,B \in M_n(\mathbb{K}),\]

ed, in particolare,

\[\det(A) = \det(LU) = \det(L)\det(U),\]

che, se usiamo la formulazione di Doolittle si riduce quindi ad essere:

\[\det(A) = \det(LU) = \det(L)\det(U) = \prod_{i=1}^{n} l_{ii} \prod_{i=1}^n u_{ii} = 1 \prod_{i=1}^n u_{ii} = \prod_{i=1}^n u_{ii}.\]

Che possiamo riportare facilmente su MATLAB come:

function res = ludet(A)
%%LUDET Calcola il determinante della matrice A sfruttando la sua fattorizzazione
%LU in forma di Doolittle e senza pivoting.
[~,U] = doolittlelu(A);
res = prod(diag(U));
end

e testare rapidamente con:

A = pascal(5,2);
fprintf('Errore relativo: %1.2e\n',abs(ludet(A) - det(A))/abs(det(A)));

che ci dovrebbe restituire un errore relativo, rispetto alla funzione det nativa di MATLAB dell’ordine di 3.33e-16.

Avvertimento

Per costruire il calcolo del determinante anche nel caso in cui si utilizzi la fattorizzazione \(LU\) con pivoting è necessario tenere traccia del numero di scambi di righe effettuati.

Esercizi

Consideriamo alcune applicazioni ingegneristiche del problema calcolare le soluzione di un sistema [Kiu15].

Esercizio

La formulazione per gli spostamenti di una trave reticolare piana è simile a quella di un sistema di masse e molle. Le differenze sono che

1. i termini di rigidità sono dati da \(k_{i} = (E A/L)_ {i}\), dove \(E\) è il modulo di elasticità, \(A\) rappresenta l’area della sezione di taglio e \(L\) è la lunghezza dell’elemento;

2. ci sono due componenti dello scostamento in ogni punto di giuntura.

Per la trave reticolare in figura:

Fig. 7 Profilo di un elemento della trave reticolare piana.

si ottiene quindi il seguente sistema di equazioni lineari

\[K \mathbf{u} = \mathbf{p}\]

dove

K = [27.58  7.004 -7.004 0 0
     7.004 29.57  -5.253 0 -24.32
    -7.004 -5.253 29.57  0 0
    0 0 0 27.58 -7.004
    0 -24.32 0 -7.004 29.57];
p = [0 0 0 0 -45]'; % kN

si determinino tutti gli scostamenti \(u_i\) per il sistema.

Esercizio

Riprendiamo l’esperimento fatto ne La fattorizzazione con pivoting con la matrice di Hilbert

\[(H_n)_{i,j} = \int_{0}^{1} x^{i+j-2}\,{\rm d}x = \frac{1}{i+j-1}, \qquad i,j=1,\ldots,n\]

che è un noto esempio di matrice malcondizionata. Si scriva un programma specializzato nella risoluzione delle equazioni \(H \mathbf{x} = \mathbf{b}\) mediante il metodo di decomposizione \(LU\) con e senza pivoting, dove \(H\) è la matrice di Hilbert di dimensione arbitraria \(n \times n\), e

\[b_i = \sum_{j=1}^{n} (H_n)_{i,j},\]

L’unico input del programma deve essere \(n\). Eseguendo il programma, si determinino i più grandi valori di \(n\) per i due algoritmi per cui la soluzione è entro 6 cifre significative della soluzione esatta.

Sistema di molle

Fig. 8 Sistema lineare di molle.

Esercizio

Il sistema mostrato in Fig. 8 è costituito da \(n\) molle lineari che supportano \(n\) masse. Le costanti elastiche delle molle sono indicate con \(k_i\), i pesi delle masse sono \(W_i\), e \(x_i\) sono gli spostamenti delle masse (misurati dalle posizioni in cui le molle sono indeformate). La cosiddetta formulazione dello spostamento si ottiene scrivendo l’equazione di equilibrio di ciascuna massa e sostituendo le forze elastiche con \(F_i = k_i (x_{i+1} - x_{i} )\).

Il risultato è il sistema simmetrico e tridiagonale di equazioni

\[\begin{split}\begin{split} (k_1 + k_2) x_1 - k_2 x_2 = & W_1,\\ - k_i x_{i-1} + (k_i + k_{i+1}) x_i - k_{i+1} x_{i+1} = & W_i, \qquad i=2,3,\ldots,n-1,\\ -k_n x_{n-1} + k_n x_n = & W_n, \end{split}\end{split}\]

si scriva un programma che risolve queste equazioni per dati valori di \(n\), \(\mathbf{k}\) e \(\mathbf{W}\). Si risolva il problema con \(n=5\) e

\[\begin{split}\begin{split} k_1 = k_2 = k_3 = 10\, \text{N}/\text{mm}\;\; k_4 = k_5 = 5\, \text{N}/\text{mm}\\ W_1 = W_3 = W_5 = 100\,\text{N},\qquad W_2 = W_4 = 50\,\text{N}. \end{split}\end{split}\]

Esercizio: Risoluzione di un Problema Termico con FEM

Con i seguenti passaggi possiamo risolvere un problema termico (stazionario) utilizzando il metodo degli elementi finiti (FEM). Non ci soffermeremo sugli aspetti che pertengono la discretizzazione, ma su quelli che hanno a che fare con la soluzione del sistema lineare risultante.

Passo 1: Creazione del modello Inizia creando un modello FEM per un’analisi termica stazionaria. Utilizza la geometria definita da crackg.

model = femodel(AnalysisType="thermalSteady", ...
  Geometry=@crackg);

Passo 2: Visualizzazione della geometria Visualizza la geometria del modello con le etichette dei bordi attive.

pdegplot(model,EdgeLabels="on");
axis equal

Passo 3: Definizione delle proprietà del materiale Definisci le proprietà del materiale, inclusa la conducibilità termica, la densità di massa e il calore specifico.

model.MaterialProperties = ...
  materialProperties(ThermalConductivity=1, ...
  MassDensity=1, ...
  SpecificHeat=1);

Passo 4: Impostazione delle condizioni al contorno Imposta le condizioni al contorno per il modello. Ad esempio, imposta la temperatura su un bordo specifico e un carico di calore su un altro.

model.EdgeBC(6) = edgeBC(Temperature=100);
model.EdgeLoad(1) = edgeLoad(Heat=-10);
model.FaceIC = faceIC(Temperature=0);

Passo 5: Generazione della mesh Genera la mesh per il modello.

model = generateMesh(model);

Passo 6: Visualizzazione della mesh Visualizza la mesh generata con elementi triangolari quadratici.

figure
pdemesh(model);
title("Mesh with Quadratic Triangular Elements")
axis equal

Passo 7: Risoluzione del modello Risolvi il modello per ottenere i risultati.

%results = solve(model);
FEMs = assembleFEMatrices(model,"stiff-spring");
us = FEMs.Ks\FEMs.Fs;

Passo 8: Visualizzazione della soluzione Visualizza la soluzione PDE sulla mesh utilizzando i valori ottenuti.

pdeplot(model.Mesh, 'XYData', us, 'ColorMap','hot', 'Mesh', 'on');
colorbar;
xlabel('X-axis');
ylabel('Y-axis');
title('PDE Solution Visualization');
axis equal

Sperimentazione

• Si verifichi la dimensione del sistema lineare con il comando size e si visualizzi la matrice con il comando spy: cosa si osserva?

• Si può modificare la finezza della discretizzazione modificando il comando: model = generateMesh(model,'Hmax',0.02);, dove Hmax rappresenta la dimensione massima del lato dei triangoli; si provi a verificare la variazione del costo al rifinire delle mesh.

• Si faccia sperimentazione usando direttamente il comando \, oppure le varianti del comando lu illustrate nella teoria.

## Bibliografia

 ```{bibliography}
 :filter: docname in docnames