Metoda lui Jacobi#

Așa cum am menționat anterior, diferența între metodele iterative este dictată de alegerea valorilor pentru matricele \(N\) și \(P\) ce îl descompun pe \(A=D-L-U\). Metoda lui Jacobi folosește următoarele valori:

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 lui Jacobi poate fi scris în două forme diferite, dar echivalente:

1. Forma matriceală. Folosind relațiile anterioare, deducem că:

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

2. Forma analitică. Pornind de la relația \(\bm{x^{(p+1)}}=G_\text{Jacobi}\bm{x^{(p)}}+\bm{c_\text{Jacobi}}\), obținem:

   \[ \boxed{x_i^{(p+1)}=\frac{b_i-\sum_{j=1,j\neq i}^na_{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{Jacobi}\bm{x^{(p)}}+\bm{c_\text{Jacobi}}\):

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

   Ținând cont de faptul că \(A=D-L-U\) și de faptul că \(D\) este o matrice diagonală (și deci \(D^{-1}\) conține exclusiv inversul elementelor de pe diagonala principală a matricei \(D\)), se obține:

   \[\begin{split} \begin{bmatrix} x_1^{(p+1)}\\ x_2^{(p+1)}\\ \vdots\\ x_n^{(p+1)} \end{bmatrix}&= \begin{bmatrix} a_{11}^{-1}&0&\dots&0\\ 0&a_{22}^{-1}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&a_{nn}^{-1}\\ \end{bmatrix} \left( \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}-\begin{bmatrix} 0&a_{12}&\dots&a_{1n}\\ a_{21}&0&\dots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&\dots&\\ \end{bmatrix}\begin{bmatrix} x_0^{(p)}\\ x_1^{(p)}\\ \vdots\\ x_n^{(p)} \end{bmatrix} \right)\\ &= \begin{bmatrix} a_{11}^{-1}&0&\dots&0\\ 0&a_{22}^{-1}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&a_{nn}^{-1}\\ \end{bmatrix}\begin{bmatrix} b_0-\sum_{j=1,j\neq 1}^na_{ij}x_j^{(p)}\\ b_1-\sum_{j=1,j\neq 2}^na_{ij}x_j^{(p)}\\ \vdots\\ b_n-\sum_{j=1,j\neq n}^na_{ij}x_j^{(p)} \end{bmatrix} \end{split}\]

   Se cristalizează deci forma analitică:

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

Exemplificarea algoritmului#

Să considerăm următorul sistem de ecuații liniare, scris direct sub formă matriceală: \(\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}\) și rulăm algoritmul lui Jacobi pentru \(6\) iterații:

• Iterația 1. Vom folosi direct formula deja cunoscută (în forma analitică) a metodei Jacobi pentru a obține valorile folosite mai departe, anume:

  \[ \begin{split} x_1^{(1)}&=\left[2-(1\cdot0+1\cdot0)\right]\\ x_2^{(1)}&=\frac{1}{6}\left[9-(-2\cdot0+1\cdot0)\right]\\ x_3^{(1)}&=\frac{1}{7}\left[-6-(-1\cdot0+1\cdot0)\right]\\ \end{split}\Rightarrow \begin{cases} x_1^{(1)}=2\\ x_2^{(1)}=1.5\\ x_3^{(1)}\approx-0.8571 \end{cases} \]

  Aceste rezultate vor fi ținute minte pentru următoarea iterație: \(\bm{x^{(1)}}=\begin{bmatrix} 2\\1.5\\-0.8571 \end{bmatrix}\).

• Iterația 2. Rezultatele precedente vor fi folosite din nou în forma analitică a metodei:

  \[ \begin{split} x_1^{(2)}&=\left[2-(1\cdot1.5-1\cdot0.8571)\right]\\ x_2^{(1)}&=\frac{1}{6}\left[9-(-2\cdot2-1\cdot0.8571)\right]\\ x_3^{(1)}&=\frac{1}{7}\left[-6-(-1\cdot2+1\cdot1.5)\right]\\ \end{split}\Rightarrow \begin{cases} x_1^{(2)}\approx1.3571\\ x_2^{(2)}\approx2.3095\\ x_3^{(2)}\approx-0.7857 \end{cases} \]

  Aceste rezultate vor fi ținute minte pentru următoarea iterație: \(\bm{x^{(2)}}=\begin{bmatrix} 1.3571\\2.3095\\-0.7857 \end{bmatrix}\).

• Iterațiile 3-6. Pentru a nu plictisi cititorul cu calcule, vom continua să ilustrăm evoluția vectorului \(\bm{x}\) pe parcursul iterațiilor:

  \[ \bm{x^{(3)}}=\begin{bmatrix} 0.4762\\2.0833\\-0.9932 \end{bmatrix}\Rightarrow \bm{x^{(4)}}=\begin{bmatrix} 0.9099\\1.8243\\-1.0867 \end{bmatrix}\Rightarrow \bm{x^{(5)}}=\begin{bmatrix} 1.2625\\1.9844\\-0.9878 \end{bmatrix}\Rightarrow \bm{x^{(6)}}=\begin{bmatrix} 1.0034\\2.0855\\-0.9603 \end{bmatrix} \]

Menționăm că soluția corectă era \(\bm{x}=\begin{bmatrix} 1\\2\\-1 \end{bmatrix}\), deci observăm cum algoritmul a făcut o aproximare pertinentă în doar \(6\) iterații.

Cod ilustrativ#

Implementarea metodei lui Jacobi sub forma funcției jacobi, va utiliza 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 = jacobi(A, b, x0, tol, max_iter)
    n = length(A);
    x = x0;

    for k = [1 : max_iter]
        for i = [1 : n]
            val_sum = (A(i,:) * x0(:)) - (A(i,i) * x0(i));
            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