Álgebra Linear Computacional

Aula 20: Métodos Iterativos Estacionários (Jacobi e Gauss-Seidel)

Heitor S. Ramos
ramosh@dcc.ufmg.br

Créditos

Important

Esses slides são fortemente baseados no livro do prof. Frederico Ferreira Campos, Filho.

Métodos iterativos estacionários

  • Gerar a partir de \(x^0\) uma sequência de vetores \[ \{x^1,x^2, \ldots, x^n\}\longrightarrow x \]

  • Uma mesma série de operações é repetida várias vezes

  • Existem várias classes de métodos iterativos

  • Seja \(M\) a matriz de iteração e \(c\) um vetor constante \[ x^{k+1} = Mx^k +c \]

  • Método iterativo é dito estacionário quando a matriz \(M\) for fixa

  • Exemplos de métodos iterativos: Jacobi, Gauss-Seidel e sobre-relaxação sucessiva

Condição de convergência

Teorema (Condição necessária):

O método iterativo converge com qualquer valor inicial \(x^0\) se, e somente se, \(\rho(M)<1\), sendo \(\rho(M)\) o raio espectral (maior autovalor em módulo) da matriz de iteração \(M\)

Teorema (Condição suficiente):

É condição suficiente para a convergência dos métodos iterativos de Jacobi e Gauss-Seidel que a matriz dos coeficientes A seja diagonal estritamente dominante, ou seja,

\[ \vert a_{ii}\vert > \sum_{\substack{j=1\\j\neq i}}^{n}\vert a_{ij}\vert, i = 1,2,\ldots\,n \]

  • A convergência não depende do valor inicial \(x^0\)

Critério de parada

  • A cada passo do método iterativo a solução é obtida com exatidão crescente

\[ \lim_{k\rightarrow \infty} x^k =x \]

  • Processo deve ser iterrompido quando algum critério de parada for satisfeito

\[ \frac{\Vert x^k - x^{k-1}\Vert}{\Vert x^k \Vert}< \epsilon \] ou \[ k \ge k_{\max} \] - \(\epsilon:\) tolerância - \(k_{\max}:\) número máximo de iterações

Critério de parada adotado

  • Usando norma-\(\infty\)

\[ \frac{{\max\limits_{i\le i \le n}}\vert x^k_i - x^{k-1}_i\vert }{\max\limits_{i\le i \le n} \vert x^k_i\vert}\le \epsilon, \] - \(x^k_i:\) \(i\)-ésimo componente do vetor \(x^k\) obtido na \(k\)-ésima iteração - A tolerância \(\epsilon\) define com qual exatidão a solução é calculada - Em aritimética de ponto flutuante a exatidão não pode ser tão grand quanto queira - Ela é limitada de acordo com o número de bytes das variáveis

Método de Jacobi

  • Decompor a matriz \(A\) de modo que \[ A = D +E + F \]

  • \(D\) é uma matriz diagonal e \(E\) e \(F\) são matrizes triangulares inferiores e superiores, respectivamente, com diagonais nulas

  • O sistema \(Ax + b\) é escrito da forma \[ (D+E+F)x = b\\ Dx = -(E+F)x +b \]

  • Que é convertida para o método iterativo

\[ x^{k+1} = (-D^{-1}(E+F))x^k + D^{-1}b\\ x^{k+1} = Jx^k+c \]

  • Matriz de iteração do método de Jacobi: \(J = -D^{-1}(E+F)\)

Forma análoga de dedução do método de Jacobi

  • Sistema de equações lineares na forma \[ a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \ldots+ a_{1n}x_n = b_1\\ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + \ldots+ a_{2n}x_n = b_2\\ a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + \ldots+ a_{3n}x_n = b_3\\ \vdots\\ a_{n1}x_1 + a_{n2}x_2 + a_{n3}x_3 + \ldots+ a_{nn}x_n = b_n\\ \]
  • Explicitar \(x_i\) na \(i\)-ésima equação
  • Equações de iterações do método de Jacobi \[ \begin{cases} x^{k+1}_1 = \frac{1}{a_{11}}\left(-a_{12}x^k_2 -a_{13}x^k_3 - \ldots -a_{1n}x^k_n + b_1\right), \\ x^{k+1}_2 = \frac{1}{a_{22}}\left(-a_{21}x^k_1 -a_{23}x^k_3 - \ldots -a_{2n}x^k_n + b_2\right), \\ x^{k+1}_3 = \frac{1}{a_{33}}\left(-a_{31}x^k_1 - a_{32}x^k_2 - \ldots -a_{3n}x^k_n + b_3\right), \\ \vdots \\ x^{k+1}_n = \frac{1}{a_{nn}}\left(-a_{n1}x^k_1 - a_{n2}x^k_2 - \ldots -a_{n,n-1}x^k_{n-1} + b_n\right) \end{cases} \]

Forma de recorrência \(x^{k+1} = Jx^k + c\)

\[ \underbrace{ \begin{bmatrix} x^{k+1}_1 \\ x^{k+1}_2 \\ x^{k+1}_3 \\ \vdots \\ x^{k+1}_n \end{bmatrix}}_{x^{k+1}} = \underbrace{ \begin{bmatrix} 0 & -\frac{a_{12}}{a_{11}} & -\frac{a_{13}}{a_{11}} & \ldots & -\frac{a_{1n}}{a_{11}}\\ -\frac{a_{21}}{a_{22}} & 0 & \frac{-a_{23}}{a_{22}} & \ldots & \frac{-a_{2n}}{a_{22}}\\ -\frac{a_{31}}{a_{33}} & -\frac{-a_{32}}{a_{33}}& 0 & \ldots & -\frac{-a_{3n}}{a_{33}}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\frac{a_{n1}}{a_{nn}} & -\frac{-a_{n2}}{a_{nn}}& \ldots & \frac{-a_{n,n-1}}{a_{nn}} & 0\\ \end{bmatrix}}_J \underbrace{ \begin{bmatrix} x^{k}_1 \\ x^{k}_2 \\ x^{k}_3 \\ \vdots \\ x^{k}_n \end{bmatrix}}_{x^{k}} + \underbrace{ \begin{bmatrix} \frac{b_1}{a_{11}} \\ \frac{b_2}{a_{22}} \\ \frac{b_3}{a_{33}} \\ \vdots \\ \frac{b_n}{a_{nn}} \end{bmatrix}}_{c} \]

  • Convergência independe do valor inicial \(x^0\)
  • Vetor inicial: \[ x^0_i = \frac{b_i}{a_{ii}} \]

Exemplo do método de Jacobi

Resolver o sistema de equações pelo método de Jacobi com \(\epsilon<10^{-5}\) e \(k_{\max}=50\) usando os critérios discutidos anteriormente \[ \begin{bmatrix} 10 & 3 & -2 \\ 2 & 8 & -1 \\ 1 & 1 & 5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 57 \\ 20 \\ -4 \end{bmatrix} \]

  • Matrix diagonal dos coeficientes estritamente dominante \[ \vert 10\vert > \vert 3 \vert + \vert -2\vert \\ \vert 8 \vert > \vert 2 \vert + \vert -1\vert \\ \vert 5 \vert > \vert 1 \vert + \vert 1\vert \\ \]

  • Equações de iterações

\[ x^{k+1}_1 = \frac{1}{10}\left(-3x^k_2 + 2x^k_3 + 57 \right), \\ x^{k+1}_2 = \frac{1}{8}\left(-2x^k_1 + x^k_3 + 20 \right), \\ x^{k+1}_3 = \frac{1}{5}\left(-x^k_1 - x^k_2 - 4 \right) \]

  • Vetor inicial \(x^0 = \begin{bmatrix}5,7 & 2,5 & -0,8\end{bmatrix}^\top\)

\[ x^{1}_1 = \frac{1}{10}\left(-3x^0_2 + 2x^0_3 + 57 \right)= \frac{1}{10}\left(-3(2,5) + 2(-0,8) + 57 \right) = 4,79, \\ x^{1}_2 = \frac{1}{8}\left(-2x^0_1 + x^0_3 + 20 \right)= \frac{1}{8}\left(-2(5,7) +(-0,8) + 20 \right) = 0,975, \\ x^{1}_3 = \frac{1}{5}\left(-x^0_1 - x^0_2 - 4 \right)= \frac{1}{5}\left(-(5,7) - (2,5) - 4 \right) = -2,44 \]

  • Vetor da primeira iteração: \(x^1 = \begin{bmatrix}4,79 & 0,975& -2,44\end{bmatrix}^\top\)

  • Critério de parada \[ \frac{\Vert x^1 - x^0\Vert_\infty}{\Vert x^1 \Vert_\infty} = \frac{\max(\vert 4,79 - 5,7\vert , \vert 0,975 - 2,5\vert, -2,44 - (-0,8)\vert)}{\max(\vert 4,79\vert, \vert 0,975\vert, \vert -2,44\vert)}\\ \frac{\Vert x^1 - x^0\Vert_\infty}{\Vert x^1 \Vert_\infty} = 0,3424 \]

  • Se repitirmos esse procedimento veremos \(x \approx x^9 = \begin{bmatrix}5,00000 & 1,00001 & -2,000000\end{bmatrix}^\top\)

Exemplo do método de Jacobi (2)

Calcular três aproximações do vetor solução do sistema abaixo usando a formulação \[ x^{k+1} = -D^{-1}(E+F)x^k + D^{-1}b = Jx^k +c \]

\[ \begin{bmatrix} 10 & 3 & -2 \\ 2 & 8 & -1 \\ 1 & 1 & 5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 57 \\ 20 \\ -4 \end{bmatrix} \]

Decompondo \(A = D+E+F\) temos \[ D = \begin{bmatrix} 10 & 0 & 0\\ 0 & 8 & 0 \\ 0 & 0 & 5 \end{bmatrix}, E = \begin{bmatrix} 0 & 0 & 0\\ 2 & 0 & 0 \\ 1 & 1 & 0 \end{bmatrix}, F = \begin{bmatrix} 0 & 3 & -2\\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{bmatrix} \]

  • Matrix \(J\) e vetores \(c\) e \(x^0\)

\[ J = -D^{-1}(E+F) = \begin{bmatrix} 0 & -0,3 & 0,2 \\ -0,25 & 0 & 0,125 \\ -0,2 & -0,2 & 0 \end{bmatrix} \] e \[ c = D^{-1}b = x^0 = \begin{bmatrix} 5,7 & 2,5 & -0,8 \end{bmatrix}^\top \]

k \(x^k_1\) \(x^k_2\) \(x^k_3\)
0 5,70000 2,50000 -0,80000
1 4,79000 0,97500 -2,44000
2 4,91950 0,9975 -1,95300
3 5,01015 1,02600 -1,98340

Método de Gauss-Seidel

  • Decompor a matriz \(A\) de forma que \[ A = D + E + F \]

  • \(D\) é uma matriz diagonal e \(E\) e \(F\) são matrizes triangulares inferior e superior, respectivamente, com diagonais nulas

  • Sistema linear \(Ax=b\) escrito na forma \[ (D+E+F)x = b\\ (D+E)x = -F x +b \]

  • Forma de iteração obtida pela recorrência \[ x^{k+1} = \left(-(D+E)^{-1}\right)x^k + (D+E)^{-1}b\\ x^{k+1} = Sx^k + d \]

  • Matriz de iteração do método de Gauss-Seidel: \(S = -(D+E)^{-1}F\)

Forma alternativa de dedução do método de Gauss-Seidel

  • Seja \[ (D + E + F)x = b \longrightarrow (D + E)x = −Fx + b \]

  • Na forma de recorrência \[ (D+E)x^{k+1} = −Fx^k + b \longrightarrow Dx^{k+1} =−Ex^{k+1} −Fx^k + b \]

  • Escrevendo a segunda equação na forma matricial

\[ Dx^{k+1} = - \underbrace{ \begin{bmatrix} 0 & 0 & 0 & \ldots & 0 \\ a_{21} & 0 & 0 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n-1,1} & a_{n-1,2}& \ldots & 0 & 0 \\ a_{n1} & a_{n2} & \ldots & a_{n,n-1} & 0 \end{bmatrix}}_{E} \underbrace{ \begin{bmatrix} x^{k+1}_1\\ x^{k+1}_2\\ x^{k+1}_3\\ \vdots \\ x^{k+1}_n \end{bmatrix}}_{x^{k+1}} - \\ \underbrace{ \begin{bmatrix} 0 & a_{12} & a_{13} & \ldots & a_{1n} \\ 0 & 0 & a_{23} & \ldots & a_{2n} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \ldots & 0 & a_{n-1,n} \\ 0 & 0 & \ldots & 0 & 0 \\ \end{bmatrix}}_{F} \underbrace{ \begin{bmatrix} x^{k}_1\\ x^{k}_2\\ x^{k}_3\\ \vdots \\ x^{k}_n \end{bmatrix}}_{x^k} + \underbrace{ \begin{bmatrix} b_1\\ b_2\\ b_3\\ \vdots \\ b_n \end{bmatrix}}_{B} \]

Equações de iterações do método Gauss-Seidel

\[\begin{cases} x^{k+1}_1 = \frac{1}{a_{11}}\left(-a_{12}x^k_2 - a_{13}x^k_3 - \ldots - a_{1n}x^k_n+b_1\right) \\ x^{k+1}_2 = \frac{1}{a_{22}}\left(-a_{21}x^{k+1}_1 - a_{23}x^k_3 - \ldots - a_{2n}x^k_n+b_2\right) \\ x^{k+1}_3 = \frac{1}{a_{33}}\left(-a_{31}x^{k+1}_1 - a_{32}x^{k+1}_2 - \ldots - a_{3n}x^k_n+b_3\right) \\ \vdots \\ x^{k+1}_n = \frac{1}{a_{nn}}\left(-a_{n1}x^{k+1}_1 - a_{n2}x^{k+1}_2 - \ldots - a_{n,n-1}x^{k+1}_{n-1}+b_n\right) \\ \end{cases} \]

  • Vetor \(x^{k+1}\) obtido a partir dos elementos mais recentes, incluindo o próprio \(x^{k+1}\) e \(x^{k}\)

  • Vetor inicial: \[ x^0_i = \frac{b_i}{a_{ii}} \]

Exemplo do método Gauss-Seidel

  • Resolver o sistema a seguir pelo método de Gauss-Seidel com \(\epsilon < 10^{-5}\) e \(k_{\max} = 50\) usando os critérios discutidos anteriormente \[ \begin{bmatrix} 10 & 3 & -2 \\ 2 & 8 & -1 \\ 1 & 1 & 5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 57 \\ 20 \\ -4 \end{bmatrix} \]

  • Matriz dos coeficientes é diagnalmente dominante

  • Equações de iterações: \[ x^{k+1}_1 = \frac{1}{10}\left(-3x^k_2 + 2x^k_3 + 57\right)\\ x^{k+1}_2 = \frac{1}{8}\left(-2x^{k+1}_1 + x^k_3 + 20\right)\\ x^{k+1}_3 = \frac{1}{5}\left(-x^{k+1}_1 - x^{k+1}_2 - 4\right)\\ \]

  • Vetor inicial: \(x^0 = \begin{bmatrix}5,7 & 2,5 & -0,8\end{bmatrix}^\top\)

\[ x^{1}_1 = \frac{1}{10}\left(-3x^0_2 + 2x^0_3 + 57\right) = \frac{1}{10}\left(-3(2,5) + 2(-0,8) + 57\right) = 4,79\\ x^{1}_2 = \frac{1}{8}\left(-2x^{1}_1 + x^0_3 + 20\right) = \frac{1}{8}\left(-2(4,79) + (-0,8) + 20\right) = 1,2025\\ x^{1}_3 = \frac{1}{5}\left(-x^{1}_1 - x^{1}_2 - 4\right) = \frac{1}{5}\left(-(4,79) - (1,2025) - 4\right) = 1,9985\\ \]

  • Vetor da primeira iteração: \(x^1 = \begin{bmatrix}4,79 & 1,2025 & -1,9985\end{bmatrix}^\top\)

  • Critério de parada \[ \frac{\Vert x^1 - x^0\Vert_\infty}{\Vert x^1 \Vert_\infty} = \frac{\max(\vert 4,79 - 5,7\vert , \vert 1,2025 - 2,5\vert, -1,9985 - (-0,8)\vert)}{\max(\vert 4,79\vert, \vert 1,2025\vert, \vert -1,9985\vert)}\\ \frac{\Vert x^1 - x^0\Vert_\infty}{\Vert x^1 \Vert_\infty} = 0,2709 \]

  • Se repitirmos esse procedimento veremos \(x \approx x^6 = \begin{bmatrix}5,00000 & 1,000010& -2,000000\end{bmatrix}^\top\)

Exemplo do método de Gauss-Seidel (2)

Calcular três aproximações do vetor solução do sistema abaixo usando a formulação \[ x^{k+1} = −(D + E)^{−1}Fx^k + (D + E)^{−1}b = Sx^k + d \]

\[ \begin{bmatrix} 10 & 3 & -2 \\ 2 & 8 & -1 \\ 1 & 1 & 5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 57 \\ 20 \\ -4 \end{bmatrix} \]

  • Decompondo \(A = D + E + F\) \[ D = \begin{bmatrix} 10 & 0 & 0\\ 0 & 8 & 0 \\ 0 & 0 & 5 \end{bmatrix}, E = \begin{bmatrix} 0 & 0 & 0\\ 2 & 0 & 0 \\ 1 & 1 & 0 \end{bmatrix}, F = \begin{bmatrix} 0 & 3 & -2\\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{bmatrix} \]

  • Matriz \((D+E)^{-1}\) calculada utilizando o esquema: \((D + E)(D + E)^{−1} = I\)

  • Matrizes \((D+E)^{−1}\) e \(S\)

\[ (D+E)^{−1} = \begin{bmatrix} 0,1 & 0 & 0 \\ -0,025 & 0,125 & 0 \\ -0,015 & -0,025 & 0,2 \end{bmatrix} \] e \[ S=−(D+E)^{−1}F= \begin{bmatrix} 0 & −0,3 & 0,2 \\ 0 & 0,075 & 0,075 \\ 0 & 0,045 & -0,055 \end{bmatrix} \]

  • Vetores \(d\) e \(x^0\)

\[ d = (D + E)^{−1}b = \begin{bmatrix}5,7 & 1,075 &−2,155\end{bmatrix}^\top \] \[ x^0 = D^{−1}b = \begin{bmatrix} 5,7 & 2,5 & −0,8\end{bmatrix}^\top \]

Primeiras aproximações

  • Jacobi
k \(x^k_1\) \(x^k_2\) \(x^k_3\)
0 5,70000 2,50000 -0,80000
1 4,79000 0,97500 -2,44000
2 4,91950 0,9975 -1,95300
3 5,01015 1,02600 -1,98340
  • Gauss-Seidel
k \(x^k_1\) \(x^k_2\) \(x^k_3\)
0 5,70000 2,50000 -0,80000
1 4,79000 1,20250 -1,99850
2 4,93955 1,01530 -1,99097
3 4,99722 1,00182 -1,99981

Método da sobre-relaxação sucessiva

  • Decompor a matriz \(A\) de modo que \[A = D + E + F \]

  • D é uma matriz diagonal e \(E\) e \(F\) são matrizes triangulares inferiores nulas

  • Multiplicando o sistema linear \(Ax=b\) por um parâmetro \(\omega\)

\[ \omega(D+E+F)x = wb \]

  • Somando o vetor nulo \((D-D)x\) ao primeiro termo \[ (D − D)x + \omega (D + E + F )x = \omega b \]

  • Rearranjando \[ (D + \omega E)x = [(1 − \omega)D − \omega F]x + \omega b \]

  • Forma de iteração \[ (D+\omega E)x^{k+1} =[(1−\omega )D− \omega F]x^k +\omega b \] \[ x^{k+1} =(D+\omega E)^{−1}[(1−\omega )D−\omega F]x^k+\omega (D+\omega E)^{−1}b \\ x^{k+1} =Rx^k + e \]

  • Matriz de iteração do método SOR: \(R = (D + \omega E)^{−1} [(1 − \omega )D − \omega F ]\)

  • SOR: Successive over-relaxation

Convergência do método SOR

  • Depende de \(\omega\)

  • Para garantir a converência: \(0<\omega <2\)

  • Usualmente: \(1<\omega <2\)

  • Para \(\omega =1\), a recorrência torna-se \[x^{k+1} = −(D + E)^{−1}Fx^k + (D + E)^{−1}b\]

  • Método de Gauss-Seidel é um caso particular de SOR

Equações de iterações do método SOR

\[\begin{cases} x^{k+1}_1 = \frac{\omega}{a_{11}}\left(-a_{12}x^k_2 - a_{13}x^k_3 - \ldots - a_{1n}x^k_n+ b_1\right) + (1-\omega) x^k_1\\ x^{k+1}_2 = \frac{\omega}{a_{22}}\left(-a_{21}x^{k+1}_1 - a_{23}x^k_3 - \ldots - a_{2n}x^k_n+b_2\right) + (1-\omega) x^k_2\\ x^{k+1}_3 = \frac{\omega}{a_{33}}\left(-a_{31}x^{k+1}_1 - a_{32}x^{k+1}_2 - \ldots - a_{3n}x^k_n+b_3\right) + (1-\omega) x^k_3\\ \vdots \\ x^{k+1}_n = \frac{\omega}{a_{nn}}\left(-a_{n1}x^{k+1}_1 - a_{n2}x^{k+1}_2 - \ldots - a_{n,n-1}x^{k+1}_{n-1}+b_n\right) + (1-\omega) x^k_n\\ \end{cases} \]

Influência de \(\omega\) no raio espectral

Resolver o sistema pelo método SOR, com \(\epsilon < 10^{−5}\) e \(k_{\max} = 500\) usando os critérios estabelecidos

\[ \begin{bmatrix} 4 & -2 & 1 & 3 & 0\\ -1 & 10 & 0 & 8 & 1 \\ -1 & 1 & 15 & 2 & 4 \\ 0 & 1 & 10 & 5 & 1 \\ 2 & -3 & 1 & 2 & 20 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} \begin{bmatrix} 15\\ 56 \\ 74 \\ 57 \\ 107 \end{bmatrix} \]

  • Influência e \(\omega\) no raio espectral \(\rho\) da matriz de iteração \[ R = (D + \omega E)^{−1} [(1 − \omega )D − \omega F] \]
\(\omega\) 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8
\(\rho(R)\) 0,9357 0,8618 0,7747 0,6674 0,5231 0,4459 0,7506 1,0660 1,3943
iterações 118 63 41 29 20 17 44 >500 >500
  • Quanto menor for \(\rho(R)\) menor será o número de iterações

Sobre o teorema da condição suficiente

  • O teorema não se aplica a SOR por conta do parâmetro \(\omega\)
  • No exemplo anterior, SOR não converge para \(\omega > 1,6\) pois \(\rho(R) > 1\)

https://heitorramos.github.io/alc.html