Use o botão abaixo para reportar erros ou dar sugestões.

Cálculo Numérico - Versão Python

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 a11x1 + a12x2 + + a1nxn = y1 (4.178) a21x1 + a22x2 + + a2nxn = y2 (4.179) (4.180) an1x1 + an2x2 + + annxn = yn (4.181)

Isolando o elemento x1 da primeira equação temos

x1(k+1) = y1 - a12x2(k) + + a 1nxn(k) a11 (4.182)

Note que utilizaremos os elementos xi(k) da iteração k (à 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.183) x2(k+1) = y2 - a21x1(k) + a 23x3(k) + + a 2nxn(k) a22 (4.184) (4.185) xn(k+1) = yn - an1x1(k) + + a n,n-2xn-2(k) + a n,n-1xn-1(k) ann (4.186)

Em notação mais compacta, o método de Jacobi consiste na iteração x(1) = aproximação inicial (4.187) xi(k+1) = y i - j=1 ji na ijxj(k) a ii (4.188)

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

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

Exemplo 4.7.2. Considere o seguinte sistema

- 3x1 + x2 + x3 = 2 2x1 + 5x2 + x3 = 5 2x1 + 3x2 + 7x3 = -17 (4.197)

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.


Aqui, cabe um código Python explicativo. Escreva você mesmo o código.
Veja como participar da escrita do livro em:
https://github.com/livroscolaborativos/CalculoNumerico

Código Python: Método de Jacobi

from __future__ import division  
import numpy as np  
from numpy import linalg  
 
def jacobi(A,b,x0,tol,N):  
    #preliminares  
    A = A.astype(double)  
    b = b.astype(double)  
    x0 = x0.astype(double)  
 
    n=np.shape(A)[0]  
    x = np.zeros(n)  
    it = 0  
    #iteracoes  
    while (it < N):  
        it = it+1  
        #iteracao de Jacobi  
        for i in np.arange(n):  
            x[i] = b[i]  
            for j in np.concatenate((np.arange(0,i),np.arange(i+1,n))):  
                x[i] -= A[i,j]*x0[j]  
            x[i] /= A[i,i]  
        #tolerancia  
        if (np.linalg.norm(x-x0,np.inf) < tol):  
            return x  
        #prepara nova iteracao  
        x0 = np.copy(x)  
    raise NameError(num. max. de iteracoes excedido.)

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 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.198)

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

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

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

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

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

from __future__ import division  
import numpy as np  
from numpy import linalg  
 
def gauss_seidel(A,b,x0,tol,N):  
    #preliminares  
    A = A.astype(double)  
    b = b.astype(double)  
    x0 = x0.astype(double)  
 
    n=np.shape(A)[0]  
    x = np.copy(x0)  
    it = 0  
    #iteracoes  
    while (it < N):  
        it = it+1  
        #iteracao de Jacobi  
        for i in np.arange(n):  
            x[i] = b[i]  
            for j in np.concatenate((np.arange(0,i),np.arange(i+1,n))):  
                x[i] -= A[i,j]*x[j]  
            x[i] /= A[i,i]  
            print(x[i],A[i,i])  
        #tolerancia  
        if (np.linalg.norm(x-x0,np.inf) < tol):  
            return x  
        #prepara nova iteracao  
        x0 = np.copy(x)  
    raise NameError(num. max. de iteracoes excedido.)

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 Ax = b, onde A = [ai,j]i,j=1n,n é a matriz (real) dos coeficientes, b = (aj)j=1n é o 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.213)

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,k 1, (4.214)

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.215)

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.216)

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.217)

Exemplo 4.7.4. Considere o seguinte sistema linear: 3x1 + x2 - x3 = 2 (4.218) -x1 - 4x2 + x3 = -10 (4.219) x1 - 2x2 - 5x3 = 10 (4.220)

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.221)

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.222)

Em Python, podemos construir as matrizes L, D e U, da seguinte forma:

>>> A = np.array([[3,1,-1],  
...               [-1,-4,1],  
...               [1,-2,5]],  
...               dtype=’double’)  
>>> D = np.diag(np.diag(A))  
>>> L = np.tril(A)-D  
>>> U=np.triu(A)-D

Iteração de Jacobi

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.223) Dx = -(L + U)x + b (4.224) x = -D-1(L + U) =:TJx + D-1b =:cJ. (4.225)

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.226)

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.227)

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.228)

Em python, podemos computar TJ e cJ da seguinte forma:

>>> TJ = -np.linalg.inv(D).dot(L+U);  
>>> cJ = np.linalg.inv(D).dot(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: Ax = b (L + D + U)x = b (4.229) (L + D)x = -Ux + b (4.230) x = -(L + D)-1U =:TGx + (L + D)-1b =:cG (4.231)

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.232)

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.233)

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.234)

Em Python, podemos computar TG e cG da seguinte forma:

-->TG = -np.linalg.inv(L+D).dot(U);  
-->cG = np.linalg.inv(L+D).dot(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 Ax = b e uma iteração:

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

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):

ρ(T) Tp, (4.236)

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.237)

Assim, temos:

lim kTk p lim k Tpm = 0. (4.238)

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

Observação 4.7.1. Observamos que:

lim kTk p = 0,,1 p , lim ktijk = 0,1 i,j n. (4.239)

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

(I - T)-1 = k=0Tk. (4.240)

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=0Tk = lim m(I -T) k=0mTk = lim m(I -Tm+1) = I, (4.241)

o que mostra que (I - T)-1 = k=0Tk.

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

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

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.243) = T2x(k-1) + (I + T)c (4.244) (4.245) = T(k)x(1) + k=0k-1Tk c. (4.246)

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

lim kx(k) = (I - T)(-1)c. (4.247)

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.248) = T(x* - x(k)) (4.249) (4.250) = T(k)(x* - x(1)) = T(k)y. (4.251)

Logo, para qualquer 1 p , temos, :

0 = lim kx* - x(k+1) = lim kT(k)y. (4.252)

Como y é arbitrário, da Observação 4.7.1 temos lim kT(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)kx(1) - x*. (4.253)

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.

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:

aii > j=1 ji n a ij ,i = 1,...,n (4.254)

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

aii j=1 ji n a ij ,i = 1,...,n (4.255)

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 dominante7 , 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:

x - x(k) x . (4.256)

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.257)

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 x1 - x2 = 1 (4.258) -x1 + 2x2 - x3 = 1 (4.259) -x2 + (2 + ε)x3 - x4 = 1 (4.260) -x3 + 2x4 - x5 = 1 (4.261) x4 - x5 = 1 (4.262)

  • Escreva na forma Ax = b e resolva usando eliminação gaussiana para ε = 10-3 no Python.
  • 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.263)

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(j10),2 j 10 x11 = x102 (4.264)

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.265) 4x1 + x3 = 6 (4.266) 2x1 + x2 + 4x3 = 12 (4.267)

Resposta.

Permute as linhas 1 e 2.

Creative Commons License 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 15/5/2019 às 15:24:50.

Informe erros ou edite você mesmo!