Como comparar aproximação de matrizes

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.

O método de Jacobi pode ser obtido a partir do sistema linear a11x1 + a12x2 + ⋯ + a1nxn = y1 (4.181) a21x1 + a22x2 + ⋯ + a2nxn = y2 (4.182) ⋮ (4.183) an1x1 + an2x2 + ⋯ + annxn = yn (4.184)

Isolando o elemento x1 da primeira equação temos

x1(k+1) = y1 − a12x2(k) + ⋯ + a 1nxn(k) a11

(4.185)

Note que utilizaremos os elementos xi(k) da iteração k (a direita da equação) para estimar o elemento x1 da próxima iteração.

Da mesma forma, isolando o elemento xi de cada equação i, para todo i = 2,...,n podemos construir a iteração x1(k+1) = y1 − a12x2(k) + ⋯ + a 1nxn(k) a11 (4.186) x2(k+1) = y2 − a21x1(k) + a 23x3(k) + ⋯ + a 2nxn(k) a22 (4.187) ⋮ (4.188) xn(k+1) = y2 − an1x1(k) + ⋯ + a n,n−2xn−2(k) + a n,n−1xn−1(k) ann (4.189)

Em notação mais compacta, o método de Jacobi consiste na iteração x(1) = aproximação inicial (4.190) xi(k+1) = y i −∑ j=1 j≠i na ijxj(k) ∕a ii (4.191)

Exemplo 4.7.1. Resolva o sistema 10x + y = 23 (4.192) x + 8y = 26 (4.193)

usando o método de Jacobi iniciando com x(1) = y(1) = 0. x(k+1) = 23 − y(k) 10 (4.194) y(k+1) = 26 − x(k) 8 (4.195) x(2) = 23 − y(1) 10 = 2,3 (4.196) y(2) = 26 − x(1) 8 = 3,25 (4.197) x(3) = 23 − y(2) 10 = 1,975 (4.198) y(3) = 26 − x(2) 8 = 2,9625 (4.199)

Exemplo 4.7.2. Considere o seguinte sistema

− 3x1 + x2 + x3 = 2 2x1 + 5x2 + x3 = 5 2x1 + 3x2 + 7x3 = −17

(4.200)

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)∥ ∞

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,5E −1
5	(−0,91,1,92, − 3,07)	2,4E −1
⋮	⋮	⋮
10	(−1,00,2,00, − 3,00)	6,0E −3

Verifique a resposta.

No GNU Octave, podemos fazer as contas com o seguinte script:

A = [-3 1 1;2 5 1;2 3 7]; b = [2,5,-17]’; x=[1 1 -1]’; for iter=1:9 x0 = x; for i=1:3 x(i) = (b(i)-A(i,[1:i-1,i+1:3])*x0([1:i-1,i+1:3]))/A(i,i); endfor x’ norm(x-x0,’inf’)

endfor

function [x]=jacobi(A,b,x,tol,N)
  n = size(A,1);
  x0 = x;
  iter = 1;
  while (iter <= N)
    #iteracao Jacobi
    for i = 1:n
      x(i) = (b(i) - A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);
    endfor
    #condicao de parada
    if (norm(x-x0,’inf’)<tol)
      return
    endif
    #prepara nova iteracao
    x0 = x;
    iter += 1;
  endwhile
  error(’num. max. iter. excedido.’)
endfunction

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

x2(k+1) = y2 − a21x1(k+1) + a 23x3(k) + ⋯ + a 2nxn(k) a22

(4.201)

Aplicando esse raciocínio, podemos construir o método de Gauss-Seidel como x1(k+1) = y1 − a12x2(k) + ⋯ + a 1nxn(k) a11 (4.202) x2(k+1) = y2 − a21x1(k+1) + a 23x3(k) + ⋯ + a 2nxn(k) a22 (4.203) ⋮ (4.204) xn(k+1) = y2 − an1x1(k+1) + ⋯ + a n(n−1)xn−1(k+1) ann (4.205)

Em notação mais compacta, o método de Gauss-Seidel consiste na iteração: x(1) = aproximação inicial (4.206) xi(k+1) = yi −∑ j=1i−1a ijxj(k+1) −∑ j=i+1na ijxj(k) aii (4.207)

Exemplo 4.7.3. Resolva o sistema 10x + y = 23 (4.208) x + 8y = 26 (4.209)

usando o método de Gauss-Seidel com condições iniciais x(1) = y(1) = 0. x(k+1) = 23 − y(k) 10 (4.210) y(k+1) = 26 − x(k+1) 8 (4.211) x(2) = 23 − y(1) 10 = 2,3 (4.212) y(2) = 26 − x(2) 8 = 2,9625 (4.213) x(3) = 23 − y(2) 10 = 2,00375 (4.214) y(3) = 26 − x(3) 8 = 2,9995312 (4.215)

function [x]=gauss_seidel(A,b,x,tol,N)
  n = size(A,1);
  x0 = x;
  iter = 1;
  while (iter <= N)
    #iteracao G-S
    for i = 1:n
      x(i) = (b(i) - A(i,[1:i-1,i+1:n])*x([1:i-1,i+1:n]))/A(i,i);
    endfor
    #condicao de parada
    if (norm(x-x0,’inf’)<tol)
      return
    endif
    #prepara nova iteracao
    x0 = x;
    iter += 1;
  endwhile
  error(’num. max. iter. excedido.’)
endfunction

Nesta seção, analisamos a convergência de métodos iterativos para solução de sistema lineares. Para tanto, consideramos um sistema linear Ax = b, onde A = [ai,j]i,j=1n,n é a matriz (real) dos coeficientes, b = (aj)j=1n é um vetor dos termos constantes e x = (xj)j=1n é 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 é:

Ax = b ⇔ x = Tx + c,

(4.216)

onde T = [ti,j]i,j=1n,n é chamada de matriz da iteração e c = (cj)j=1n 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) = Tx(k) + c,n ≥ 1,

(4.217)

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:

onde D é a matriz diagonal D = diag (a11,a22,…,ann), isto é:

D := a11 0 0 ⋯ 0 0 a22 0 ⋯ 0 0 0 a33 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ ann ,

(4.219)

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

L := 0 0 0 ⋯ 0 a 21 0 0 ⋯ 0 a 31 a32 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ an1 an2 an3 ⋯ 0 ,U := 0 a12 a13 ⋯ a1n 0 0 a23 ⋯ a2n 0 0 0 ⋯ a3n ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 0 .

(4.220)

Exemplo 4.7.4. Considere o seguinte sistema linear: 3x1 + x2 − x3 = 2 (4.221) −x1 − 4x2 + x3 = −10 (4.222) x1 − 2x2 − 5x3 = 10 (4.223)

Escreva o sistema na sua forma matricial Ax = 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 é Ax = b, onde:

A = 3 1 −1 −1 −4 1 1 −2 −5 ,x = x1 x2 x3 eb = 2 −10 10 .

(4.224)

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

3 1 −1 −1 −4 1 1 −2 −5 ︸ A = 0 0 0 −1 0 0 1 −2 0 ︸ L+ 3 0 0 0 −4 0 0 0 −5 ︸ D+ 0 1 −1 0 0 1 0 0 0 ︸ U.

(4.225)

Vamos, agora, usar a decomposição discutida acima para construir a matriz de iteração TJ e o vetor de iteração cJ associado ao método de Jacobi. Neste caso, temos: Ax = b ⇔ (L + D + U)x = b (4.226) ⇔ Dx = −(L + U)x + b (4.227) ⇔ x = −D−1(L + U)︸ =:TJx + D−1b︸ =:cJ. (4.228)

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

x(k+1) = T Jx(k) + c J,k ≥ 1,

(4.229)

com x(1) uma aproximação inicial dada, sendo TJ := −D−1(L + U) a matriz de iteração e cJ = D−1b o vetor da iteração.

Exemplo 4.7.5. Construa a matriz de iteração TJ e o vetor de iteração cJ do método de Jacobi para o sistema dado no Exemplo 4.7.4.

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

TJ := −D−1(L+U) = − 1 3 0 0 0 −1 4 0 0 0 −1 5 ︸ D−1 0 1 −1 −1 0 1 1 2 0 ︸ (L+U) = 0 −1 3 1 3 −1 4 0 1 4 1 5 2 5 0 .

(4.230)

O vetor da iteração de Jacobi é:

cJ := D−1b = 1 3 0 0 0 −1 4 0 0 0 −1 5 ︸ D−1 2 −10 10 ︸ b = 2 3 5 2 −2 .

(4.231)

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: Ax = b ⇔ (L + D + U)x = b (4.232) ⇔ (L + D)x = −Ux + b (4.233) ⇔ x = −(L + D)−1U︸ =:TGx + (L + D)−1b︸ =:cG (4.234)

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

x(k+1) = T Gx(k) + c G,k ≥ 1,

(4.235)

com x(1) uma aproximação inicial dada, sendo TG := −(L + D)−1U a matriz de iteração e cJ = (L + D)−1b o vetor da iteração.

Exemplo 4.7.6. Construa a matriz de iteração TG e o vetor de iteração cG do método de Gauss-Seidel para o sistema dado no Exemplo 4.7.4.

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

TG = −(L+D)−1U = − 3 0 0 −1 −4 0 1 −2 −5 −1︸ (L+D)−1 0 1 −1 0 0 1 0 0 0 ︸ U = 0 −1 3 1 3 0 1 12 1 6 0 −1 10 0 .

(4.236)

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

cG := (L+D)−1b = 3 0 0 −1 −4 0 1 −2 −5 −1︸ (L+D)−1 2 −10 10 ︸ b = 2 3 7 3 −28 10 .

(4.237)

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

x(k+1) = Tx(k) + c,k ≥ 1,

(4.238)

x(1) dado, estabelecemos condições nas quais x(k) → x∗, onde x∗ é a solução do sistema dado, isto é, x∗ = Tx∗ + c ou, equivalentemente, Ax∗ = b.

Lema 4.7.1. Seja T uma matriz real n × n. O limite lim k→∞ Tk p = 0, 1 ≤ p ≤∞, 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 Tp < 1, 1 ≤ p ≤∞. Como (veja [8], lema 2, pág. 12):

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

Agora, suponhamos que ρ(T) < 1 e seja 0 < ϵ < 1 − ρ(T). Então, existe 1 ≤ p ≤∞ tal que (veja [8], Teorema 3, página 12):

Tp ≤ ρ(T) + ϵ < 1.

(4.240)

Assim, temos:

lim k→∞∥Tk∥ p ≤ lim k→∞ Tpm = 0.

(4.241)

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

Observação 4.7.1. Observamos que:

lim k→∞∥Tk∥ p = 0,,1 ≤ p ≤∞, ⇔ lim k→∞tijk = 0,1 ≤ i,j ≤ n.

(4.242)

Lema 4.7.2. Se ρ(T) < 1, então existe (I − T)−1 e:

(I − T)−1 = ∑ k=0∞Tk.

(4.243)

Demonstração. Primeiramente, provamos a existência de (I − T)−1. Seja λ um autovalor de T e x um autovetor associado, isto é, Tx = λx. Então, (I − T)x = (1 − λ)x. Além disso, temos |λ| < ρ(T) < 1, logo (1 − λ)≠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) ∑ k=0∞Tk = lim m→∞(I −T) ∑ k=0mTk = lim m→∞(I −Tm+1) = I,

(4.244)

o que mostra que (I − T)−1 = ∑ k=0∞Tk.

Teorema 4.7.1. A sequência recursiva {x(k)} k∈ℕ dada por:

x(k+1) = Tx(k) + c

(4.245)

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

Demonstração. Primeiramente, assumimos que ρ(T) < 1. Observamos que: x(k+1) = Tx(k) + c = T(Tx(k−1) + c) + c (4.246) = T2x(k−1) + (I + T)c (4.247) ⋮ (4.248) = T(k)x(1) + ∑ k=0k−1Tk c. (4.249)

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

lim k→∞x(k) = (I − T)(−1)c.

(4.250)

Ora, se x∗ é a solução de x = Tx + c, então (I − T)x∗ = c, isto é, x∗ = (I − T)−1c. Logo, temos demonstrado que x(k) converge para a solução de x = Tx + c, para qualquer escolha de x(1).

Agora, suponhamos que x(k) converge para x∗ solução de x = Tx + c, para qualquer escolha de x(1). Seja, então, y um vetor arbitrário e x(1) = x∗ − y. Observamos que: x∗ − x(k+1) = (Tx∗ + c) − (Tx(k) + c) (4.251) = T(x∗ − x(k)) (4.252) ⋮ (4.253) = T(k)(x∗ − x(1)) = T(k)y. (4.254)

Logo, para qualquer 1 ≤ p ≤∞, temos, :

0 = lim k→∞x∗ − x(k+1) = lim k→∞T(k)y.

(4.255)

Como y é arbitrário, da Observação 4.7.1 temos lim k→∞∥T(k)∥ p = 0, 1 ≤ p ≤∞. 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∗∥≈ ρ(T)k∥x(1) − x∗∥.

(4.256)

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 ρ(TJ) < 1 e ρ(TG) < 1. Computando estes raios espectrais, obtemos ρ(TJ) ≈ 0,32 e ρ(TG) ≈ 0,13. Isto mostra que ambos os métodos serão convergentes.

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:

aii > ∑ j=1 j≠i n a ij ,i = 1,...,n

(4.257)

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

aii ≥∑ j=1 j≠i n a ij ,i = 1,...,n

(4.258)

e para ao menos um i, aii é 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:

Em geral não temos x, entretanto podemos estimar o vetor resíduo r(k) = b − Ax(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])

∥x − x(k)∥ x ≤ κ(A)∥r∥ ∥b∥

(4.260)

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.

E 4.7.1. Considere o problema de 5 incógnitas e cinco equações dado por x1 − x2 = 1 (4.261) −x1 + 2x2 − x3 = 1 (4.262) −x2 + (2 + ε)x3 − x4 = 1 (4.263) −x3 + 2x4 − x5 = 1 (4.264) x4 − x5 = 1 (4.265)

Escreva na forma Ax = b e resolva usando eliminação gaussiana para ε = 10−3 no Octave.
Obtenha o vetor incógnita x com ε = 10−3 usando Jacobi com tolerância 10−2. Compare o resultado com o resultado obtido no item d.
Obtenha o vetor incógnita x com ε = 10−3 usando Gauss-Seidel com tolerância 10−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.

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

5x1 + x2 + x3 = 50 − x1 + 3x2 − x3 = 10 x1 + 2x2 + 10x3 = −30

(4.266)

Use como critério de paragem tolerância inferior a 10−3 e inicialize com x0 = y0 = z0 = 0.

E 4.7.3. Refaça o Exercício ?? 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: x1 − x2 = 0 −xj−1 + 5xj − xj+1 = cos(j∕10),2 ≤ j ≤ 10 x11 = x10∕2 (4.267)

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 10−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

E 4.7.5. Faça uma permutação de linhas no sistema abaixo e resolva pelos métodos de Jacobi e Gauss-Seidel: x1 + 10x2 + 3x3 = 27 (4.268) 4x1 + x3 = 6 (4.269) 2x1 + x2 + 4x3 = 12 (4.270)

Resposta.

Permute as linhas 1 e 2.