Factorizarea Doolittle#

Înainte de a trece în revistă acest subcapitol, asigurați-vă că ați parcurs cu mare atenție noțiunile introductive.

Explicația algoritmului#

Pornim 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.

Pentru factorizarea Doolittle, impunem suplimentar ca matricea $L$ să fie, pe lângă inferior triunghiulară, unitriunghiulară:

\[ l_{ii}=1,\,\forall i\in\overline{1,n}\,\Rightarrow\, A = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & u_{nn} \end{bmatrix} \]

Înainte să formalizăm calculul valorilor ce alcătuiesc matricele $L$ și $U$, cred că ar fi benefică parcurgerea unui exemplu.

Exemplificarea algoritmului#

Să considerăm matricea $A=\begin{bmatrix}2 & -1 & 3 \\4 & 5 & 1 \\2 & 1 & 2\end{bmatrix}$. Fie $A=LU$ descompunerea LU de tip Doolittle a matricei alese; așadar, se formează ecuația:

\[ \begin{bmatrix} 2 & -1 & 3 \\ 4 & 5 & 1 \\ 2 & 1 & 2 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} \]

Pentru rezolvarea pe hârtie, cea mai simplă modalitate de a ajunge la valorile ce constituie matricele $L$ și $U$ este de a crea un sistem de ecuații liniare în care se execută înmulțirea propriu-zisă:

\[ \begin{bmatrix} 2 & -1 & 3 \\ 4 & 5 & 1 \\ 2 & 1 & 2 \end{bmatrix} = \begin{bmatrix} u_{11} &u_{12} & u_{13} \\ l_{21}u_{11} & l_{21}u_{12} + u_{22} & l_{21}u_{13} + u_{23} \\ l_{31}u_{11} & l_{31}u_{12} + l_{32}u_{22} &l_{31}u_{13} + l_{32}u_{23} + u_{33} \end{bmatrix} \]

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

\[ \begin{bmatrix} 2 & -1 & 3 \\ 4 & 5 & 1 \\ 2 & 1 & 2 \end{bmatrix}=\begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & \frac{2}{7} & 1 \end{bmatrix}\begin{bmatrix} 2 & -1 & 3 \\ 0 & 7 & -5 \\ 0 & 0 & \frac{3}{7} \end{bmatrix} \]

Deducerea algoritmului#

Termenul general

Amintim că, în subcapitolul precedent, am prezentat următoarea formulă pentru termenul general al unei factorizări LU:

\[ \boxed{ a_{ij}=\sum_{k=1}^{\min{\{i,j\}}}l_{ik}u_{kj} }\,,\,\,1\leq i,j\leq n \]

Folosind această ecuație, se cristalizează un algoritm. Să împărțim munca pe mai multe iterații, unde, la fiecare iterație $k\in\overline{1,n}$, ne propunem să calculăm întâi $u_{kk},\dots,u_{kn}$, apoi $l_{k+1,n},\dots,l_{nk}$.

De ce așa? Să urmărim pentru un moment descompunerea din exemplul anterior:

\[ \begin{bmatrix} u_{11} &u_{12} & u_{13} \\ l_{21}u_{11} & l_{21}u_{12} + u_{22} & l_{21}u_{13} + u_{23} \\ l_{31}u_{11} & l_{31}u_{12} + l_{32}u_{22} &l_{31}u_{13} + l_{32}u_{23} + u_{33} \end{bmatrix} \]

Observăm cum putem să calculăm foarte rapid $u_{11}$, $u_{12}$ și $u_{13}$, întrucât acestea nu depind de nicio altă variabilă.

Putem apoi să ne îndreptăm atenția către $l_{21}$ și $l_{31}$, care necesita cunoștințe despre $u_{11}$ și $u_{12}$, dobândite la pasul anterior. Calcularea lui $u_{22}$ și a lui $u_{23}$ este acum posibilă.

Într-un final, putem calcula și $l_{32}$, urmat de $u_{33}$. Ordinea finală în care pot fi obținute elementele este descrisă mai jos:

\[ u_{11},u_{12},u_{13}\rightarrow l_{21}, l_{31}\rightarrow u_{22},u_{23}\rightarrow l_{32}\rightarrow u_{33} \]

Pentru cei ce au fost atenți, cu siguranță ați realizat că aceata nu este singura ordine a calculelor ce se poate utiliza pentru a ajunge la un rezultat corect - de exemplu, se poate să calculăm $l_{31}$ înaintea lui $u_{12}$.

Generalizarea propusă însă de algoritm are grijă ca respectarea unui set cât mai simplu de reguli să conducă la un rezultat bun.

Pentru a aprofunda algoritmul, voi face o pseudo-notație menită să vă ajute să vizualizați mai bine algoritmul: matricea pe care o voi scrie acum va conține elemente de forma "Pas Xy: element", unde X este iterația ($i$-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ă:

\[ \begin{bmatrix} \,\,\color{#000000}\colorbox{#c0e9ef}{\textbf{Pas 1a: }$u_{11}$} & \color{#000000}\colorbox{#c0e9ef}{\textbf{Pas 1a: }$u_{12}$} & \color{#000000}\colorbox{#c0e9ef}{\textbf{Pas 1a: }$u_{13}$}\,\,\\ \,\,\color{#000000}\colorbox{#c4d8f3}{\textbf{Pas 1b: }$l_{21}$} & \color{#000000}\colorbox{#c8c7f7}{\textbf{Pas 2a: }$u_{22}$} & \color{#000000}\colorbox{#c8c7f7}{\textbf{Pas 2a: }$u_{23}$}\,\,\\ \,\,\color{#000000}\colorbox{#c4d8f3}{\textbf{Pas 1b: }$l_{31}$} & \color{#000000}\colorbox{#ccb5fb}{\textbf{Pas 2b: }$l_{32}$} & \color{#000000}\colorbox{#d0a4ff}{\textbf{Pas 3a: }$u_{33}$}\,\, \end{bmatrix} \]

Așadar: linie, coloană, linie, coloană, ...

Sintetizarea algoritmului#

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

Se calculează mai întâi valorile $u_{kj}$, $j\in\overline{k,n}$, pornind de la formula termenului general:
\[ \boxed{u_{kj}=a_{kj}-\sum_{t=1}^{k-1}l_{kt} u_{tj}} \]
Utilizând aceste valori, se vor calcula $l_{ik}$, $i\in\overline{k+1,n}$, pornind tot de la formula termenului general:
\[ \boxed{l_{ik}=\frac{1}{u_{kk}}\left(a_{ik}-\sum_{t=1}^{k-1}l_{it} u_{tk}\right)} \]

Demonstrarea formulelor

Aceste formule nu au fost scoase de nicăieri - se poate dovedi corectitudinea matematică a acestora, însă este destul de lungă demonstrația.

Pentru cei cărora le face plăcere, acesta poate fi un exercițiu bun.

Cod ilustrativ#

Implementarea Doolittle (sub forma funcției Doolittle) poate fi acum ușor realizată. 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 unitriunghiulară;
U - o matrice superior triunghiulară.

Fiind o factorizare LU, se respectă proprietatea că $A=LU$.

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

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

Licență#

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