calc_num_09_sistemas_lineares_iterativos

Prévia do material em texto

Instituto Federal Minas Gerais - Campus Bambuí
Departamento de Engenharia e Computação - DEC
Curso de Engenharia de Computação - ENGCOMP
Cálculo Numérico
09 – Sistemas Lineares - Métodos Iterativos
Marcos Roberto Ribeiro
2021
Introdução
Os métodos diretos para resolução de sistemas lineares são adequados para sistemas de
porte pequeno (até 30 equações) a médio (até 50 equações)
Para sistemas maiores, os sistemas diretos começam a enfrentar problemas de desempenho
Nessas situações, é mais conveniente utilizar métodos iterativos
Os métodos iterativos iniciam com um aproximação inicial (que pode ser ruim) e vão
melhorando essa aproximação através de sucessivas iterações
IFMG - Campus Bambuí - DEC - ENGCOMP 2/26
Método de Jacobi I
Considere um sistema linear no seguinte formato:
a1,1x1 + a1,2x2 + ...+ a1,nxn = b1
a2,1x1 + a2,2x2 + ...+ a2,nxn = b2
...
am,1x1 + am,2x2 + ...+ am,nxn = bm
Isolando o x1 da primeira equação temos:
x (k+1)1 =
b1 − (a1,2x (k)2 + ...+ a1,nx
(k)
n )
a1,1
Da mesma forma, isolando o elemento xi de cada equação i , para todo i ∈ {2, ..., n}
podemos construir um processo iterativo
IFMG - Campus Bambuí - DEC - ENGCOMP 3/26
Método de Jacobi II
Basicamente, temos
x (1) = aproximação inicial
x (k+1)i =
bi − n∑
j=1
j 6=i
ai ,jx (k)j
 /ai ,i
IFMG - Campus Bambuí - DEC - ENGCOMP 4/26
Exemplo com o método de Jacobi
Como exemplo, vamos resolver o seguinte sistema em Python:{
10x0 + x1 = 23
x0 + 8x1 = 26
Vamos considerar x (0)0 = x
(0)
1 = 0
IFMG - Campus Bambuí - DEC - ENGCOMP 5/26
Resolvendo o exemplo em Python
Código
1 import numpy as np
2
3 A = np.array([[10, 1],[1, 8]],
dtype='double')↪→
4 b = np.array([23, 26], dtype='double')
5 x = np.array([0, 0], dtype='double')
6 xant = np.array([0, 0], dtype='double')
7
8 print('k:', 'x(k)')
9 for k in range(10):
10 print(k, ':', x)
11 x[0] = (b[0] - A[0,1] * xant[1])/A[0,0]
12 x[1] = (b[1] - A[1,0] * xant[0])/A[1,1]
13 xant = np.copy(x)
Resultado
k: x(k)
0 : [ 0. 0.]
1 : [ 2.3 3.25]
2 : [ 1.975 2.9625]
3 : [ 2.00375 3.003125]
4 : [ 1.9996875 2.99953125]
5 : [ 2.00004687 3.00003906]
6 : [ 1.99999609 2.99999414]
7 : [ 2.00000059 3.00000049]
8 : [ 1.99999995 2.99999993]
9 : [ 2.00000001 3.00000001]
IFMG - Campus Bambuí - DEC - ENGCOMP 6/26
Critério de parada
No caso dos métodos iterativos, podemos definir um número máximo de iterações ou um
erro menor do que uma tolerância como critérios de parada
Calculamos o erro usando norma da diferença entre a solução da iteração atual e da iteração
anterior:
||x (k) − x (k−1)||∞
Agora vamos ver uma implementação completa em Python para resolver os sistema
−3x0 + x1 + x2 = 2
2x0 + 5x1 + x2 = 5
2x0 + 3x1 + 7x2 = −17
Vamos considerar x (0) = (1, 1,−1), máximo de 20 iterações e tolerância de 10−4
IFMG - Campus Bambuí - DEC - ENGCOMP 7/26
Método de Jacobi em Python I
1 #!/usr/bin/python3
2 # -*- coding: utf-8 -*-
3
4 import numpy as np
5
6
7 def jacobi(A, b, x0, tol, iteracoes):
8 n = np.alen(A) # n é o número de linhas
9 x = np.zeros(n, dtype='double') # Solução atual
10 xant = x0 # Solução da iteração anterior
11 for k in range(iteracoes): # Iterações do método
12 for i in range(n): # Iterações para cada incógnita
13 x[i] = b[i] # Termo constante
14 for j in range(i): # Incógnitas anteriores
15 x[i] -= A[i,j]*xant[j]
16 for j in range(i+1, n): # Incógnitas posteriores
17 x[i] -= A[i,j]*xant[j]
IFMG - Campus Bambuí - DEC - ENGCOMP 8/26
Método de Jacobi em Python II
18 x[i] /= A[i,i] # Divisão pelo coeficiente da incógnita atual
19 erro = np.linalg.norm(x-xant, np.inf) # Cálculo do erro
20 print("Iteração {k:3d}: ".format(k=k+1) +
21 "x={x}, ".format(x=np.round(x,8)) +
22 "Erro={e:+5.8f}".format(e=erro))
23 if (erro < tol): # Testa se erro é menor que a tolerância
24 return x
25 # Copia solução atual para ser a anterior na próxima iteração
26 xant = np.copy(x)
27
28
29 # Entrada: M, b
30 M = np.array([
31 [ -3, 1, 1],
32 [ 2, 5, 1],
33 [ 3, 3, 7]], dtype='double')
34 b = np.array([2, 5, -17], dtype='double')
IFMG - Campus Bambuí - DEC - ENGCOMP 9/26
Método de Jacobi em Python III
35 # Aproximação inicial
36 x0 = np.array([1, 1, -1], dtype='double')
37 # Executa o método de Jacobi
38 x = jacobi(M, b, x0, 0.0001, 20)
39 print('\nSolução aproximada encontrada')
40 print('x =', x)
IFMG - Campus Bambuí - DEC - ENGCOMP 10/26
Método de Jacobi em Python – Execução
Iteração 1: x=[-0.66666667 0.8 -3.28571429], Erro=+2.28571429
Iteração 2: x=[-1.4952381 1.92380952 -2.48571429], Erro=+1.12380952
Iteração 3: x=[-0.85396825 2.0952381 -2.6122449 ], Erro=+0.64126984
Iteração 4: x=[-0.83900227 1.86403628 -2.96054422], Erro=+0.34829932
Iteração 5: x=[-1.03216931 1.92770975 -2.86787172], Erro=+0.19316704
Iteração 6: x=[-0.98005399 1.98644207 -2.81237447], Erro=+0.05873232
Iteração 7: x=[-0.94197747 1.95449649 -2.85988061], Erro=+0.04750613
Iteração 8: x=[-0.96846137 1.94876711 -2.86250815], Erro=+0.02648390
Iteração 9: x=[-0.97124701 1.95988618 -2.84870246], Erro=+0.01380569
Iteração 10: x=[-0.96293876 1.9582393 -2.85227393], Erro=+0.00830825
Iteração 11: x=[-0.96467821 1.95563029 -2.8551288 ], Erro=+0.00285487
Iteração 12: x=[-0.9664995 1.95689704 -2.85326518], Erro=+0.00186362
Iteração 13: x=[-0.96545604 1.95725284 -2.85302752], Erro=+0.00104346
Iteração 14: x=[-0.96525823 1.95678792 -2.8536272 ], Erro=+0.00059968
Iteração 15: x=[-0.96561309 1.95682873 -2.85351273], Erro=+0.00035487
Iteração 16: x=[-0.96556133 1.95694778 -2.85337813], Erro=+0.00013460
Iteração 17: x=[-0.96547678 1.95690016 -2.85345134], Erro=+0.00008455
Solução aproximada encontrada
x = [-0.96547678 1.95690016 -2.85345134]
IFMG - Campus Bambuí - DEC - ENGCOMP 11/26
Método de Gauss-Seidel
Assim como o método de Jacobi, o método de Gauss-Seidel também isola o elemento xi da
equação i
Porém, podemos notar que x (k)i depende dos demais elementos da iteração k − 1
Intuitivamente, podemos pensar em usar x (k−1)j com j < k para calcular x
(k)
i
Basicamente, o método de Gauss-Seidel consiste em:
x (0) = aproximação inicial
x (k)i =
bi − i−1∑
j=1
ai ,jx (k)j
n∑
j=i+1
ai ,jx (k−1)j
 /ai ,i
IFMG - Campus Bambuí - DEC - ENGCOMP 12/26
Exemplo com o método Gauss-Seidel
Como exemplo, vamos considerar novamente o sistema:{
10x0 + x1 = 23
x0 + 8x1 = 26
Vamos considerar x (0)0 = x
(0)
1 = 0
Buscaremos por uma solução aproximada em Python com método Gauss-Seidel
IFMG - Campus Bambuí - DEC - ENGCOMP 13/26
Resolvendo o exemplo em Python
Código
1 import numpy as np
2
3 A = np.array([[10, 1],[1, 8]],
dtype='double')↪→
4 b = np.array([23, 26], dtype='double')
5 x = np.array([0, 0], dtype='double')
6 xant = np.array([0, 0], dtype='double')
7
8 print('k:', 'x(k)')
9 for k in range(10):
10 print(k, ':', x)
11 x[0] = (b[0] - A[0,1] * xant[1])/A[0,0]
12 x[1] = (b[1] - A[1,0] * x[0])/A[1,1]
13 xant = np.copy(x)
Resultado
k: x(k)
0 : [0. 0.]
1 : [2.3 2.9625]
2 : [2.00375 2.99953125]
3 : [2.00004687 2.99999414]
4 : [2.00000059 2.99999993]
5 : [2.00000001 3. ]
6 : [2. 3.]
7 : [2. 3.]
8 : [2. 3.]
9 : [2. 3.]
IFMG - Campus Bambuí - DEC - ENGCOMP 14/26
Método Gauss-Seidel em Python
Agora vamos ver uma implementação completa em Python para resolver o sistema:
−3x0 + x1 + x2 = 2
2x0 + 5x1 + x2 = 5
2x0 + 3x1 + 7x2 = −17
Vamos considerar x (0) = (1, 1,−1), máximo de 20 iterações e tolerância de 10−4
IFMG - Campus Bambuí - DEC - ENGCOMP 15/26
Método de Gauss-Seidel em Python I
1 #!/usr/bin/python3
2 # -*- coding: utf-8 -*-
3
4 import numpy as np
5
6
7 def gauss_seidel(A, b, x0, tol, iteracoes):
8 n = np.alen(A) # n é o número de linhas
9 x = np.zeros(n, dtype='double') # Solução atual
10 xant = x0 # Solução da iteração anterior
11 for k in range(iteracoes): # Iterações do método
12 for i in range(n): # Iterações para cada incógnita
13 x[i] = b[i] # Termo constante
14 for j in range(i): # Incógnitas anteriores
15 x[i] -= A[i,j]*x[j]
16 for j in range(i+1, n): # Incógnitas posteriores
17 x[i] -= A[i,j]*xant[j]
IFMG - Campus Bambuí - DEC - ENGCOMP16/26
Método de Gauss-Seidel em Python II
18 x[i] /= A[i,i] # Divisão pelo coeficiente da incógnita atual
19 erro = np.linalg.norm(x-xant, np.inf) # Cálculo do erro
20 print("Iteração {k:3d}: ".format(k=k+1) +
21 "x={x}, ".format(x=np.round(x,8)) +
22 "Erro={e:+5.8f}".format(e=erro))
23 if (erro < tol): # Testa se erro é menor que a tolerância
24 return x
25 # Copia solução atual para ser a anterior na próxima iteração
26 xant = np.copy(x)
27
28
29 # Entrada: M, b
30 M = np.array([
31 [ -3, 1, 1],
32 [ 2, 5, 1],
33 [ 3, 3, 7]], dtype='double')
34 b = np.array([2, 5, -17], dtype='double')
IFMG - Campus Bambuí - DEC - ENGCOMP 17/26
Método de Gauss-Seidel em Python III
35 # Aproximação inicial
36 x0 = np.array([1, 1, -1], dtype='double')
37 # Executa o método de Gauss-Seidel
38 x = gauss_seidel(M, b, x0, 0.0001, 20)
39 print('\nSolução aproximada encontrada')
40 print('x =', x)
IFMG - Campus Bambuí - DEC - ENGCOMP 18/26
Método de Gauss-Seidel em Python - Execução
Iteração 1: x=[-0.66666667 1.46666667 -2.77142857], Erro=+1.77142857
Iteração 2: x=[-1.1015873 1.99492063 -2.81142857], Erro=+0.52825397
Iteração 3: x=[-0.93883598 1.93782011 -2.85670748], Erro=+0.16275132
Iteração 4: x=[-0.97296246 1.96052648 -2.85181315], Erro=+0.03412648
Iteração 5: x=[-0.96376222 1.95586752 -2.85375941], Erro=+0.00920024
Iteração 6: x=[-0.96596396 1.95713747 -2.85336007], Erro=+0.00220174
Iteração 7: x=[-0.96540753 1.95683503 -2.85346893], Erro=+0.00055643
Iteração 8: x=[-0.96554463 1.95691164 -2.853443 ], Erro=+0.00013710
Iteração 9: x=[-0.96551045 1.95689278 -2.85344957], Erro=+0.00003418
Solução aproximada encontrada
x = [-0.96551045 1.95689278 -2.85344957]
IFMG - Campus Bambuí - DEC - ENGCOMP 19/26
Decomposição da Matriz A
Os métodos iterativos podem ser construídos como uma iteração de ponto fixo
No caso dos sistemas lineares, podemos reescrever da seguinte maneira:
Ax = b ⇔ x = Tx + c
T é a matriz de iteração e c o vetor de iteração
Assim, o método iterativo consiste em computar a iteração:
x (k+1) = Tx (k) + c, para k ≥ 1
A matriz A de um sistema Ax = b, pode ser decomposta como A = L + D + U, onde onde
D é a matriz diagonal, L e U são as matrizes triangular inferior e superior
a1,1 a1,2 ... a1,n
a2,1 a2,2 ... a2,n
...
...
. . .
...
an,1 an,2 ... an,n

︸ ︷︷ ︸
A
=

0 0 ... 0
a2,1 0 ... 0
...
...
. . .
...
an,1 an,2 ... 0

︸ ︷︷ ︸
L
+

a1,1 0 ... 0
0 a2,2 ... 0
...
...
. . .
...
0 0 ... an,n

︸ ︷︷ ︸
D
+

0 a1,2 ... a1,n
0 0 ... a2,n
...
...
. . .
...
0 0 ... 0

︸ ︷︷ ︸
U
IFMG - Campus Bambuí - DEC - ENGCOMP 20/26
Exemplo de decomposição

3x1 + x2 − x3 = 2
−x1 − 4x2 + x3 = −10
x1 − 2x2 − 5x3 = 10
A =
 3 1 −1−1 −4 1
1 −2 −5
 , x =
x1x2
x3
 e b =
 2−10
10

 3 1 −1−1 −4 1
1 −2 −5

︸ ︷︷ ︸
A
=
 0 0 0−1 0 0
1 −2 0

︸ ︷︷ ︸
L
+
3 0 00 −4 0
0 0 −5

︸ ︷︷ ︸
D
+
0 1 −10 0 1
0 0 0

︸ ︷︷ ︸
U
IFMG - Campus Bambuí - DEC - ENGCOMP 21/26
Decomposição em Python
1 ...
2 A = np.array([[3,1,-1],
3 [-1,-4,1],
4 [1,-2,5]],
5 dtype='double')
6 D = np.diag(np.diag(A))
7 L = np.tril(A)-D
8 U=np.triu(A)-D
IFMG - Campus Bambuí - DEC - ENGCOMP 22/26
Iteração de Jacobi
Ax = b ⇔ (L + D + U)x = b
⇔ Dx = −(L + U)x + b
⇔ x = −D−1(L + U)︸ ︷︷ ︸
TJ
x + D−1b︸ ︷︷ ︸
cJ
.
Assim, a iteração do método de Jacobi é x (k+1) = TJx (k) + cJ , com x (1) como aproximação
inicial, TJ = −D−1(L + U) a matriz de iteração e cJ = D−1b o vetor da iteração
TJ e cJ em Python
1 ...
2 TJ = -np.linalg.inv(D).dot(L+U);
3 cJ = np.linalg.inv(D).dot(b);
IFMG - Campus Bambuí - DEC - ENGCOMP 23/26
Iteração de Gauss-Seidel
Ax = b ⇔ (L + D + U)x = b
⇔ (L + D)x = −Ux + b
⇔ x = −(L + D)−1U︸ ︷︷ ︸
TG
x + (L + D)−1b︸ ︷︷ ︸
cG
A iteração do método de Gauss-Seidel é x (k+1) = TGx (k) + cG , com x (1) como aproximação
inicial, TG = −(L + D)−1U a matriz de iteração e cG = (L + D)−1b o vetor da iteração
TG e cG em Python
1 ...
2 TG = -np.linalg.inv(L+D).dot(U);
3 cG = np.linalg.inv(L+D).dot(b);
IFMG - Campus Bambuí - DEC - ENGCOMP 24/26
Condições de convergência
A garantia de convergência dos métodos iterativos pode ser verificada pelo critério do raio
espectral
Seja ρ(M) o raio espectral de M (maior autovalor absoluto de M), os métodos interativos
convergem se o raio espectral da matriz de iteração (T ) for menor ou igual a 1, ρ(T ) ≤ 1
Cálculo do Raio Espectral em Python
1 ...
2 av, _ = np.linalg.eig(T) # Autovalores de T
3 raio_espectral = max(abs(av))
IFMG - Campus Bambuí - DEC - ENGCOMP 25/26
Referências
FRANCO, N. B. Cálculo Numérico. São Paulo: Pearson Prentice Hall, 2006. p. 505.
JUSTO, D. A. R. et al. Cálculo Numérico: um livro colaborativo versão python. Porto Alegre,
2020. Disponível em: <https://www.ufrgs.br/reamat>.
IFMG - Campus Bambuí - DEC - ENGCOMP 26/26
https://www.ufrgs.br/reamat
	Introdução
	Método de Jacobi
	Método de Gauss-Seidel
	Convergência dos métodos iterativos
	Referências