Métodos Iterativos de Jacobi e Gauss-Seidel - Resolução de Sistemas Lineares

O programa abaixo foi desenvolvido no compilador Dev-Pascal. Cuidado com os critérios de convergência! Pois esses métodos nem sempre convergem, gerando erros de "Overflow" quando utilizados para resolução de sistemas que não convergem.

Download do Código em Pascal e Executável do Programa:  
http://www.4shared.com/file/BuIURAFH/Iterativo_Jacobi_Gauss_Seidel.html

Download do Compilador Dev-Pascal: http://www.bloodshed.net/devpascal.html

Código:

Program iterativo_Jacobi_Gauss_Seidel;
//Programa para resolucao de sistemas de equacoes lineares por metodos iterativos de Jacobi e Gauss-Seidel//
Uses crt;
 Var
 a : Array[1..100,1..100] Of Real;
 x, x0, b, r, er, xg : Array[1..100] Of Real;
 i, j, n, m, k, ki, met, teste, convergencia : Integer;  //sao variaveis auxiliares//
 resposta : char;
 nome : String;
 erro, e, soma : Real;
 resultado: Text;

Begin
Repeat //Repeat utilizado para reiniciar o programa//
 repeat
 Clrscr;
 Textcolor (White);
 Writeln ('Resolucao de sistemas de equacoes lineares por metodos iterativos');
 Writeln;
 Write ('Digite a ordem do sistema (n), maximo 100: ');
 Readln (n);
 until n >= 2;

 //le os coeficientes da matriz A//
 for j:=1 to n do
 begin
 Repeat
 clrscr;
 Writeln ('AX = B');
 Writeln;
 Writeln ('Matriz ampliada do sistema');
 Writeln;
 Writeln ('Digite os elementos da coluna ',j,' da matriz A (coeficientes de x',j,'):');
 Writeln;
 for i:=1 to n do
 begin
 Write('a[',i,'x',j,']= ');
 readln(a[i,j]);
 end;
 writeln;
 writeln;
  Textcolor (White);
  Write ('Confira os valores digitados, estao corretos? (digite s=sim, n=nao): ');
  Readln (resposta);
 Until (resposta = 'S') Or (resposta = 's');
 end;

 //le os coeficientes da matriz B//
 repeat
 clrscr;
 Writeln ('AX = B');
 Writeln;
 Writeln ('Matriz ampliada do sistema');
 Writeln;
 Writeln ('Digite os elementos da matriz B (termos independentes)');
 Writeln;
 for i:=1 to n do
 begin
 Write('b[',i,']= ');
 readln(b[i]);
 end;
 writeln;
 writeln;
 Write ('Confira os valores digitados, estao corretos? (digite s=sim, n=nao): ');
 Readln (resposta);
 Until (resposta = 'S') Or (resposta = 's');

 Clrscr;
 teste:= 0;    //Zera as variaveis envolvidas no teste de convergencia//
 convergencia:= 0;

 For i:= 1 To n Do //Testa se existe elementos a[i,i] iguais a zero//
 Begin                
   If a[i,i] = 0 Then
    Begin
    Writeln('Erro:');
    Writeln;
    Writeln ('existem elementos da diagonal principal, a[i,i], iguais a zero (i=1,2...n).');
    Writeln;
    Writeln('a[i,i] diferente de zero para todo i e uma condicao necessaria');
    Writeln('para a aplicacao dos metodos iterativos de Jacobi e Gauss-Seidel.');
    Writeln('Reagrupe as equacoes para que essa condicao seja satisfeita.');
    Writeln;
    teste:= 1;
    Break;
    End;
 End;

  Soma:= 0;

 For i:= 1 To n Do //Testa o criterio de convergencia das linhas//
 Begin
   For j:= 1 To n Do
   soma := soma + ABS(a[i,j]);
 
   soma := soma - ABS(a[i,i]);
 
   If soma >= ABS(a[i,i]) Then
   Begin
     convergencia:= 1;
     Break;
   End;
   Soma:= 0;
 End;

  Soma:= 0;

 For j:= 1 To n Do //Testa o criterio de convergencia das colunas//
 Begin
   For i:= 1 To n Do
   soma := soma + ABS(a[i,j]);
 
   soma := soma - ABS(a[j,j]);
 
   If soma >= ABS(a[j,j]) Then
   Begin
     convergencia:= convergencia + 1;
     Break;
   End;
   Soma:= 0;
 End;

 If convergencia = 2 Then
 Begin
 Writeln ('O uso dos criterios das linhas e colunas detectou que talvez neo haja');
 Writeln ('convergencia. Se isso ocorrer tente utilizar metodos diretos, como o');
 Writeln ('de Gauss e Jordan com Retrosubstituicao.');
 Writeln;
  If teste <> 1 Then
   Begin
    Write ('Deseja continuar mesmo assim? (digite s=sim, n=nao): ');
    Readln (resposta);
    Clrscr;
     If (resposta = 'S') Or (resposta = 's') Then
       convergencia := 0;
   End;
 End;

 If (teste <> 1) AND (convergencia <> 2) Then //so faz os calculos se houver convergencia, se nao houver leva ate a pergunta para reinicializacao do programa//
                                             //convergencia = 0 ou = 1 significa que ira convergir, convergencia = 2 significa que nao ira//
 Begin                                      //teste = 1 significa que nao passou no teste da diagonal principal livre de a[i,i] = 0//

 Write ('Digite o valor do erro aceitavel: ');
 Readln (erro);
 Writeln;
 Write ('Digite o valor maximo de iteracoes (ate 999, digite 0 para ilimitada): ');
 Readln (m);
 Writeln;
 Write ('Gostaria de digitar uma solucao inicial? (digite s=sim, n=nao): ');
 Readln (resposta);

 If (resposta = 'S') Or (resposta = 's') Then
 Begin
 Writeln;
 Writeln ('Digite a solucao:');
 Writeln;
 For ki:= 1 To n Do
  Begin
   Write ('X',ki,' = ');
   Readln (x0[ki]);
  End;
 End;
 Writeln;

  Repeat
  Writeln ('Digite o numero do metodo iterativo a ser utilizado:');
  Writeln;
  Writeln ('1) Metodo de Jacobi');
  Writeln ('2) Metodo de Gauss-Seidel');
  Writeln;
  Write ('Metodo: ');
  Readln (met);
  Clrscr;
  Until (met=1) Or (met=2);

 Textcolor (White);


 i:= 0;
 j := 0;  //Zera para evitar valores anteriores//
 k:=0;

 Repeat

  //Iteracao Jacobi//

  If met=1 Then
  Begin
   For i:= 1 To n Do
   Begin
    For j:= 1 To n Do
     x[i] := x[i] + (-1* a[i,j] * x0[j]);
                                      
    x[i] := x[i] + (a[i,i]*x0[i]) + b[i];
    x[i] := x[i]/a[i,i];
   End;
  End;

  //Iteracao Gauss-Seidel//

  If met=2 Then
  Begin
  For ki:=1 To n Do
    xg[ki]:= x0[ki];
   For i:= 1 To n Do
   Begin
    For j:= 1 To n Do
     x[i] := x[i] + (-1* a[i,j] * xg[j]);
                                      
    x[i] := x[i] + (a[i,i]*xg[i]) + b[i];
    x[i] := x[i]/a[i,i];
    xg[i]:= x[i]
   End;
  End;
 

  //Testa o erro//

  teste:= 0;
  For ki:=1 To n Do
  Begin
   e:= Abs(x[ki] - x0[ki]); //Calcula o erro//
   er[ki] := e; //Guarda o erro para ser impresso//
   If e <= Abs(erro) Then
   teste := teste + 1; //Ou seja: uma condicao de parada e quando todos os erros forem menores ou iguais ao
  End;                //digitado pelo usuario. Isso se da quando teste = n (teste ganha +1 para cada erro aceitavel)
                     //n e o numero de variaveis do sistema, a sua ordem, logo serao calculados n erros, um para cada variavel.


  For ki:=1 To n Do  //Guarda os resultados da solucao do sistema//
  r[ki] := x[ki];

 //Prepara x0 e x para uma nova iteracao//

   For ki:=1 To n Do
    x0[ki]:= x[ki];
   For ki:=1 To n Do
    x[ki]:= 0;


  k:= k+1; //Conta o numeros de iteracoes//

   Textcolor (White);
   Writeln ('Resultados (com ', k, ' iteracoes):'); //Imprime os resultados na tela//
   Textcolor (Yellow);
   For ki:=1 To n Do                 
    Begin
    Writeln('X',ki, ' = ', r[ki], '  ','Erro = ',er[ki]);
    Writeln;
    End;

 Until (teste = n) Or (k = m) Or (k = 1000); //teste=n significa todos os erros aceitaveis, k=m que o maximo de iteracoes
                                            //que o usuario escolheu foi atingido e k=1000 uma parada de seguranca para evitar a nao convergencia

 Textcolor (White);

 If k = 1000 Then
 Writeln ('Sistema nao convergiu apos 1000 iteracoes!');
 Writeln;

 For ki:=1 To n Do   //Zera x0 para uma possivel reinicializacao do programa, evitando erros na proxima iteracao//
 x0[ki]:= 0;

                 
                 //Grava o resultado em um arquivo de texto//
 
  If k < 1000 Then
  Begin
   Write ('Deseja salvar os resultados em um arquivo de texto? (digite s=sim, n=nao): '); 
   Readln (resposta);
   If (resposta = 's') Or (resposta = 'S') Then
   Begin
   Writeln;
   Write ('Digite o nome do arquivo a ser salvo: ');
   Readln  (nome);
   nome := (nome + '.txt');
   Assign(resultado, nome);
   Rewrite (resultado);
   Writeln (resultado, 'valores de x (incognitas):');
   Writeln (resultado);
   //escreve os valores de x//
   For ki:=1 To n Do                 
   Writeln (resultado,'X',ki, ' = ', r[ki], '  ','Erro = ',er[ki]);
  
    //escreve os valores de b original e os calculados com as variaveis x encontradas//
     For i:= 1 To n Do
    Begin
     Writeln (resultado);
     Writeln (resultado, 'b[',i,'] original = ', b[i]);
     b[i]:=0;
     For j:= 1 To n Do
      Begin
      b[i]:= b[i] + (a[i,j]*r[j]); //Calcula os termos da matriz B com as variaveis encontradas//
     End;
     Writeln (resultado, 'b[',i,'] calculado = ', b[i]);
    End;
     Close (resultado);
     Writeln;
      Writeln ('Arquivo salvo!');
   End;
  End;

 Writeln;

 End; //End usado na condicao "If (teste <> 1) AND (convergencia <> 2)", que faz com que o programa pule os calculos se nao houver convergencia detectada//

 Write ('Deseja reiniciar o programa? (digite s=sim, n=nao): ');   //Reinicia o programa//
 Readln (resposta);


Until (resposta <> 's') And (resposta <> 'S');
End.


Se desejarem fazer qualquer melhoria, correção de algum eventual erro, fiquem a vontade! ; )