Métodos iterativos para sistemas lineares

Cálculo Numérico - Versão Scilab

4.7 Métodos iterativos para sistemas lineares

Na seção anterior, tratamos de métodos diretos para a resolução de sistemas lineares. Em um método direto (por exemplo, solução via fatoração LU) obtemos uma aproximação da solução depois de realizarmos um número finito de operações (só teremos a solução ao final do processo).

Veremos nessa seção dois métodos iterativos básicos para obter uma aproximação para a solução de um sistema linear. Geralmente em um método iterativo, iniciamos com uma aproximação para a solução (que pode ser ruim) e vamos melhorando essa aproximação através de sucessivas iterações.

4.7.1 Método de Jacobi

O método de Jacobi pode ser obtido a partir do sistema linear $\begin{array}{rcl} a_{11} x_{1} + a_{12} x_{2} + \dots + a_{1 n} x_{n} & = & y_{1} & (4.199) \\ a_{21} x_{1} + a_{22} x_{2} + \dots + a_{2 n} x_{n} & = & y_{2} & (4.200) \\ ⋮ & (4.201) \\ a_{n 1} x_{1} + a_{n 2} x_{2} + \dots + a_{n n} x_{n} & = & y_{n} & (4.202) \end{array}$

Isolando o elemento $x_{1}$ da primeira equação temos

x_{1}^{(k + 1)} = \frac{y_{1} - (a_{12} x_{2}^{(k)} + \dots + a_{1 n} x_{n}^{(k)})}{a_{11}}

(4.203)

Note que utilizaremos os elementos $x_{i}^{(k)}$ da iteração $k$ (a direita da equação) para estimar o elemento $x_{1}$ da próxima iteração.

Da mesma forma, isolando o elemento $x_{i}$ de cada equação $i$ , para todo $i = 2, . . ., n$ podemos construir a iteração $\begin{array}{rcl} x_{1}^{(k + 1)} & = & \frac{y_{1} - (a_{12} x_{2}^{(k)} + \dots + a_{1 n} x_{n}^{(k)})}{a_{11}} & (4.204) \\ x_{2}^{(k + 1)} & = & \frac{y_{2} - (a_{21} x_{1}^{(k)} + a_{23} x_{3}^{(k)} + \dots + a_{2 n} x_{n}^{(k)})}{a_{22}} & (4.205) \\ ⋮ & (4.206) \\ x_{n}^{(k + 1)} & = & \frac{y_{2} - (a_{n 1} x_{1}^{(k)} + \dots + a_{n, n - 2} x_{n - 2}^{(k)} + a_{n, n - 1} x_{n - 1}^{(k)})}{a_{n n}} & (4.207) \end{array}$

Em notação mais compacta, o método de Jacobi consiste na iteração $\begin{array}{rcl} x^{(1)} & = & aproximação inicial & (4.208) \\ x_{i}^{(k + 1)} & = & (y_{i} - \sum_{\begin{array}{c} j = 1 \\ j \neq i \end{array}}^{n} a_{i j} x_{j}^{(k)}) ∕ a_{i i} & (4.209) \end{array}$

Exemplo 4.7.1. Resolva o sistema $\begin{array}{rcl} 10 x + y & = & 23 & (4.210) \\ x + 8 y & = & 26 & (4.211) \end{array}$

usando o método de Jacobi iniciando com $x^{(1)} = y^{(1)} = 0$ . $\begin{array}{rcl} x^{(k + 1)} & = & \frac{23 - y^{(k)}}{10} & (4.212) \\ y^{(k + 1)} & = & \frac{26 - x^{(k)}}{8} & (4.213) \\ x^{(2)} & = & \frac{23 - y^{(1)}}{10} = 2, 3 & (4.214) \\ y^{(2)} & = & \frac{26 - x^{(1)}}{8} = 3, 25 & (4.215) \\ x^{(3)} & = & \frac{23 - y^{(2)}}{10} = 1, 975 & (4.216) \\ y^{(3)} & = & \frac{26 - x^{(2)}}{8} = 2, 9625 & (4.217) \end{array}$

Exemplo 4.7.2. Considere o seguinte sistema

\begin{matrix} \begin{aligned} - 3 x_{1} + x_{2} + x_{3} & = 2 \\ 2 x_{1} + 5 x_{2} + x_{3} & = 5 \\ 2 x_{1} + 3 x_{2} + 7 x_{3} & = - 17 \end{aligned} \end{matrix}

(4.218)

Usando o método de Jacobi com aproximação inicial $x^{(1)} = (1, 1, - 1)$ , obtemos os seguintes resultados:


$k$	$x^{(k)}$	$∥ x^{(k)} - x^{(k - 1)} ∥_{\infty}$

1	$(1, 1, 1)$	-x-
2	$(- 0, 67, 0, 80, - 3, 14)$	$2, 1$
3	$(- 1, 45, 1, 90, - 2, 58)$	$1, 1$
4	$(- 0, 90, 2, 10, - 2, 83)$	$5, 5 E - 1$
5	$(- 0, 91, 1, 92, - 3, 07)$	$2, 4 E - 1$
⋮	⋮	⋮
10	$(- 1, 00, 2, 00, - 3, 00)$	$6, 0 E - 3$

Verifique a resposta.

Aqui, cabe um código Scilab explicativo. Escreva você mesmo o código.

Veja como participar da escrita do livro em:

https://github.com/livroscolaborativos/CalculoNumerico

Código Scilab: Método de Jacobi

function [x,deltax]=jacobi(A,b,x,tol,N)
n=size(A,1)
xnew     =x
convergiu=%F                //FALSE;

k=1
while k<=N & ~convergiu
  xnew(1)=(b(1) - A(1,2:n)*x(2:n))/A(1,1)
  for i=2:n-1
    xnew(i)=(b(i) -A(i,1:i-1)*x(1:i-1) ...
                  -A(i,i+1:n)*x(i+1:n) )/A(i,i)
  end
  xnew(n)=  (b(n) -A(n,1:n-1)*x(1:n-1) )/A(n,n)

  deltax=max( abs(x-xnew) )
  if deltax<tol then
     convergiu=%T       //TRUE
  end
  k=k+1
  x=xnew                 // atualiza x
  disp([k,x’,deltax])    // depuracao
end
if ~convergiu then
    error(’Nao convergiu’)
end

endfunction

4.7.2 Método de Gauss-Seidel

Assim, como no método de Jacobi, no método de Gauss-Seidel também isolamos o elemento $x_{i}$ da equação $i$ . Porém perceba que a equação para $x_{2}^{(k + 1)}$ depende de $x_{1}^{(k)}$ na iteração $k$ . Intuitivamente podemos pensar em usar $x_{1}^{(k + 1)}$ que acabou de ser calculado e temos

x_{2}^{(k + 1)} = \frac{y_{2} - (a_{21} x_{1}^{(k + 1)} + a_{23} x_{3}^{(k)} + \dots + a_{2 n} x_{n}^{(k)})}{a_{22}}

(4.219)

Aplicando esse raciocínio, podemos construir o método de Gauss-Seidel como $\begin{array}{rcl} x_{1}^{(k + 1)} & = & \frac{y_{1} - (a_{12} x_{2}^{(k)} + \dots + a_{1 n} x_{n}^{(k)})}{a_{11}} & (4.220) \\ x_{2}^{(k + 1)} & = & \frac{y_{2} - (a_{21} x_{1}^{(k + 1)} + a_{23} x_{3}^{(k)} + \dots + a_{2 n} x_{n}^{(k)})}{a_{22}} & (4.221) \\ ⋮ & (4.222) \\ x_{n}^{(k + 1)} & = & \frac{y_{2} - (a_{n 1} x_{1}^{(k + 1)} + \dots + a_{n (n - 1)} x_{n - 1}^{(k + 1)})}{a_{n n}} & (4.223) \end{array}$

Em notação mais compacta, o método de Gauss-Seidel consiste na iteração: $\begin{array}{rcl} x^{(1)} & = & aproximação inicial & (4.224) \\ x_{i}^{(k + 1)} & = & \frac{y_{i} - \sum_{j = 1}^{i - 1} a_{i j} x_{j}^{(k + 1)} - \sum_{j = i + 1}^{n} a_{i j} x_{j}^{(k)}}{a_{i i}} & (4.225) \end{array}$

Exemplo 4.7.3. Resolva o sistema $\begin{array}{rcl} 10 x + y & = & 23 & (4.226) \\ x + 8 y & = & 26 & (4.227) \end{array}$

usando o método de Gauss-Seidel com condições iniciais $x^{(1)} = y^{(1)} = 0$ . $\begin{array}{rcl} x^{(k + 1)} & = & \frac{23 - y^{(k)}}{10} & (4.228) \\ y^{(k + 1)} & = & \frac{26 - x^{(k + 1)}}{8} & (4.229) \\ x^{(2)} & = & \frac{23 - y^{(1)}}{10} = 2, 3 & (4.230) \\ y^{(2)} & = & \frac{26 - x^{(2)}}{8} = 2, 9625 & (4.231) \\ x^{(3)} & = & \frac{23 - y^{(2)}}{10} = 2, 00375 & (4.232) \\ y^{(3)} & = & \frac{26 - x^{(3)}}{8} = 2, 9995312 & (4.233) \end{array}$

Código Scilab: Método de Gauss-Seidel

function [x,deltax]=gauss_seidel(A,b,x,tol,N)
n=size(A,1)
xnew     =x
convergiu=%F                //FALSE;

k=1
while k<=N & ~convergiu
  xnew(1)=(b(1) - A(1,2:n)*x(2:n))/A(1,1)
  for i=2:n-1
    xnew(i)=(b(i) -A(i,1:i-1)*xnew(1:i-1) ...
                  -A(i,i+1:n)*x(i+1:n) )/A(i,i)
  end
  xnew(n)=(b(n) -A(n,1:n-1)*xnew(1:n-1) )/A(n,n)

  deltax=max( abs(x-xnew) )
  if deltax<tol then
     convergiu=%T       //TRUE
  end
  k=k+1
  x=xnew                 // atualiza x
  disp([k,x’,deltax])    // depuracao
end
if ~convergiu then
    error(’Nao convergiu’)
end

endfunction

4.7.3 Análise de convergência

Nesta seção, analisamos a convergência de métodos iterativos para solução de sistema lineares. Para tanto, consideramos um sistema linear $A x = b$ , onde $A = {[a_{i, j}]}_{i, j = 1}^{n, n}$ é a matriz (real) dos coeficientes, $b = {(a_{j})}_{j = 1}^{n}$ é um vetor dos termos constantes e $x = {(x_{j})}_{j = 1}^{n}$ é o vetor incógnita. No decorrer, assumimos que $A$ é uma matriz não singular.

Geralmente, métodos iterativos são construídos como uma iteração de ponto fixo. No caso de um sistema linear, reescreve-se a equação matricial em um problema de ponto fixo equivalente, isto é:

A x = b \Leftrightarrow x = T x + c,

(4.234)

onde $T = {[t_{i, j}]}_{i, j = 1}^{n, n}$ é chamada de matriz da iteração e $c = {(c_{j})}_{j = 1}^{n}$ de vetor da iteração. Construídos a matriz $T$ e o vetor $c$ , o método iterativo consiste em computar a iteração:

x^{(k + 1)} = T x^{(k)} + c, n \geq 1,

(4.235)

onde $x^{(1)}$ é uma aproximação inicial dada.

A fim de construirmos as matrizes e os vetores de iteração do método de Jacobi e de Gauss-Seidel, decompomos a matriz $A$ da seguinte forma:

A = L + D + U,

(4.236)

onde $D$ é a matriz diagonal $D = diag (a_{11}, a_{22}, \dots, a_{n n})$ , isto é:

D : = [\begin{matrix} a_{11} & 0 & 0 & \dots & 0 \\ 0 & a_{22} & 0 & \dots & 0 \\ 0 & 0 & a_{33} & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & a_{n n} \end{matrix}],

(4.237)

e, respectivamente, $L$ e $U$ são as seguintes matrizes triangular inferior e superior:

L : = [\begin{matrix} 0 & 0 & 0 & \dots & 0 \\ a_{21} & 0 & 0 & \dots & 0 \\ a_{31} & a_{32} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n 1} & a_{n 2} & a_{n 3} & \dots & 0 \end{matrix}], U : = [\begin{matrix} 0 & a_{12} & a_{13} & \dots & a_{1 n} \\ 0 & 0 & a_{23} & \dots & a_{2 n} \\ 0 & 0 & 0 & \dots & a_{3 n} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & 0 \end{matrix}] .

(4.238)

Exemplo 4.7.4. Considere o seguinte sistema linear: $\begin{array}{rcl} 3 x_{1} + x_{2} - x_{3} & = & 2 & (4.239) \\ - x_{1} - 4 x_{2} + x_{3} & = & - 10 & (4.240) \\ x_{1} - 2 x_{2} - 5 x_{3} & = & 10 & (4.241) \end{array}$

Escreva o sistema na sua forma matricial $A x = b$ identificando a matriz dos coeficientes $A$ , o vetor incógnita $x$ e o vetor dos termos constantes $b$ . Em seguida, faça a decomposição $A = L + D + U$ .

Solução. A forma matricial deste sistema é $A x = b$ , onde:

A = [\begin{matrix} 3 & 1 & - 1 \\ - 1 & - 4 & 1 \\ 1 & - 2 & - 5 \end{matrix}], x = [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}] e b = [\begin{matrix} 2 \\ - 10 \\ 10 \end{matrix}] .

(4.242)

A decomposição da matriz $A$ nas matrizes $L$ triangular inferior, $D$ diagonal e $U$ triangular superior é:

\underset{A}{\underset{︸}{[\begin{matrix} 3 & 1 & - 1 \\ - 1 & - 4 & 1 \\ 1 & - 2 & - 5 \end{matrix}]}} = \underset{L}{\underset{︸}{[\begin{matrix} 0 & 0 & 0 \\ - 1 & 0 & 0 \\ 1 & - 2 & 0 \end{matrix}]}} + \underset{D}{\underset{︸}{[\begin{matrix} 3 & 0 & 0 \\ 0 & - 4 & 0 \\ 0 & 0 & - 5 \end{matrix}]}} + \underset{U}{\underset{︸}{[\begin{matrix} 0 & 1 & - 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{matrix}]}} .

(4.243)

No Scilab, podemos construir as matrizes $L$ , $D$ e $U$ , da seguinte forma:

-->A = [3 1 -1;-1 -4 1;1 -2 -5];
-->D = eye(A).*A;
-->L = tril(A)-D;
-->U=triu(A)-D;

$♢$

Iteração de Jacobi

Vamos, agora, usar a decomposição discutida acima para construir a matriz de iteração $T_{J}$ e o vetor de iteração $c_{J}$ associado ao método de Jacobi. Neste caso, temos: $\begin{array}{rcl} A x = b & \Leftrightarrow & (L + D + U) x = b & (4.244) \\ \Leftrightarrow & D x = - (L + U) x + b & (4.245) \\ \Leftrightarrow & x = \underset{= : T_{J}}{\underset{︸}{- D^{- 1} (L + U)}} x + \underset{= : c_{J}}{\underset{︸}{D^{- 1} b}} . & (4.246) \end{array}$

Ou seja, a iteração do método de Jacobi escrita na forma matricial é:

x^{(k + 1)} = T_{J} x^{(k)} + c_{J}, k \geq 1,

(4.247)

com $x^{(1)}$ uma aproximação inicial dada, sendo $T_{J} : = - D^{- 1} (L + U)$ a matriz de iteração e $c_{J} = D^{- 1} b$ o vetor da iteração.

Exemplo 4.7.5. Construa a matriz de iteração $T_{J}$ e o vetor de iteração $c_{J}$ do método de Jacobi para o sistema dado no Exemplo 4.7.4.

Solução. A matriz de iteração é dada por:

T_{J} : = - D^{- 1} (L + U) = - \underset{D^{- 1}}{\underset{︸}{[\begin{matrix} \frac{1}{3} & 0 & 0 \\ 0 & - \frac{1}{4} & 0 \\ 0 & 0 & - \frac{1}{5} \end{matrix}]}} \underset{(L + U)}{\underset{︸}{[\begin{matrix} 0 & 1 & - 1 \\ - 1 & 0 & 1 \\ 1 & 2 & 0 \end{matrix}]}} = [\begin{matrix} 0 & - \frac{1}{3} & \frac{1}{3} \\ - \frac{1}{4} & 0 & \frac{1}{4} \\ \frac{1}{5} & \frac{2}{5} & 0 \end{matrix}] .

(4.248)

O vetor da iteração de Jacobi é:

c_{J} : = D^{- 1} b = \underset{D^{- 1}}{\underset{︸}{[\begin{matrix} \frac{1}{3} & 0 & 0 \\ 0 & - \frac{1}{4} & 0 \\ 0 & 0 & - \frac{1}{5} \end{matrix}]}} \underset{b}{\underset{︸}{[\begin{matrix} 2 \\ - 10 \\ 10 \end{matrix}]}} = [\begin{matrix} \frac{2}{3} \\ \frac{5}{2} \\ - 2 \end{matrix}] .

(4.249)

No Scilab, podemos computar $T_{J}$ e $c_{J}$ da seguinte forma:

-->TJ = -inv(D)*(L+U);
-->cJ = inv(D)*b;

$♢$

Iteração de Gauss-Seidel

A forma matricial da iteração do método de Gauss-Seidel também pode ser construída com base na decomposição $A = L + D + U$ . Para tando, fazemos: $\begin{array}{rcl} A x = b & \Leftrightarrow & (L + D + U) x = b & (4.250) \\ \Leftrightarrow & (L + D) x = - U x + b & (4.251) \\ \Leftrightarrow & x = \underset{= : T_{G}}{\underset{︸}{- {(L + D)}^{- 1} U}} x + \underset{= : c_{G}}{\underset{︸}{{(L + D)}^{- 1} b}} & (4.252) \end{array}$

Ou seja, a iteração do método de Gauss-Seidel escrita na forma matricial é:

x^{(k + 1)} = T_{G} x^{(k)} + c_{G}, k \geq 1,

(4.253)

com $x^{(1)}$ uma aproximação inicial dada, sendo $T_{G} : = - {(L + D)}^{- 1} U$ a matriz de iteração e $c_{J} = {(L + D)}^{- 1} b$ o vetor da iteração.

Exemplo 4.7.6. Construa a matriz de iteração $T_{G}$ e o vetor de iteração $c_{G}$ do método de Gauss-Seidel para o sistema dado no Exemplo 4.7.4.

Solução. A matriz de iteração é dada por:

T_{G} = - {(L + D)}^{- 1} U = - \underset{{(L + D)}^{- 1}}{\underset{︸}{{[\begin{matrix} 3 & 0 & 0 \\ - 1 & - 4 & 0 \\ 1 & - 2 & - 5 \end{matrix}]}^{- 1}}} \underset{U}{\underset{︸}{[\begin{matrix} 0 & 1 & - 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{matrix}]}} = [\begin{matrix} 0 & - \frac{1}{3} & \frac{1}{3} \\ 0 & \frac{1}{12} & \frac{1}{6} \\ 0 & - \frac{1}{10} & 0 \end{matrix}] .

(4.254)

O vetor da iteração de Gauss-Seidel é:

c_{G} : = {(L + D)}^{- 1} b = \underset{{(L + D)}^{- 1}}{\underset{︸}{{[\begin{matrix} 3 & 0 & 0 \\ - 1 & - 4 & 0 \\ 1 & - 2 & - 5 \end{matrix}]}^{- 1}}} \underset{b}{\underset{︸}{[\begin{matrix} 2 \\ - 10 \\ 10 \end{matrix}]}} = [\begin{matrix} \frac{2}{3} \\ \frac{7}{3} \\ - \frac{28}{10} \end{matrix}] .

(4.255)

No Scilab, podemos computar $T_{G}$ e $c_{G}$ da seguinte forma:

-->TG = -inv(L+D)*U;
-->cG = inv(L+D)*b;

$♢$

Condições de convergência

Aqui, vamos discutir condições necessárias e suficientes para a convergência de métodos iterativos. Isto é, dado um sistema $A x = b$ e uma iteração:

x^{(k + 1)} = T x^{(k)} + c, k \geq 1,

(4.256)

$x^{(1)}$ dado, estabelecemos condições nas quais $x^{(k)} \to x^{*}$ , onde $x^{*}$ é a solução do sistema dado, isto é, $x^{*} = T x^{*} + c$ ou, equivalentemente, $A x^{*} = b$ .

Lema 4.7.1. Seja $T$ uma matriz real $n \times n$ . O limite $lim_{k \to \infty} {∥T^{k}∥}_{p} = 0$ , $1 \leq p \leq \infty$ , se, e somente se, $ρ (T) < 1$ .

Demonstração. Aqui, fazemos apenas um esboço da demonstração. Para mais detalhes, veja [8], Teorema 4, pág. 14.

Primeiramente, suponhamos que ${∥T∥}_{p} < 1$ , $1 \leq p \leq \infty$ . Como (veja [8], lema 2, pág. 12):

ρ (T) \leq {∥T∥}_{p},

(4.257)

temos $ρ (T) < 1$ , o que mostra a implicação.

Agora, suponhamos que $ρ (T) < 1$ e seja $0 < ϵ < 1 - ρ (T)$ . Então, existe $1 \leq p \leq \infty$ tal que (veja [8], Teorema 3, página 12):

{∥T∥}_{p} \leq ρ (T) + ϵ < 1 .

(4.258)

Assim, temos:

lim_{k \to \infty} ∥ T^{k} ∥_{p} \leq lim_{k \to \infty} {∥T∥}_{p}^{m} = 0 .

(4.259)

Da equivalência entre as normas segue a recíproca.

$■$

Observação 4.7.1. Observamos que:

lim_{k \to \infty} ∥ T^{k} ∥_{p} = 0,, 1 \leq p \leq \infty, \Leftrightarrow lim_{k \to \infty} t_{i j}^{k} = 0, 1 \leq i, j \leq n .

(4.260)

Lema 4.7.2. Se $ρ (T) < 1$ , então existe ${(I - T)}^{- 1}$ e:

{(I - T)}^{- 1} = \sum_{k = 0}^{\infty} T^{k} .

(4.261)

Demonstração. Primeiramente, provamos a existência de

{(I - T)}^{- 1}

. Seja

λ

um autovalor de

T

x

um autovetor associado, isto é,

T x = λ x

. Então,

(I - T) x = (1 - λ) x

. Além disso, temos

| λ | < ρ (T) < 1

, logo

(1 - λ) \neq 0

, o que garante que

(I - T)

é não singular. Agora, mostramos que

{(I - T)}^{- 1}

admite a expansão acima. Do Lema 4.7.1 e da Observação 4.7.1 temos:

(I - T) \sum_{k = 0}^{\infty} T^{k} = lim_{m \to \infty} (I - T) \sum_{k = 0}^{m} T^{k} = lim_{m \to \infty} (I - T^{m + 1}) = I,

(4.262)

o que mostra que ${(I - T)}^{- 1} = \sum_{k = 0}^{\infty} T^{k}$ .

$■$

Teorema 4.7.1. A sequência recursiva ${x^{(k)}}_{k \in ℕ}$ dada por:

x^{(k + 1)} = T x^{(k)} + c

(4.263)

converge para solução de $x = T x + c$ para qualquer escolha de $x^{(1)}$ se, e somente se, $ρ (T) < 1$ .

Demonstração. Primeiramente, assumimos que

ρ (T) < 1

. Observamos que:

\begin{array}{rcl} x^{(k + 1)} & = T x^{(k)} + c = T (T x^{(k - 1)} + c) + c & (4.264) \\ = T^{2} x^{(k - 1)} + (I + T) c & (4.265) \\ ⋮ & (4.266) \\ = T^{(k)} x^{(1)} + (\sum_{k = 0}^{k - 1} T^{k}) c . & (4.267) \end{array}

Daí, do Lema 4.7.1 e do Lema 4.7.2 temos:

lim_{k \to \infty} x^{(k)} = {(I - T)}^{(- 1)} c .

(4.268)

Ora, se $x^{*}$ é a solução de $x = T x + c$ , então $(I - T) x^{*} = c$ , isto é, $x^{*} = {(I - T)}^{- 1} c$ . Logo, temos demonstrado que $x^{(k)}$ converge para a solução de $x = T x + c$ , para qualquer escolha de $x^{(1)}$ .

Agora, suponhamos que $x^{(k)}$ converge para $x^{*}$ solução de $x = T x + c$ , para qualquer escolha de $x^{(1)}$ . Seja, então, $y$ um vetor arbitrário e $x^{(1)} = x^{*} - y$ . Observamos que: $\begin{array}{rcl} x^{*} - x^{(k + 1)} & = (T x^{*} + c) - (T x^{(k)} + c) & (4.269) \\ = T (x^{*} - x^{(k)}) & (4.270) \\ ⋮ & (4.271) \\ = T^{(k)} (x^{*} - x^{(1)}) = T^{(k)} y . & (4.272) \end{array}$

Logo, para qualquer $1 \leq p \leq \infty$ , temos, :

0 = lim_{k \to \infty} x^{*} - x^{(k + 1)} = lim_{k \to \infty} T^{(k)} y .

(4.273)

Como $y$ é arbitrário, da Observação 4.7.1 temos $lim_{k \to \infty} ∥ T^{(k)} ∥_{p} = 0$ , $1 \leq p \leq \infty$ . Então, o Lema 4.7.1 garante que $ρ (T) < 1$ .

$■$

Observação 4.7.2. Pode-se mostrar que tais métodos iterativos tem taxa de convergência super linear com:

∥ x^{(k + 1)} - x^{*} ∥ \approx ρ {(T)}^{k} ∥ x^{(1)} - x^{*} ∥ .

(4.274)

Para mais detalhes, veja [8], pág. 61-64.

Exemplo 4.7.7. Mostre que, para qualquer escolha da aproximação inicial, ambos os métodos de Jacobi e Gauss-Seidel são convergentes quando aplicados ao sistema linear dado no Exemplo 4.7.4.

Solução. Do Teorema 4.7.1, vemos que é necessário e suficiente que $ρ (T_{J}) < 1$ e $ρ (T_{G}) < 1$ . Computando estes raios espectrais, obtemos $ρ (T_{J}) \approx 0, 32$ e $ρ (T_{G}) \approx 0, 13$ . Isto mostra que ambos os métodos serão convergentes.

$♢$

Condição suficiente

Uma condição suficiente porém não necessária para que os métodos de Gauss-Seidel e Jacobi convirjam é a que a matriz seja estritamente diagonal dominante.

Definição 4.7.1. Uma matriz $A$ é estritamente diagonal dominante quando:

|a_{i i}| > \sum_{\begin{array}{c} j = 1 \\ j \neq i \end{array}}^{n} |a_{i j}|, i = 1, . . ., n

(4.275)

Definição 4.7.2. Uma matriz $A$ é diagonal dominante quando

|a_{i i}| \geq \sum_{\begin{array}{c} j = 1 \\ j \neq i \end{array}}^{n} |a_{i j}|, i = 1, . . ., n

(4.276)

e para ao menos um $i$ , $a_{i i}$ é estritamente maior que a soma dos elementos fora da diagonal.

Teorema 4.7.2. Se a matriz $A$ for diagonal dominante⁹ , então os métodos de Jacobi e Gauss-Seidel serão convergentes independente da escolha inicial $x^{(1)}$ .

Se conhecermos a solução exata $x$ do problema, podemos calcular o erro relativo em cada iteração como:

\frac{∥ x - x^{(k)} ∥}{∥x∥} .

(4.277)

Em geral não temos $x$ , entretanto podemos estimar o vetor resíduo $r^{(k)} = b - A \tilde{x^{(k)}}$ . Note que quando o erro tende a zero, o resíduo também tende a zero.

Teorema 4.7.3. O erro relativo e o resíduo estão relacionados como (veja [3])

\frac{∥ x - x^{(k)} ∥}{∥x∥} \leq κ (A) \frac{∥ r ∥}{∥ b ∥}

(4.278)

onde $k (A)$ é o número de condicionamento.

Exemplo 4.7.8. Ambos os métodos de Jacobi e Gauss-Seidel são convergentes para o sistema dado no Exemplo 4.7.4, pois a matriz dos coeficientes deste é uma matriz estritamente diagonal dominante.

Exercícios

E 4.7.1. Considere o problema de 5 incógnitas e cinco equações dado por $\begin{array}{rcl} x_{1} - x_{2} & = & 1 & (4.279) \\ - x_{1} + 2 x_{2} - x_{3} & = & 1 & (4.280) \\ - x_{2} + (2 + ε) x_{3} - x_{4} & = & 1 & (4.281) \\ - x_{3} + 2 x_{4} - x_{5} & = & 1 & (4.282) \\ x_{4} - x_{5} & = & 1 & (4.283) \end{array}$

Escreva na forma $A x = b$ e resolva usando eliminação gaussiana para $ε = 1 0^{- 3}$ no Scilab.
Obtenha o vetor incógnita $x$ com $ε = 1 0^{- 3}$ usando Jacobi com tolerância $1 0^{- 2}$ . Compare o resultado com o resultado obtido no item d.
Obtenha o vetor incógnita $x$ com $ε = 1 0^{- 3}$ usando Gauss-Seidel com tolerância $1 0^{- 2}$ . Compare o resultado com o resultado obtido no item d.
Discuta com base na relação esperada entre tolerância e exatidão conforme estudado na primeira área para problemas de uma variável.

Resposta.

epsilon=1e-3;

A=[1 -1 0 0 0; -1 2 -1 0 0; 0 -1 (2+epsilon) -1 0; 0 0 -1 2 -1; 0 0 0 1 -1]

v=[1 1 1 1 1]’
xgauss=gauss([A v])

function x=q_Jacobi()
    x0=[0 0 0 0 0]’

    i=0
    controle=0
    while controle<3 & i<1000
    i=i+1

    x(1)=1+x0(2)
    x(2)=(1+x0(3)+x0(1))/2
    x(3)=(1+x0(2)+x0(4))/(2+epsilon)
    x(4)=(1+x0(3)+x0(5))/2
    x(5)=x0(4)-1

    delta=norm(x-x0,2)
    if delta<1e-6 then
        controle=controle+1
    else
        controle=0
    end
    mprintf(’i=%d, x1=%f, x5=%f, tol=%.12f\n’, i, x(1), x(5), delta)
    x0=x;
    end

endfunction

function x=q_Gauss_Seidel()
    x0=[0 0 0 0 0]’

    i=0
    controle=0
    while controle<3 & i<15000
    i=i+1

    x(1)=1+x0(2)
    x(2)=(1+x0(3)+x(1))/2
    x(3)=(1+x(2)+x0(4))/(2+epsilon)
    x(4)=(1+x(3)+x0(5))/2
    x(5)=x(4)-1

    delta=norm(x-x0,2)
    if delta<1e-2 then
        controle=controle+1
    else
        controle=0
    end
    mprintf(’i=%d, x1=%f, x5=%f, tol=%.12f\n’,i,x(1),x(5),delta)
    x0=x;
    end

endfunction

E 4.7.2. Resolva o seguinte sistema pelo método de Jacobi e Gauss-Seidel:

\{\begin{matrix} 5 x_{1} + x_{2} + x_{3} & = 50 \\ - x_{1} + 3 x_{2} - x_{3} & = 10 \\ x_{1} + 2 x_{2} + 10 x_{3} & = - 30 \end{matrix}

(4.284)

Use como critério de paragem tolerância inferior a $1 0^{- 3}$ e inicialize com $x^{0} = y^{0} = z^{0} = 0$ .

E 4.7.3. Refaça o Exercício 4.1.6 construindo um algoritmo que implemente os métodos de Jacobi e Gauss-Seidel.

E 4.7.4. Considere o seguinte sistema de equações lineares: $\begin{array}{rcl} x_{1} - x_{2} & = & 0 \\ - x_{j - 1} + 5 x_{j} - x_{j + 1} & = & cos (j ∕ 10), 2 \leq j \leq 10 \\ x_{11} & = & x_{10} ∕ 2 & (4.285) \end{array}$

Construa a iteração para encontrar a solução deste problema pelos métodos de Gauss-Seidel e Jacobi. Usando esses métodos, encontre uma solução aproximada com erro absoluto inferior a $1 0^{- 5}$ . Veja também Exercício 4.5.2

Resposta.

$0, 324295$ , $0, 324295$ , $0, 317115$ , $0, 305943$ , $0, 291539$ , $0, 274169$ , $0, 253971$ , $0, 230846$ , $0, 203551$ , $0, 165301$ , $0, 082650$

Exemplos de rotinas:

function x=jacobi()
    x0=zeros(11,1)
    k=0;
    controle=0;
    while controle<3 & k<1000
        k=k+1
        x(1)=x0(2)
        for j=2:10
        x(j)=(cos(j/10)+x0(j-1)+x0(j+1))/5
        end
        x(11)=x0(10)/2

        delta=norm(x-x0) //norma 2
        if delta<1e-5 then
            controle=controle+1
        else
            controle=0;
        end
        mprintf(’k=%d, x=[%f,%f,%f], tol=%.12f\n’,k,x(1),x(2),x(3),delta)
        x0=x;
    end

endfunction

function x=gs()
    x0=zeros(11,1)
    k=0;
    controle=0;
    while controle<3 & k<1000
        k=k+1
        x(1)=x0(2)
        for j=2:10
        x(j)=(cos(j/10)+x(j-1)+x0(j+1))/5
        end
        x(11)=x0(10)/2

        delta=norm(x-x0) //norma 2
        if delta<1e-5 then
            controle=controle+1
        else
            controle=0;
        end
        mprintf(’k=%d, x=[%f,%f,%f], tol=%.12f\n’,k,x(1),x(2),x(3),delta)
        x0=x;
    end
endfunction

E 4.7.5. Faça uma permutação de linhas no sistema abaixo e resolva pelos métodos de Jacobi e Gauss-Seidel: $\begin{array}{rcl} x_{1} + 10 x_{2} + 3 x_{3} = 27 & (4.286) \\ 4 x_{1} + x_{3} = 6 & (4.287) \\ 2 x_{1} + x_{2} + 4 x_{3} = 12 & (4.288) \end{array}$

Resposta.

Permute as linhas 1 e 2.

Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 19/8/2020 às 17:36:32.

Recursos

Projeto

IME - UFRGS
- Página do IME
- Página da UFRGS

UFRGS - IME - Recursos Educacionais Abertos de Matemática. Contato: reamat@ufrgs.br.

Informe erros ou edite você mesmo!