Resolvendo um sistema de equações lineares iterativamente

Em algumas situações, é mais interessante resolver um sistema de equações lineares Ax = b usando algum método iterativo do que um método direto, especialmente para o caso de sistemas de grande porte e com muitos elementos nulos. O método iterativo mais usual é Gauss-Seidel que pode ser expresso por:
\[ x_i = \frac{b_i - \sum_{j=1,j\neq i}^{N}a_{ij}x_{j}}{a_{ii}} \]
com $i = 1, 2, ..., N$. Em geral, a solução inicial é o vetor nulo: x = 0.

Como é um método iterativo, devemos adotar algum critério de parada. Os critérios de parada mais comuns são: número máximo de iterações e o erro menor que uma certa tolerância. O erro pode ser calculado como a menor diferença entre $x_i^k$ e $x_i^{k-1}$, onde $k$ significa o número da iteração.

Um outro ponto que devemos considerar é que a ordem das equações no sistema é importante. Podemo trocar linhas ou colunas (sem afetar o resultado efetivo do sistema) de tal forma que a diagonal principal se torne mais dominante, isto é, os elementos da diagonal principal devem ser, em termos absolutos, maiores que os coeficientes da mesma linha.

Usaremos o sistema $8 \times 8$ abaixo para ilustrar esses conceitos:

Esse sistema é não singular nem mesmo mal condicionado, mas está longe de ser o tipo de sistema ideal para resolvermos de forma iterativa. O código abaixo calcula a sua solução:

clc; close; // Limpando o console e fechando janelas

M = [2.4 2.2 1.3 1.2 1.2 1.1 0.5 0.3
     -1 3.3 -1 1 1 0.5 0.25 0.2
     -1 -1 2.5 -1.4 -1.7 -1  0.1 0.2
      1  0.5 2.0  2.6 1  0.5 0.2 0.1 
     0.5 0.5 -1.2 -1 2.1 1 -0.2 0.1
     0.5 1.5 -1 1 -1 2 -0.1 -0.2
     0.25 -0.5 1.5 -1.3 1 -1 2.1 -0.5
     -0.3 0.6 0.5 1.5 -1.2 1 -1.3 -2]
b = M*[1;2;3;4;5;6;-0.5;-2]; // vetor 
N = 8; // ordem da matriz

Mi = M;
xk = 0*b;  // solução inicial
bi = b;
for k=1:N 
    Mi(k,k) = 0;  Mi(k,:) = Mi(k,:)/M(k,k);
    bi(k) = bi(k)/M(k,k);
end

nmax = 90; // número máximo de iterações
erro = 1; // erro inicial
nk = 2;  // número da iteração
mvx = zeros(N,nmax); // guardando os valores de xk
mvx(:,1) = xk;
while erro>0.005
    erro = xk;
    for n=1:N
        s = Mi(n,:)*xk;
        xk(n) = bi(n) - s; 
    end    
    mvx(:,nk) = xk;
    erro = max(abs(erro - xk));
    printf('%d & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f %1.4f \n',...
    nk,xk(1),xk(2),xk(3),xk(4),xk(5),erro);
    nk = nk + 1;
    if nk > nmax then 
        erro =0; 
    end; // número máximo de iterações
end
figure;
plot(1:nk-1,mvx(:,1:nk-1),'-o'); xgrid;
xlabel('Iteração k'); ylabel('$ x_i \$');
legend('x1','x2','x3','x4','x5','x6','x7','x8');

O resultado podemos ver no gráfico abaixo:

A convergência só ocorre com quase 90 iterações. Muito demorado!

Percebemos que o total de iterações para ocorrer a convergência é muito longo, pois são necessárias quase 90 iterações. O que podemos fazer para melhorar isso? Trocar linhas/colunas e iniciar com um ponto de partida melhor. Um ponto de partida possivelmente melhor é fazer $x_i = b_i/a_{ii}$, aplicando esse ganho obtemos uma redução de algumas iterações, é um ganho pequeno. Um ganho mais significativo é se trocarmos as incógnitas $x_1$ com $x_2$, isso implica em trocar a 2a. com a 3a. linha e, depois, a segunda coluna com a primeira coluna. Nesse caso, o ganho é muito significativo, o número de iterações cai para apenas 19, mantendo o mesmo critério de erro. Se, adicionalmente, somarmos a 3a com a 4a. coluna, fazemos com a diagonal fique mais dominante, acelerando ainda mais a convergência. O código com essas alterações é:

clc; close; // Limpando o console e fechando janelas

M = [2.4 2.2 1.3 1.2 1.2 1.1 0.5 0.3
     -1 3.3 -1 1 1 0.5 0.25 0.2
     -1 -1 2.5 -1.4 -1.7 -1  0.1 0.2
      1  0.5 2.0  2.6 1  0.5 0.2 0.1 
     0.5 0.5 -1.2 -1 2.1 1 -0.2 0.1
     0.5 1.5 -1 1 -1 2 -0.1 -0.2
     0.25 -0.5 1.5 -1.3 1 -1 2.1 -0.5
     -0.3 0.6 0.5 1.5 -1.2 1 -1.3 -2]
b = M*[1;2;3;4;5;6;-0.5;-2]; // vetor 
N = 8; // ordem da matriz

// troca ...
aux = M(1,:); ba = b(1);
M(1,:) = M(2,:); b(1) = b(2);
M(2,:) = aux; b(2) = ba;
aux = M(:,1); M(:,1) = M(:,2); M(:,2) = aux;

//Soma 3a. linha com 4a. linha, resultado na 3a. linha:
M(3,:) = M(3,:) + M(4,:); b(3) = b(3) + b(4);
 
Mi = M;
xk = 0*b;  // solução inicial
bi = b;
for k=1:N 
  xk(k) = sign((b(k)/M(k,k))); //*sqrt(abs((b(k)/M(k,k))));
end
for k=1:N 
    Mi(k,k) = 0;  Mi(k,:) = Mi(k,:)/M(k,k);
    bi(k) = bi(k)/M(k,k);
end

nmax = 90; // número máximo de iterações
erro = 1; // erro inicial
nk = 2;  // número da iteração
mvx = zeros(N,nmax); // guardando os valores de xk
mvx(:,1) = xk;
while erro>0.005
    erro = xk;
    for n=1:N
        s = Mi(n,:)*xk;
        xk(n) = bi(n) - s; 
    end    
    mvx(:,nk) = xk;
    erro = max(abs(erro - xk));
    printf('%d & %1.4f & %1.4f & %1.4f & %1.4f & %1.4f %1.4f \n',...
    nk,xk(1),xk(2),xk(3),xk(4),xk(5),erro);
    nk = nk + 1;
    if nk > nmax then 
        erro =0; 
    end; // número máximo de iterações
end
figure;
plot(1:nk-1,mvx(:,1:nk-1),'-o'); xgrid;
xlabel('Iteração k'); ylabel('$ x_i \$');
legend('x1','x2','x3','x4','x5','x6','x7','x8');

O resultado é:

O mesmo sistema, mas com alguns ajustes "técnicos" - redução significativa no número de iterações.