Metoda Gauss-Seidel#

Dacă ar exista un algoritm prin care să folosim o parte din rezultatele pe care le calculăm în iterația curentă, se poate demonstra matematic că rata de convergență va crește. Acest algoritm există și se numește Metoda Gauss-Seidel și folosește:

Algoritmul explicit#

Pornim de la sistemul de ecuații liniare \(A\bm{x}=\bm{b}\), unde \(A\in\mathbb{R}^{n\times n}\), \(\bm{x},\bm{b}\in\mathbb{R}^n\), \(n\in\mathbb{N}^*\). Algoritmul Gauss-Seidel poate fi scris tot în două forme diferite, dar echivalente:

• Forma matriceală. În formă matriceală, metoda Gauss-Seidel se reduce la:

  \[ \begin{cases} N=D-L\\ P=U \end{cases}\Rightarrow \boxed{\begin{cases} G_\text{Gauss-Seidel}=N^{-1}P=(D-L)^{-1}U\\ \bm{c_\text{Gauss-Seidel}}=N^{-1}\bm{b}=(D-L)^{-1}\bm{b} \end{cases}} \]

• Forma analitică. Din punct de vedere al calculelor, acestea se pot desfășura pentru a se obține următoarea recurență:

  \[ \boxed{x_i^{(p+1)}=\frac{b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(p+1)}-\sum_{j=i+1}^{n}a_{ij}x_j^{(p)}}{a_{ii}}}\,,\,\,\forall i\in\overline{1,n} \]
  Demonstrație

  Să desfășurăm ecuația \(\bm{x^{(p+1)}}=G_\text{GS}\bm{x^{(p)}}+\bm{c_\text{GS}}\):

  \[\begin{split} \bm{x^{(p+1)}}&=(D-L)^{-1}U\bm{x^{(p)}}+(D-L)^{-1}\bm{b}\\ (D-L)\bm{x^{(p+1)}}&=U\bm{x^{(p)}}+\bm{b} \end{split}\]

  Ținând cont de faptul că \(A=D-L-U\), se obține:

  \[\begin{split} \begin{bmatrix} a_{11}&0&\dots&0\\ a_{21}&a_{22}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&\dots&a_{nn}\\ \end{bmatrix} \begin{bmatrix} x_1^{(p+1)}\\ x_2^{(p+1)}\\ \vdots\\ x_n^{(p+1)} \end{bmatrix}&= \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}-\begin{bmatrix} 0&a_{12}&\dots&a_{1n}\\ 0&0&\dots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&0\\ \end{bmatrix}\begin{bmatrix} x_1^{(p)}\\ x_2^{(p)}\\ \vdots\\ x_n^{(p)} \end{bmatrix} \end{split}\]

  Așadar, se obține un sistem inferior triunghiular ale cărei ecuații se rezumă la:

  \[\forall i\in\overline{1,n}:\,\,\sum_{j=1}^ia_{ij}x_i^{(p+1)}=b_i-\sum_{j=i+1}^na_{ij}x_j^{(p)}\]

  Dacă folosim SIT (\(b_i\rightarrow b_i-\sum_{j=i+1}^na_{ij}x_j^{(p)}\)), obținem chiar:

  \[ \boxed{x_i^{(p+1)}=\frac{b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(p+1)}-\sum_{j=i+1}^{n}a_{ij}x_j^{(p)}}{a_{ii}}}\,,\,\,\forall i\in\overline{1,n} \]

Exemplificarea algoritmului#

Pentru a exemplifica de ce acest algoritm este mai bun, vom folosi același sistem: \(\begin{bmatrix} 1 & 1 & 1\\ -2 & 6 & 1\\ -1 & 1 & 7 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}=\begin{bmatrix} 2 \\ 9 \\ -6 \end{bmatrix}\). Considerăm totodată soluția inițială \(\bm{x}^{(0)}=\begin{bmatrix} 0\\0\\0 \end{bmatrix}\). Vom ilustra Gauss-Seidel în \(4\) iterații:

• Iterația 1. Folosim direct formula deja cunoscută (în forma analitică) a metodei Gauss-Seidel pentru a obține valorile ce vor fi utilizate mai departe, și anume:

  \[ \begin{split} x_1^{(1)}&=\left[2-(1\cdot0+1\cdot0)\right]=2\\ x_2^{(1)}&=\frac{1}{6}\left[9-(-2\cdot2)-(1\cdot0)\right]\approx2.1667\\ x_3^{(1)}&=\frac{1}{7}\left[-6-(-1\cdot2+1\cdot2.1667)\right]\approx-0.8810\\ \end{split}\Rightarrow \bm{x^{(1)}}= \begin{bmatrix} 2\\2.1667\\-0.8810 \end{bmatrix} \]

• Iterația 2. Folosind formula analitică, se obține:

  \[ \begin{split} x_1^{(2)}&=\left[2-(1\cdot2.1667-1\cdot0.8810)\right]\approx0.7143\\ x_2^{(2)}&=\frac{1}{6}\left[9-(-2\cdot0.7143)+(1\cdot0.8810)\right]\approx1.8849\\ x_3^{(2)}&=\frac{1}{7}\left[-6-(-1\cdot0.7143+1\cdot1.8849)\right]\approx-1.0244\\ \end{split}\Rightarrow \bm{x^{(2)}}= \begin{bmatrix} 0.7143\\1.8849\\-1.0244 \end{bmatrix} \]

• Iterațiile 3-4. Vom trece repede prin rezultate, fără a mai ilustra calculele, acestea fiind făcute în aceeași manieră:

  \[ \bm{x^{(3)}}=\begin{bmatrix} 1.1395\\2.0505\\-0.9873 \end{bmatrix}\Rightarrow \bm{x^{(4)}}=\begin{bmatrix} 0.9368\\1.9768\\-1.0057 \end{bmatrix} \]

Amintim cititorilor că soluția corectă era \(\bm{x}=\begin{bmatrix} 1\\2\\-1 \end{bmatrix}\), iar aproximarea Jacobi după \(4\) iterații obținea \(\bm{x^{(4)}_{\text{Jacobi}}}=\begin{bmatrix} 0.9099\\1.8243\\-1.0867 \end{bmatrix}\). Se observă deci cum rata de convergență a fost îmbunătățită.

Convergență#

O matrice de coeficienți care este dominant diagonală garantează convergența metodei Gauss-Seidel. Propoziția inversă nu este neapărat adevărată, adică convergența metodei nu implică o matrice diagonal dominantă.

Cod ilustrativ#

Codul de mai jos va implementa algoritmul Gauss-Seidel sub forma funcției gauss_seidel, presupunând notațiile:

• A - matricea coeficienților sistemului;

• b - vectorul coloană din partea dreaptă a egalității \(A\bm{x}=\bm{b}\);

• x0 - aproximația inițială pentru soluția sistemului, echivalentul lui \(x_0\);

• tol - toleranța acceptată, adică diferența minimă dintre soluțiile a două iterații consecutive pentru ca algoritmul să continue;

• max_iter - numărul maxim de iterații pe care îl poate executa algoritmul înainte să renunțe.

Funcția returnează x, vectorul coloană soluție a sistemului de ecuații liniare.

function x = gauss_seidel(A, b, x0, tol, max_iter)
    n = length(A);
    x = x0;

    for k = [1 : max_iter]
        for i = [1 : n]
            val_sum = (A(i,1:i-1) * x(1:i-1)) + (A(i,i+1:n) * x0(i+1:n));
            x(i) = (b(i) - val_sum) / A(i,i);
        end
        if norm(x - x0) < tol
            break
        end
        x0 = x;
    end
end

Licență#

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