Factorizarea Crout#

Această factorizare este asemnănătoare cu Doolittle - din același motiv, ar fi de preferat să fi parcurs cu atenție secțiunea introductivă.

Explicația algoritmului#

Vom începe tot cu o matrice pătratică \(A\in\mathbb{R}^{n\times n}\) cu \(n\) linii și \(n\) coloane, \(n\in\mathbb{N}^*\), a cărei factorizare LU, \(A=LU\), vrem să o calculăm.

Cerând suplimentar ca matricea \(U\) să fie, pe lângă superior triunghiulară, unitriunghiulară, vom obține factorizarea Crout:

Asemănător factorizării Doolittle, prin forțarea setării valorilor aflate pe diagonala principală a matricei \(U\), se pot realiza calculele aferente celorlalte valori. Acest rezultat va fi demonstrat printr-un exemplu concret.

Exemplificarea algoritmului#

Pentru a ilustra diferența între cele două metode prezentate până acum, vom factoriza aceeași matrice \(A=\begin{bmatrix}2 & -1 & 3 \\4 & 5 & 1 \\2 & 1 & 2\end{bmatrix}\):

Rezolvarea pe hârtie presupune executarea înmulțirii propriu-zise a matricelor \(L\) și \(U\) pentru a se ajunge la un sistem asemănător acestuia:

Rezolvarea acestui sistem are ca efect soluționarea tuturor valorilor necesare; sărind peste calcul, se obține următoarea scriere Crout a matricei \(A\):

Pentru completitudine, amintim că factorizarea Doolittle a obținut matricele \(L_D=\begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & \frac{2}{7} & 1 \end{bmatrix}\) și \(U_D=\begin{bmatrix} 2 & -1 & 3 \\ 0 & 7 & -5 \\ 0 & 0 & \frac{3}{7} \end{bmatrix}\).

Deducerea algoritmului#

Folosind această ecuație, precum și intuiția dobândită în subcapitolul precedent, suntem inspirați să creăm următorul algoritm: împărțim munca pe mai multe iterații, unde, la fiecare iterație \(j\in\overline{1,n}\), ne propunem să calculăm întâi \(l_{jj},\dots,l_{nj}\), apoi \(u_{j,j+1},\dots,u_{jn}\).

De ce așa? Analog cu dezbaterea făcută pentru Doolittle, putem privi pentru început descompunerea din exemplul precedent:

Se poate vedea cu ochiul liber că valorile aflate pe prima coloană se pot obține imediat, urmate de valorile de pe prima linie, a doua coloană, a doua linie și, într-un final, a treia coloană/linie:

Dacă ne uităm atent la descompunere, putem constata de unde este optim să scoatem aceste valori. Voi face o pseudo-notație în încercarea de a vă ajuta să vizualizați mai bine algoritmul: matricea pe care o scriu acum va conține elemente de forma "Pas Xy: element", unde X este iterația (\(j\)-ul) la care se calculează elementul element, iar y va fi "a" sau "b", depinzând de ce valoare se calculează, \(u\) ("a") sau \(l\) ("b"). Poziția oricărei notații corespunde cu poziția din care se va calcula elementul în cauză:

Așadar, invers față de Doolittle: coloană, linie, coloană, linie, ...

Sintetizarea algoritmului#

Pe scurt, algoritmul, la un pas \(k\in\overline{1,n}\), poate fi descris astfel:

• Se calculează mai întâi valorile \(l_{ik}\), \(i\in\overline{k,n}\), pornind de la formula termenului general:

  \[ \boxed{l_{ik}=a_{ik}-\sum_{t=1}^{k-1}l_{it}u_{tk}} \]

• Utilizând aceste valori, se vor calcula \(u_{kj}\), \(j\in\overline{k+1,n}\), pornind tot de la formula termenului general:

  \[ \boxed{u_{kj}=\frac{1}{l_{kk}}\left(a_{kj}-\sum_{t=1}^{k-1}l_{kt} u_{tj}\right)} \]

Cod ilustrativ#

Implementarea Crout (sub forma funcției Crout) este simplă, întrucât știm acum toate formulele necesare. Presupunând că notăm cu A matricea pe care încercăm să o factorizăm, vom folosi codul de mai jos.

Funcția returnează perechea [L, U], formată din:

• L - o matrice inferior triunghiulară;

• U - o matrice superior unitriunghiulară.

Fiind o factorizare LU, se respectă proprietatea că \(A=LU\).

function [L, U] = Crout(A)
    [n, n] = size(A);
    L = eye(n);
    U = zeros(n);

    for k = [1 : n]
        % Partea I - calcularea valorilor l
        for i = [k : n]
            S = 0;
            for t = [1 : k-1]
                S += L(i,t) * U(t,k);
            end
            L(i,k) = A(i,k) - S;
        end
        % Partea II - calcularea valorilor u
        for j = [k+1 : n]
            S = 0;
            for t = [k+1 : n]
                S += L(k,t) * U(t,j);        
            end
            U(i,k) = (A(k,j) - S) / L(k,k);
        end
    end
end

Licență#

The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0