Algoritmul lui Thomas#

Algoritmul lui Thomas reprezintă de fapt o variantă mai eficientă a algoritmului de eliminare gaussiană ce poate fi folosită pentru matrice tridiagonale. Algoritmul este atât de eficient încât are complexitate liniară!

Intuiție algoritmică#

Să considerăm următorul sistem de ecuații liniare, în care matricea coeficienților este tridiagonală:

Se poate observa că am notat elementele de pe diagonala principală cu \(b_i\), cele aflate deasupra acesteia cu \(c_i\), iar cele aflate sub aceasta cu \(a_i\), unde \(i\) reprezintă indicele liniei pe care se află valoarea în speță.

Algoritmul lui Thomas își propune să rezolve astfel de sisteme, unde \(n\in\mathbb{N}^*\), iar \(a_k,b_k,c_k,d_k,x_k\in\mathbb{R}\) (pentru fiecare valoare validă a lui \(k\)).

Pentru început, să considerăm doar primele două ecuații: \(\begin{cases} b_1x_1+c_1x_2&=d_1\\ a_2x_1+b_2x_2+c_2x_3&=d_2 \end{cases}\). Ca să simplificăm a doua ecuație cât se poate de mult, se poate aplica următoarea transformare ce nu va afecta soluțiile sistemului:

Ecuația rezultantă poate fi scrisă drept \((b_2b_1-c_1a_2)x_2+c_2b_1x_3=d_2b_1-d_1a_2\), ceea ce înseamnă că operația a reușit în a elimina elementul aflat pe prima coloană, imediat sub diagonala principală.

Continuând acest raționament, ne putem îndrepta atenția către a doua, respectiv a treia ecuație (în urma primei operații), adică \(\begin{cases} (b_2b_1-c_1a_2)x_2+c_2b_1x_3&=d_2b_1-d_1a_2\\ a_3x_2+b_3x_3+c_3x_4&=d_3 \end{cases}\) - de data aceasta, se poate aplica următoarea transformare:

Efectul (poate puțin greu de digerat) este următorul:

Așadar, s-a reușit eliminarea lui \(x_2\).

Prin îmbunătățirea treptată a rezultatului, se va obține o matrice tridiagonală ce nu conține valori nenule sub diagonala principală, ceea ce ar însemna un sistem foarte ușor de rezolvat:

Explicația algoritmică#

Ar fi foarte delicată implementarea algoritmului așa cum a fost explicat anterior, însă cu un set de notații ajutătoare, sarcina devine realizabilă.

Fie \(i\in\overline{1,n}\) pasul curent ce va transforma a \(i\)-a linie a matricei \(A\in\mathbb{R}^{n\times n}\) astfel încât să conțină cel mult două intrări nenule, \(n\in\mathbb{N}^*\). Chiar dacă pentru \(i\in\{1,n\}\) nu trebuie făcut nimic, notațiile rămân valabile pentru celelalte variabile.

Fie \(\tilde{a}_i, \tilde{b}_i, \tilde{c}_i\) și \(\tilde{d}_i\) valorile ce se vor găsi în matrice după al \(i\)-lea pas (spre exemplu, \(\tilde{c}_3\) va fi calculat după a treia operație și va deveni constanta \(c_3\) pentru iterațiile \(4\), \(5\) etc.). Aceste valori pot fi calculate utilizând relațiile:

Aceasta este doar o adaptare a algoritmului propus de Thomas, întrucât folosește valori nescalate pentru \(\tilde{b}_i\). Thomas propune normalizarea valorilor \(\tilde{b}_i\) la fiecare pas, ceea ce înseamnă că ecuațiile finale vor fi de forma:

Aceste valori vor descrie următorul sistem de ecuații liniare:

Rezolvat, necunoscutele primesc valori în funcție de aceste formule (evident, algoritmic, considerăm că \(i\) descrește):

Formula anterioară reprezintă o adaptare a soluționării unui SST.

Exemplificarea algoritmului#

Să considerăm sistemul descris matriceal în forma compactă următoare: \(\overline{A}=\begin{bmatrix} 3 & 1 & 0 & | & 5\\ -1 & 3 & -2 & | & -7\\ 0 & 4 & 3 & | & -1 \end{bmatrix}\).

Începem acum să urmăm algoritmul Thomas:

• Pasul 1. Se transformă matricea \(\overline{A}\) astfel încât să aibă sub diagonala principală numai elemente nule, iar pe aceasta numai elemente de \(1\):

  \[ \begin{cases} \tilde{a}_2=\tilde{a}_3=0\\ \tilde{b}_1=\tilde{b}_2=\tilde{b}_3=1\\ \tilde{c}_1=\dfrac{1}{3}\\ \tilde{c}_2=\dfrac{-2}{3+\frac{1}{3}}=-\dfrac{3}{5}\\ \tilde{d}_1=\dfrac{5}{3}\\ \tilde{d}_2=\dfrac{-7+\frac{5}{3}}{3+\frac{1}{3}}=-\dfrac{8}{5}\\ \tilde{d}_3=\dfrac{-1+\frac{32}{5}}{3+\frac{12}{5}}=1\\ \end{cases}\Rightarrow\begin{bmatrix} 1 & 0.(3) & 0\\ 0 & 1 & -0.6\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3 \end{bmatrix}=\begin{bmatrix} 1.(6)\\-1.6\\1 \end{bmatrix} \]

• Pasul 2. Se soluționează sistemul utilizând formula pentru necunoscute: \(\begin{cases} x_3=1\\ x_2=-1.6+0.6=-1\\ x_1=1.(6)+0.(3)=2 \end{cases}\).

Cod ilustrativ#

În implementarea algoritmului Thomas (funcția Thomas), vom folosi următoarele notații:

• a - vectorul de valori aflate sub diagonala principală a matricei coeficienților sistemului;

• b - vectorul de valori aflate pe diagonala principală a matricei coeficienților sistemului;

• c - vectorul de valori aflate deasupra diagonalei principale a matricei coeficienților sistemului;

• d - vectorul coloană din partea dreaptă a egalității.

Funcția returnează x, vectorul ce soluționează ecuația liniară în formă matriceală.

function x = Thomas(a, b, c, d)
    n = length(d);
    a = [0; a];
    c(1) = c(1) / b(1);
    d(1) = d(1) / b(1);
    for i = [2 : n-1]
        temp = b(i) - a(i) * c(i-1);
        c(i) = c(i) / temp;
        d(i) = (d(i) - a(i) * d(i-1)) / temp;
    end
    d(n) = (d(n) - a(n) * d(n-1)) / (b(n) - a(n) * c(n-1));
    x(n) = d(n);
    for i = [n-1 : -1 : 1]
        x(i) = d(i) - c(i) * x(i+1);
    end
    x = x';
end

Concluzie#

Algoritmul este foarte, foarte rapid comparativ cu toate celelalte eliminări gaussiene (\(O(n)\lt \lt O(n^3)\)), însă acesta se poate folosi doar atunci când matricea \(A\) este pătratică și tridiagonală.

Suplimentar, algoritmul este foarte greu de înțeles / reținut, ceea ce îl face mult mai greu de implementat corect. Citirea trebuie de asemenea făcută într-un format codat, pentru că altfel ea se va realiza în \(O(n^2)\) pași.

Licență#

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