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.