Metoda suprarelaxării#

Cum putem îmbunătăți Gauss-Seidel? Forțat!

Pornim de la ecuația \(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}^*\). Fie \(\omega\in(0,2)\) un număr real, o constantă, introdusă pe cale pur artificială.

Folosind noile notații, obținem:

Această metodă poartă denumirea de Metoda suprarelaxării (prescurtat, SOR) și reprezintă modalitatea optimă de soluționare iterativă a unui sistem de ecuații liniar,

Algoritmul explicit#

Puitem acum să discutăm despre cele două forme ale algoritmului:

1. Forma matriceală. Revenind la \(G\) și la \(\bm{c}\), am obține:

   \[ \begin{cases} G_\text{SOR}=N^{-1}P=(D-\omega L)^{-1}\Big[(1-\omega)D+\omega U\Big]\\ \bm{c_\text{SOR}}=(D-\omega L)^{-1}(\omega\bm{b}) \end{cases} \]
   Legătura GS-SOR

   Dacă impunem \(\omega=1\), obținem Gauss-Seidel.

2. Forma analitică. Din punct de vedere al calculelor, acestea se pot desfășura pentru a se obține următoarea formulă de recurență ce converge cel mai rapid posibil:

   \[ \boxed{x_i^{(p+1)}=(1-\omega)x_{i}^{(p)}+\frac{\omega}{a_{ii}}\left(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(p+1)}-\sum_{j=i+1}^{n}a_{ij}x_j^{(p)}\right)}\,,\,\,\forall i\in\overline{1,n} \]
   Viteza de convergență

   Metoda converge mai rapid decât celelalte două numai în cazul în care valoarea lui \(\omega\) este aleasă corespunzător!

   Demonstrație

   Vom nota \(\overline{\omega}=1-\omega\). Desfășurăm \(\bm{x^{(p+1)}}=G_\text{SOR}\bm{x^{(p)}}+\bm{c_\text{SOR}}\):

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

   Partea stângă denotă un sistem inferior triunghiular:

   \[\begin{split} \begin{bmatrix} a_{11}&0&\dots&0\\ \omega a_{21}&a_{22}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ \omega a_{n1}&\omega a_{n2}&\dots&a_{nn}\\ \end{bmatrix} \begin{bmatrix} x_1^{(p+1)}\\ x_2^{(p+1)}\\ \vdots\\ x_n^{(p+1)} \end{bmatrix}=\dots \end{split}\]

   Partea dreaptă poate fi rescrisă ca:

   \[\begin{split} \dots= \begin{bmatrix} \omega b_1\\ \omega b_2\\ \vdots\\ \omega b_n \end{bmatrix}-\begin{bmatrix} \overline{\omega} a_{11}&\omega a_{12}&\dots&\omega a_{1n}\\ 0&\overline{\omega} a_{22}&\dots&\omega a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&\overline{\omega} a_{nn}\\ \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}:\,\,\omega\sum_{j=1}^{i-1}a_{ij}x_i^{(p+1)}+a_{ii}x_i^{(p+1)}=b_i-\omega\sum_{j=i+1}^na_{ij}x_j^{(p)}-\overline{\omega}a_{ii}x_i^{(p)}\]

   Dacă folosim SIT, obținem chiar:

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

Alegerea termenului introdus#

Nu există o formulă concretă general valabilă pentru calcularea lui \(\omega\). Cu toate acestea, David Young a demonstrat[52] că, în anumite situații, valoarea lui \(\omega\) poate fi calculată în funcție de \(\mu_{\max}\):

unde prin \(\mu_{\max}\) se înțelege valoarea proprie maximală a matricei \(G_{\text{Jacobi}}\) (obținută așa cum a fost descrisă anterior); detalii suplimentare despre ce este aceea o valoare proprie pot fi găsite în capitolul dedicat.

Exemplificarea algoritmului#

Pentru a exemplifica performanțele atinse de algoritm, vom refolosi sistemul: \(\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}\), cu soluția inițială \(\bm{x}^{(0)}=\begin{bmatrix} 0\\0\\0 \end{bmatrix}\) și \(\omega=0.85\). Vom ilustra SOR în \(2\) iterații:

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

  \[ \begin{split} x_1^{(1)}&=0.15\cdot0+0.85\left[2-(1\cdot0+1\cdot0)\right]=1.7\\ x_2^{(1)}&=0.15\cdot0+\frac{0.85}{6}\left[9-(-2\cdot1.7)-(1\cdot0)\right]\approx1.7567\\ x_3^{(1)}&=0.15\cdot0+\frac{0.85}{7}\left[-6-(-1\cdot1.7+1\cdot1.7567)\right]\approx-0.7355\\ \end{split} \] \[\Rightarrow \bm{x^{(1)}}= \begin{bmatrix} 1.7\\1.7567\\-0.7355 \end{bmatrix} \]

• Iterația 2. Folosind rezultatele precedente în forma analitică a metodei, obținem:

  \[ \begin{split} x_1^{(2)}&=0.15\cdot1.7+0.85\left[2-(1\cdot1.7-1\cdot0.7355)\right]\approx1.0870\\ x_2^{(2)}&=0.15\cdot1.7567+\frac{0.85}{6}\left[9-(-2\cdot1.0870)-(1\cdot0.7355)\right]\approx1.9507\\ x_3^{(2)}&=-0.15\cdot0.7355+\frac{0.85}{7}\left[-6-(-1\cdot1.0870+1\cdot1.9507)\right]\approx-0.9438\\ \end{split} \] \[\Rightarrow \bm{x^{(2)}}= \begin{bmatrix} 1.0870\\1.9507\\-0.9438 \end{bmatrix} \]

În doar două iterații, am aproximat soluția corectă, \(\bm{x}=\begin{bmatrix} 1\\2\\-1 \end{bmatrix}\). Amintim că, după două iterații, metoda Jacobi aproxima soluția cu \(\bm{x^{(2)}_{\text{Jacobi}}}=\begin{bmatrix} 1.3571\\2.3095\\-0.7857 \end{bmatrix}\), iar Gauss-Seidel o credea egală cu \(\bm{x^{(2)}_{\text{Gauss-Seidel}}}=\begin{bmatrix} 0.7143\\1.8849\\-1.0244 \end{bmatrix}\).

Cod ilustrativ#

Pentru a implementa metoda suprarelaxării în funcția sor, vom considera:

• 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\);

• omega - constanta introdusă pe cale artificială, echivalentul lui \(\omega\);

• 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 = sor(A, b, x0, omega, 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) = (1-omega) * x0(i) + omega * (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