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

Cálculo Numérico - Versão Scilab

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.199) a21x1 + a22x2 + + a2nxn = y2 (4.200) (4.201) an1x1 + an2x2 + + annxn = yn (4.202)

Isolando o elemento x1 da primeira equação temos

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

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.204) x2(k+1) = y2 a21x1(k) + a 23x3(k) + + a 2nxn(k) a22 (4.205) (4.206) xn(k+1) = y2 an1x1(k) + + a n,n2xn2(k) + a n,n1xn1(k) ann (4.207)

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

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

usando o método de Jacobi iniciando com x(1) = y(1) = 0. x(k+1) = 23 y(k) 10 (4.212) y(k+1) = 26 x(k) 8 (4.213) x(2) = 23 y(1) 10 = 2,3 (4.214) y(2) = 26 x(1) 8 = 3,25 (4.215) x(3) = 23 y(2) 10 = 1,975 (4.216) y(3) = 26 x(2) 8 = 2,9625 (4.217)

Exemplo 4.7.2. Considere o seguinte sistema

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

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(k1)



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 Scilab explicativo. Escreva você mesmo o código.
Veja como participar da escrita do livro em:
https://github.com/livroscolaborativos/CalculoNumerico

Código Scilab: Método de Jacobi

function [x,deltax]=jacobi(A,b,x,tol,N)  
n=size(A,1)  
xnew     =x  
convergiu=%F                //FALSE;  
 
k=1  
while k<=N & ~convergiu  
  xnew(1)=(b(1) - A(1,2:n)*x(2:n))/A(1,1)  
  for i=2:n-1  
    xnew(i)=(b(i) -A(i,1:i-1)*x(1:i-1) ...  
                  -A(i,i+1:n)*x(i+1:n) )/A(i,i)  
  end  
  xnew(n)=  (b(n) -A(n,1:n-1)*x(1:n-1) )/A(n,n)  
 
  deltax=max( abs(x-xnew) )  
  if deltax<tol then  
     convergiu=%T       //TRUE  
  end  
  k=k+1  
  x=xnew                 // atualiza x  
  disp([k,x,deltax])    // depuracao  
end  
if ~convergiu then  
    error(Nao convergiu)  
end  
 
endfunction  

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

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

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

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

usando o método de Gauss-Seidel com condições iniciais x(1) = y(1) = 0. x(k+1) = 23 y(k) 10 (4.228) y(k+1) = 26 x(k+1) 8 (4.229) x(2) = 23 y(1) 10 = 2,3 (4.230) y(2) = 26 x(2) 8 = 2,9625 (4.231) x(3) = 23 y(2) 10 = 2,00375 (4.232) y(3) = 26 x(3) 8 = 2,9995312 (4.233)

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

function [x,deltax]=gauss_seidel(A,b,x,tol,N)  
n=size(A,1)  
xnew     =x  
convergiu=%F                //FALSE;  
 
k=1  
while k<=N & ~convergiu  
  xnew(1)=(b(1) - A(1,2:n)*x(2:n))/A(1,1)  
  for i=2:n-1  
    xnew(i)=(b(i) -A(i,1:i-1)*xnew(1:i-1) ...  
                  -A(i,i+1:n)*x(i+1:n) )/A(i,i)  
  end  
  xnew(n)=(b(n) -A(n,1:n-1)*xnew(1:n-1) )/A(n,n)  
 
  deltax=max( abs(x-xnew) )  
  if deltax<tol then  
     convergiu=%T       //TRUE  
  end  
  k=k+1  
  x=xnew                 // atualiza x  
  disp([k,x,deltax])    // depuracao  
end  
if ~convergiu then  
    error(Nao convergiu)  
end  
 
endfunction  

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

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

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

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

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

Exemplo 4.7.4. Considere o seguinte sistema linear: 3x1 + x2 x3 = 2 (4.239) x1 4x2 + x3 = 10 (4.240) x1 2x2 5x3 = 10 (4.241)

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

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

No Scilab, podemos construir as matrizes L, D e U, da seguinte forma:

-->A = [3 1 -1;-1 -4 1;1 -2 -5];  
-->D = eye(A).*A;  
-->L = tril(A)-D;  
-->U=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.244) Dx = (L + U)x + b (4.245) x = D1(L + U) =:TJx + D1b =:cJ. (4.246)

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

com x(1) uma aproximação inicial dada, sendo TJ := D1(L + U) a matriz de iteração e cJ = D1b 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 := D1(L+U) = 1 3 0 0 0 1 4 0 0 0 1 5 D1 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.248)

O vetor da iteração de Jacobi é:

cJ := D1b = 1 3 0 0 0 1 4 0 0 0 1 5 D1 2 10 10 b = 2 3 5 2 2 . (4.249)

No Scilab, podemos computar TJ e cJ da seguinte forma:

-->TJ = -inv(D)*(L+U);  
-->cJ = inv(D)*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.250) (L + D)x = Ux + b (4.251) x = (L + D)1U =:TGx + (L + D)1b =:cG (4.252)

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

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

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

No Scilab, podemos computar TG e cG da seguinte forma:

-->TG = -inv(L+D)*U;  
-->cG = inv(L+D)*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.256)

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

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

Assim, temos:

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

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

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

(I T)1 = k=0Tk. (4.261)

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

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

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(k1) + c) + c (4.264) = T2x(k1) + (I + T)c (4.265) (4.266) = T(k)x(1) + k=0k1Tk c. (4.267)

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

lim kx(k) = (I T)(1)c. (4.268)

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.269) = T(x x(k)) (4.270) (4.271) = T(k)(x x(1)) = T(k)y. (4.272)

Logo, para qualquer 1 p , temos, :

0 = lim kx x(k+1) = lim kT(k)y. (4.273)

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

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

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

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

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 dominante9 , 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.277)

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

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.279) x1 + 2x2 x3 = 1 (4.280) x2 + (2 + ε)x3 x4 = 1 (4.281) x3 + 2x4 x5 = 1 (4.282) x4 x5 = 1 (4.283)

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

Resposta.
epsilon=1e-3;  
 
A=[1 -1 0 0 0; -1 2 -1 0 0; 0 -1 (2+epsilon) -1 0; 0 0 -1 2 -1; 0 0 0 1 -1]  
 
v=[1 1 1 1 1]’  
xgauss=gauss([A v])  
 
function x=q_Jacobi()  
    x0=[0 0 0 0 0]’  
 
    i=0  
    controle=0  
    while controle<3 & i<1000  
    i=i+1  
 
    x(1)=1+x0(2)  
    x(2)=(1+x0(3)+x0(1))/2  
    x(3)=(1+x0(2)+x0(4))/(2+epsilon)  
    x(4)=(1+x0(3)+x0(5))/2  
    x(5)=x0(4)-1  
 
    delta=norm(x-x0,2)  
    if delta<1e-6 then  
        controle=controle+1  
    else  
        controle=0  
    end  
    mprintf(’i=%d, x1=%f, x5=%f, tol=%.12f\n’, i, x(1), x(5), delta)  
    x0=x;  
    end  
 
endfunction  
 
function x=q_Gauss_Seidel()  
    x0=[0 0 0 0 0]’  
 
    i=0  
    controle=0  
    while controle<3 & i<15000  
    i=i+1  
 
    x(1)=1+x0(2)  
    x(2)=(1+x0(3)+x(1))/2  
    x(3)=(1+x(2)+x0(4))/(2+epsilon)  
    x(4)=(1+x(3)+x0(5))/2  
    x(5)=x(4)-1  
 
    delta=norm(x-x0,2)  
    if delta<1e-2 then  
        controle=controle+1  
    else  
        controle=0  
    end  
    mprintf(’i=%d, x1=%f, x5=%f, tol=%.12f\n’,i,x(1),x(5),delta)  
    x0=x;  
    end  
 
endfunction

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

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

E 4.7.3. Refaça o Exercício 4.1.6 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 xj1 + 5xj xj+1 = cos(j10),2 j 10 x11 = x102 (4.285)

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 105. 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

Exemplos de rotinas:

function x=jacobi()  
    x0=zeros(11,1)  
    k=0;  
    controle=0;  
    while controle<3 & k<1000  
        k=k+1  
        x(1)=x0(2)  
        for j=2:10  
        x(j)=(cos(j/10)+x0(j-1)+x0(j+1))/5  
        end  
        x(11)=x0(10)/2  
 
 
        delta=norm(x-x0) //norma 2  
        if delta<1e-5 then  
            controle=controle+1  
        else  
            controle=0;  
        end  
        mprintf(’k=%d, x=[%f,%f,%f], tol=%.12f\n’,k,x(1),x(2),x(3),delta)  
        x0=x;  
    end  
 
 
endfunction  
 
function x=gs()  
    x0=zeros(11,1)  
    k=0;  
    controle=0;  
    while controle<3 & k<1000  
        k=k+1  
        x(1)=x0(2)  
        for j=2:10  
        x(j)=(cos(j/10)+x(j-1)+x0(j+1))/5  
        end  
        x(11)=x0(10)/2  
 
        delta=norm(x-x0) //norma 2  
        if delta<1e-5 then  
            controle=controle+1  
        else  
            controle=0;  
        end  
        mprintf(’k=%d, x=[%f,%f,%f], tol=%.12f\n’,k,x(1),x(2),x(3),delta)  
        x0=x;  
    end  
endfunction

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.286) 4x1 + x3 = 6 (4.287) 2x1 + x2 + 4x3 = 12 (4.288)

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 19/8/2020 às 17:36:32.

Informe erros ou edite você mesmo!