| | | |

Matemática numérica I

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

3.3 Métodos de Jacobi e de Gauss-Seidel

Nesta seção, estudamos os métodos iterativos de Jacobi777Carl Gustav Jakob Jacobi, 1804 - 1851, matemático alemão. Fonte: Wikipédia: Carl Gustav Jakob Jacobi. e de Gauss888Johann Carl Friedrich Gauss, 1777 - 1855, matemático alemão. Fonte: Wikipédia: Carl Friedrich Gauss.-Seidel999Philipp Ludwig von Seidel, 1821 - 1896, matemático alemão. Fonte: Wikipédia: Philipp Ludwig von Seidel. para a aproximação da solução de sistemas lineares

A⁢𝒙=𝒃, (3.232)

onde A=[ai,j]i,j=1n,n, n≥1, é uma dada matriz dos coeficientes, 𝒙=(x1,x2,…,xn) é o vetor das incógnitas e 𝒃=(b1,b2,…,bn) é um dado vetor dos termos constantes.

Métodos iterativos para sistemas lineares têm a forma

𝒙(0)=aprox. inicial, (3.233)
𝒙(k+1)=T⁢𝒙(k)+𝒄, (3.234)

onde T=[ti,j]i,j=1n,n é a matriz de iteração e 𝒄=(c1,c2,…,cn) é o vetor de iteração.

3.3.1 Método de Jacobi

Consideramos a seguinte decomposição da matriz A=L+D+U:

A=[a11a12a13…a1⁢na21a22a23…a2⁢na31a32a33…a3⁢n⋮⋮⋮…⋮an⁢1an⁢2an⁢3…an⁢n]
=[000…0a2100…0a31a320…0⋮⋮⋮…⋮an⁢1an⁢2an⁢3…0]⏟L
+[a1100…00a220…000a33…0⋮⋮⋮…⋮000…an⁢n]⏟D
+[0a12a13…a1⁢n00a23…a2⁢n00a33…a3⁢n⋮⋮⋮…⋮000…an⁢n]⏟U. (3.235)

Isto é, a matriz A decomposta como a soma de sua parte triangular inferior L, de sua diagonal D e de sua parte triangular superior U.

Desta forma, podemos reescrever o sistema como segue:

A⁢𝒙=𝒃 (3.236)
(L+D+U)⁢𝒙=𝒃 (3.237)
(L+U)⁢𝒙+D⁢𝒙=𝒃 (3.238)
D⁢𝒙=−(L+U)⁢𝒙+𝒃 (3.239)
𝒙=−D−1⁢(L+U)⁢𝒙+D−1⁢𝒃. (3.240)

Ou seja, resolver o sistema A⁢𝒙=𝒃 é equivalente a resolver o problema de ponto fixo

𝒙=TJ⁢𝒙+𝒄J, (3.241)

onde TJ=−D−1⁢(L+U) é chamada de matriz de Jacobi e 𝒄J=D−1⁢𝒃 é chamado de vetor de Jacobi.

Exemplo 3.3.1.

Consideramos o sistema linear A⁢𝒙=𝒃 com

A=[−42−1−2521−1−3], (3.242)
𝒃=[−11−70]. (3.243)

Este sistema tem solução 𝒙=(2,−1,1). Neste caso, temos a decomposição A=L+D+U com

L=[000−2001−10], (3.244)
D=[−40005000−3], (3.245)
U=[02−1002000]. (3.246)

Ainda, observamos que

TJ⁢𝒙+𝒄J=−D−1⁢(L+U)⁢𝒙+D−1⁢𝒃 (3.247)
=[01/21/42/50−2/51/3−1/30]⏟TJ[2−11]⏟𝒙+[11/4−7/50]⏟𝒄J (3.248)
=[2−11]⏟𝒙. (3.249)

conforme esperado.

1import numpy as np
2import numpy.linalg as npla
3
4# sistema
5A = np.array([[-4., 2., -1.],
6 [-2., 5., 2.],
7 [1., -1., -3.]])
8b = np.array([-11., -7., 0.])
9
10# A = L + D + U
11L = np.tril(A, -1)
12D = np.diag(np.diag(A))
13U = np.triu(A, 1)
14
15# matriz de Jacobi
16T = -npla.inv(D) @ (L + U)
17# vetor de Jacobi
18c = npla.inv(D) @ b

A forma matricial das iterações de Jacobi consiste em

𝒙(0)=aprox. inicial, (3.250)
𝒙(k+1)=TJ⁢𝒙(k)+𝒄J, (3.251)

onde 𝒙(k)=(x1(k),x2(k),…,xn(k)) é a k-ésima aproximação da solução do sistema, k=0,1,2,….

Equivalentemente, tem-se a forma algébrica das iterações de Jacobi

xi(k+1)=bi−∑j≠ij=1nai⁢j⁢xj(k)ai⁢i, (3.252)

com i=1,2,…,n. Por não requerer as computações da matriz TJ e vetor 𝒄J, esta é a forma mais usada em implementações computacionais.

Exemplo 3.3.2.

Consideramos o sistema A⁢𝒙=𝒃 com

A=[−42−1−2521−1−3], (3.253)
𝒃=[−11−70]. (3.254)

Aplicando o método de Jacobi com aproximação inicial 𝒙(1)=(0,0,0) obtemos os resultados da Tabela 3.1.

k 𝒙(k) ‖A⁢𝒙(k)−𝒃‖
0 (0,0,0) 2.4⁢e+1
1 (2.75,−1.4,0) 7.4⁢e+0
2 (2.05,−0.3,1.38) 4.6⁢e+0
3 (2.25,−1.13,0.78) 2.2⁢e+0
4 (1.99,−0.81,1.13) 1.4⁢e+0
5 (2.06,−1.06,0.93) 6.9⁢e−1
⋮ ⋮ ⋮
29 (2,−1,1) 5.7⁢e−7
Tabela 3.1: Resultados referentes ao Exemplo 3.3.2.
Código 8: jacobi.py
1import numpy as np
2import numpy.linalg as npla
3
4def jacobi(A, b, x0, maxiter = 100,
5 tol=4.9e-8, atol=4.9e-8):
6
7 n = b.size
8 info = -1
9
10 x = np.empty_like(x0)
11 nres = npla.norm(A@x - b)
12 print(f'\n{0}: {x0}, {nres:.1e}')
13
14 # iterações
15 for k in range(maxiter):
16 for i in range(n):
17 x[i] = b[i]
18 # x[i] -= np.dot(A[i,:i], x0[:i])
19 for j in range(i):
20 x[i] -= A[i,j]*x0[j]
21 # x[i] -= np.dot(A[i,i+1:], x0[i+1:])
22 for j in range(i+1,n):
23 x[i] -= A[i,j]*x0[j]
24 x[i] /= A[i,i]
25 # critério de parada
26 nres = npla.norm(A@x - b)
27 print(f'{k+1}: {x}, {nres:.1e}')
28 if (nres <= max(tol*npla.norm(b), atol)):
29 info = 0
30 break
31 x0 = x.copy()
32
33 return x, info

A aplicação de Jacobi pode, então, ser feita como segue:

1# sistema
2A = np.array([[-4., 2., -1.],
3 [-2., 5., 2.],
4 [1., -1., -3.]])
5
6b = np.array([-11., -7., 0.])
7
8x0 = np.zeros_like(b)
9x, info = jacobi(A, b, x0)

3.3.2 Método de Gauss-Seidel

Como acima, começamos considerando um sistema linear A⁢𝒙=𝒃 e a decomposição A=L+D+U, onde L é a parte triangular inferior de A, D é sua parte diagonal e U sua parte triangular superior. Então, observamos que

A⁢𝒙=𝒃 (3.255)
(L+D+U)⁢𝒙=𝒃 (3.256)
(L+D)⁢𝒙=−U⁢𝒙+𝒃 (3.257)
𝒙=−(L+D)−1⁢U⁢𝒙+(L+D)−1⁢𝒃. (3.258)

Isto nos leva a forma matricial iteração de Gauss-Seidel

𝒙(1)=aprox. inicial, (3.259)
𝒙(k+1)=TG⁢𝒙(k)+𝒄G, (3.260)

onde

TG=−(L+D)−1⁢U, (3.261)
𝒄G=(L+D)−1⁢𝒃, (3.262)

são a matriz e o vetor de Gauss-Seidel.

Equivalentemente e mais adequada para implementação computacional, temos a forma algébrica da iteração de Gauss-Seidel

xi(k+1)=bi−∑j=1i−1ai⁢j⁢xj(k+1)−∑j=i+1nai⁢j⁢xj(k)ai⁢i, (3.263)

para i=1,2,…,n.

Exemplo 3.3.3.

Consideramos o sistema A⁢𝒙=𝒃 com

A=[−42−1−2521−1−3], (3.264)
𝒃=[−11−70]. (3.265)

Aplicando o método de Gauss-Seidel com aproximação inicial 𝒙(1)=(0,0,0) obtemos os resultados da Tabela 3.2.

Tabela 3.2: Resultados referentes ao Exemplo 3.3.3.
k 𝒙(k) ‖A⁢𝒙(k)−𝒃‖
0 (0,0,0) 2.4⁢e+1
1 (2.75,−0.3,1.02) 2.6⁢e+0
2 (2.35,−0.87,1.07) 1.2⁢e+0
3 (2.05,−1.01,1.02) 2.5⁢e−1
4 (1.99,−1.01,1) 4.0⁢e−2
5 (1.99,−1,1) 2.0⁢e−2
⋮ ⋮ ⋮
13 (2,−1,1) 2.3⁢e−7
Código 9: gs.py
1import numpy as np
2import numpy.linalg as npla
3
4def gs(A, b, x0, maxiter = 100,
5 tol=4.9e-8, atol=4.9e-8):
6
7 n = b.size
8 info = -1
9
10 x = np.empty_like(x0)
11 nres = npla.norm(A@x - b)
12 print(f'\n{0}: {x0}, {nres:.1e}')
13
14 # iterações
15 for k in range(maxiter):
16 for i in range(n):
17 x[i] = b[i]
18 # x[i] -= np.dot(A[i,:i], x[:i])
19 for j in range(i):
20 x[i] -= A[i,j]*x[j]
21 # x[i] -= np.dot(A[i,i+1:], x0[i+1:])
22 for j in range(i+1,n):
23 x[i] -= A[i,j]*x0[j]
24 x[i] /= A[i,i]
25 # critério de parada
26 nres = npla.norm(A@x - b)
27 print(f'{k+1}: {x}, {nres:.1e}')
28 if (nres <= max(tol*npla.norm(b), atol)):
29 info = 0
30 break
31 x0 = x.copy()
32
33 return x, info

3.3.3 Análise numérica

Observamos que ambos os métodos de Jacobi e de Gauss-Seidel consistem de iterações da forma

𝒙(k+1)=T⁢𝒙(k)+𝒄, (3.266)

para k=1,2,…, com x(0) uma aproximação inicial dada, T e c a matriz e o vetor de iteração, respectivamente. O seguinte teorema nos fornece uma condição suficiente e necessária para a convergência de tais métodos.

Teorema 3.3.1.

Para qualquer 𝒙(0)∈ℝn, temos que a sequência {𝒙(k)}k=0∞ dada por

𝒙(k+1)=T⁢𝒙(k)+𝒄, (3.267)

converge para a solução única de 𝒙=T⁢𝒙+𝒄 se, e somente se, ρ⁢(T)<1101010ρ⁢(T) é o raio espectral da matriz T, i.e. o máximo dos módulos dos autovalores de T..

Demonstração.

Veja [2, Cap. 7, Sec. 7.3]. ∎

Observação 3.3.1.(Taxa de convergência)

Para uma iteração da forma (3.266), temos a seguinte estimativa

‖𝒙(k)−𝒙‖≈ρ⁢(T)k+1⁢‖𝒙(0)−𝒙‖𝟏, (3.268)

onde 𝒙 é a solução de 𝒙=T⁢𝒙+𝒄.

Exemplo 3.3.4.

Consideramos o sistema A⁢𝒙=𝒃 com

A=[−42−1−2521−1−3], (3.269)
𝒃=[−11−70]. (3.270)

Nos Exemplo 3.3.2 e Exemplo 3.3.3 vimos que ambos os métodos de Jacobi e de Gauss-Seidel são convergentes para este sistema. Este convergiu aproximadamente duas vezes mais rápido que esse. Isto é confirmado pelos raios espectrais das respectivas matrizes de iteração

ρ⁢(TJ) ≈0.56, (3.271)
ρ⁢(TG) ≈0.26. (3.272)
1import numpy as np
2import numpy.linalg as npla
3
4# matriz dos coefs
5A = np.array([[-4., 2., -1.],
6 [-2., 5., 2.],
7 [1., -1., -3.]])
8
9# A = L + D + U
10L = np.tril(A, -1)
11D = np.diag(np.diag(A))
12U = np.triu(A, 1)
13
14# matriz de Jacobi
15TJ = -npla.inv(D) @ (L + U)
16rho_TJ = max(np.abs(npla.eigvals(TJ)))
17print(f'rho(T_J) = {rho_TJ:.2f}')
18
19# matriz de Gauss-Seidel
20TG = -npla.inv(L+D) @ U
21rho_TG = max(np.abs(npla.eigvals(TG)))
22print(f'rho(T_G) = {rho_TG:.2f}')
Observação 3.3.2.(Matriz estritamente diagonal dominante)

Pode-se mostrar que se A é uma matriz estritamente diagonal dominante, i.e. se

|ai⁢i|>∑j≠ij=1n|ai⁢j|, (3.273)

para todo i=1,2,…,n, então ambos os métodos de Jacobi e de Gauss-Seidel são convergentes.

3.3.4 Exercícios

E. 3.3.1.

Considere o seguinte sistema linear

−4⁢x1+x2+x3−x4=−1 (3.274)
5⁢x2−x3+2⁢x4=3 (3.275)
−x1+4⁢x3−2⁢x4=−2 (3.276)
x1−x2−5⁢x4=1 (3.277)

Compute a aproximação 𝒙(5) obtida da aplicação do Método de Jacobi com aproximação inicial 𝒙(0)=(1,1,−1,−1). Também, compute ‖A⁢𝒙(5)−b‖.


𝒙(5)=(0.325164,0.593096,−0.546128,−0.253043); ‖A⁢𝒙(5)−b‖=7.2⁢e−3

E. 3.3.2.

Considere o sistema linear do Exercício . Compute a aproximação 𝒙(5) obtida da aplicação do Método de Gauss-Seidel com aproximação inicial 𝒙(0)=(1,1,−1,−1). Também, compute ‖A⁢𝒙(5)−b‖.


𝒙(5)=(0.325208,0.592417,−0.545397,−0.253442); ‖A⁢𝒙(5)−b‖=7.1⁢e−4

E. 3.3.3.

Verifique que o Método de Jacobi, com 𝒙(0)=𝟎, é divergente para o sistema A⁢𝒙=𝒃 com

A=[−301−11−14002100−11−5], (3.278)
b=[−2−112−19] (3.279)

Então, escreva um sistema equivalente para o qual o Método de Jacobi seja convergente para qualquer escolha de aproximação inicial.


A=[−301−102101−1400−11−5], (3.280)
b=[−22−11−19] (3.281)
E. 3.3.4.

Considere o sistema linear

−x1+2⁢x2−2⁢x3=6 (3.282)
3⁢x1−4⁢x2+x3=−11 (3.283)
x1−5⁢x2+3⁢x3=−10. (3.284)

Empregando a aproximação inicial 𝒙(0)=𝟎, compute a solução com o:

  1. a)

    Método de Jacobi.

  2. b)

    Método de Gauss-Seidel.


a) divergente. b) (−2,1,−1)

E. 3.3.5.

Considere o sistema linear A⁢𝒙=𝒃 com

A=[10000−12−1000−12−1000−12−100001], (3.285)
b=[00.50.50.50]. (3.286)

Compute a solução empregando o:

  1. a)

    Método de Jacobi.

  2. b)

    Método de Gauss-Seidel.


𝒙=(0,0.75,1,0.75,0)

E. 3.3.6.

Considere o seguinte sistema de equações

2.1⁢x1−x2+0.9⁢x3=−1.3 (3.287)
2⁢x1+2⁢x2+2.1⁢x3=2.1 (3.288)
−1.2⁢x1−x2+2⁢x3=−π (3.289)

Usando a aproximação inicial 𝒙(0)=𝟎, verifique que o método de Jacobi não converge para sua solução, enquanto que o método de Gauss-Seidel converge. Por quê?


ρ⁢(TJ)=1.12⁢e+0>1, ρ⁢(TG)=5.48⁢e−1<1.

E. 3.3.7.

Considere o seguinte sistema de equações

1.1⁢x1+2⁢x2−1.9⁢x3=−1.3 (3.290)
x1+x2+x3=2.1 (3.291)
2.1⁢x1+1.9⁢x2+x3=−π (3.292)

Usando a aproximação inicial 𝒙(0)=𝟎, verifique que o método de Jacobi converge para sua solução, enquanto que o método de Gauss-Seidel não converge. Por quê?


ρ⁢(TJ)=8.5⁢e−1<1, ρ⁢(TG)=1.95⁢e+0>1.


Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen
| | | |